首先演示 $$\begin{align}&[\dots, 0, 0, 1,\hphantom{-}1, 0, 0, \dots] \text{ 和}\\ &[\点, 0, 0, 1, -1, 0, 0, \dots]\end{align}$$

$$\begin{align}&[\dots, 0, 0, 1,\hphantom{-}1, 0, 0, \dots] \text{ and}\\ &[\dots, 0, 0, 1, -1, 0, 0, \dots]\end{align}$$

平等的

$$[\dots, 0, 0, 1, \hphantom{-}1, 0, 0, \dots]\hphantom{\text{ and}}$$

但它们的 sinc 插值的平方不同（图 1）：

图 1. $[1, 1]$（蓝色）和 $[1, -1]$（红色）的 sinc 插值平方。$[1, 1]$ (blue) and $[1, -1]$ (red).

这表明通常不可能从在带限信号的临界采样频率处获取的均匀样本中恢复带限信号的平方。

---

让我们在 Octave 中测试不同的插值方法。带限插值的黄金标准是 DFT–zero-pad–DFT：

>> format free
>> x = [1 2 3];
>> y = [1 2 3  0 0 0];
>> abs(fft(x))
ans =

 6 1.73205 1.73205

>> abs(fft(y))
ans =

 6 4.3589 1.73205 2 1.73205 4.3589

最后一组数字是从完美插值的频域 bin 计算的幅度值。让我们尝试对幅度进行插值：

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)))), [0])))
ans =

 6 4.57735 1.73205 0.309401 1.73205 4.57735

好像有点离谱。现在让我们尝试插值平方幅度：

>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)).^2)), [0])))
ans =

 36 25 3 -8 3 25

>> sqrt(ans)
ans =

 (6,0) (5,0) (1.73205,0) (0,2.82843) (1.73205,0) (5,0)

它不仅关闭，而且还sqrt()返回一个负内插值的复数。那么插入 bin 值的唯一有效方法是什么？

让我们再给一次这个机会，尝试对上采样 2 倍的频域数据进行插值。由于偶数长度变换，这需要复制“奈奎斯特样本”，如果代码变得难以阅读，请道歉。

>> z = [1 2 3  0 0 0  0 0 0  0 0 0];
>> abs(fft(z))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

以上就是我们想要的。让我们尝试插值2x 上采样幅度：

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)))), [0 0 0]))
ans =

 3.36365 1.10447 0.318175 0 0 0 0 0 0 -0.208949 0.318175 1.10447

>> ans(4) = ans(length(ans)-2)
ans =

 3.36365 1.10447 0.318175 -0.208949 0 0 0 0 0 -0.208949 0.318175 1.10447

>> fft(ans)
ans =

 5.79105 5.59483 4.56785 2.7273 1.5231 1.76882
 2.20895 1.76882 1.5231 2.7273 4.56785 5.59483

那仍然关闭。让我们尝试对2x 上采样平方幅度进行插值：

>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)).^2)), [0 0 0]))
ans =

 14 8 3 0 0 0 0 0 0 1.18424e-15 3 8

>> ans(4) = ans(length(ans)-2)
ans =

 14 8 3 1.18424e-15 0 0 0 0 0 1.18424e-15 3 8

>> sqrt(fft(ans))
ans =

 6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485

现在它完美无缺！带回家的信息是在尝试在频域中插值之前至少将（时域零填充）上采样两倍，并插值平方幅度而不是幅度。它之所以有效，是因为取平方幅度与将每个 bin 值乘以其复共轭相同。复共轭保留了数据表示的带限函数的带宽，所以乘法使“时域带宽”翻倍，因为它相当于时域卷积。请注意，在选择插值方法时，2x 上采样平方幅度仍然是严格采样的，因此进一步的过采样应该会使精确插值更容易。

我大约回答了我自己的问题，但会接受进一步的见解作为答案！

PS刚刚发现还有interpft，它用较少的语法进行插值。

利用附加信息

使用有关数据的附加信息，插值变得更容易甚至更精确，例如，它是一个严格采样的平方时移 sinc。在这种情况下，给定两个样本 $\alpha$ 就在最大样本之前和 $\gamma$ 就在最大样本之后，峰值的时间 $-1<d<1$ 可以通过以下方式计算：$\alpha$ just before and $\gamma$ just after the largest sample, the time $-1<d<1$ of the peak can be calculated by:

$$d = \begin{cases}0&\text{if }\alpha = \gamma,\\ \frac{1 - \sqrt{1 - p^2}}{p},\, p = \frac{γ - α}{α + γ}&\text{otherwise,}\end{cases}$$

$d = 0$ 意味着 sinc 正好移动到最大样本的时间。可以在半采样率下进行相同的插值，这是基础函数的临界采样频率，其幅度等于时移 sinc 的幅度，两个连续的最大值样本 $\alpha$ 和 $\beta $ 基础函数绝对值的平方：$d = 0$ meaning that the sinc is shifted to exactly the time of the largest sample. The same interpolation can be done at half sample rate, which is the critical sampling frequency of the underlying function the magnitude of which equals that of a time-shifted sinc, with the two successive largest-valued samples $\alpha$ and $\beta$ of the square of the absolute value of the underlying function:

$$0 \le d \le 1\\ d = \begin{cases}0&\text{if }\beta = 0,\\ 1&\text{if }\alpha = 0,\\ \frac{1 + p - \sqrt{1 - p^2}}{2p},\, p = \frac{\beta - α}{α + \beta}&\text{otherwise,}\end{cases}$$

公式不受数据幅度缩放的影响。对于频率估计，您很少有纯粹的实时时移 sinc，但如果有，则有精确的插值公式。

可以使用峰值区域周围的几个样本与预先计算的插值向量的点积来计算 DFT 的插值点。插值向量由所需插值样本的位置确定，同时考虑到所需的零填充量等。

此技术和计算插值向量的方法包含在本文档的附录 B 中：

http://ericjacobsen.org/FTinterp.pdf

我希望这会有所帮助。