从 FFT 恢复正弦信号的相位

Question

我有一个简单的正弦函数 sin(2*pift+phi)。我想获得相位信号phi。 我尝试使用 FFT 来计算 phi。在 matlab 中，我执行以下操作

f=200; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
phase = 3/5*pi; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=sin(2*pi*f*t+phase); %replace with cos if a cosine wave is desired

NFFT=1024; %NFFT-point DFT
X=fft(x,NFFT); %compute DFT using FFT
XX=2*abs(X(1:NFFT/2+1));
[tt ind]=max(XX);
phase_Estimate=angle(X(ind);

这个结果对我来说几乎没有意义。例如当phi=0.523时，得到phase_Estimate-0.98。

Answer 1

仅当正弦曲线的周期恰好是 FFT 长度的整数约数时，使用非插值 FFT 结果相位才有效。在您的示例中，正弦波在孔径中不是整数周期。

如果不是，您将需要对相位进行插值以获得更好的估计。这是获得更好插值相位的一种方法：

在执行 FFT 之前，首先 fftshift（旋转 N/2）数据以将零相位参考点移动到窗口的中心。 （这是为了防止相邻 FFT 结果箱之间的相位翻转/交替。*）

然后进行 FFT 并通过抛物线或更好的 Sinc 插值估计正弦曲线的频率。

然后使用估计的频率对最近的两个 FFT 结果 bin 相位之间的相位进行线性插值。更新：或者更好的是，分别使用 FFT 结果的实部和虚部的 Sinc 插值，然后在插值的 IQ 分量上使用 atan2 以获得插值相位。

然后使用窗口中心的估计频率和相位来计算其他点的相位，例如 FFT 窗口的开始。

还要注意，正弦波的相位与余弦波的相位相差 pi/2。 atan(im,re) 返回余弦相位。

(* 作为预先 fftshit-ing 数据的一种替代方法，也可以对奇数 FFT 结果箱的相位进行后翻转。)

Answer 2

这个问题实际上比最初看起来更难回答。 @hotpaw2 给出的答案是完全正确的，并且比我找到的任何其他资源都更好地拼写出来，但它仍然只是一个大纲，我花了几个小时才把所有的肉都放在骨头上。

希望其他人也能找到相关的问题（并供我自己将来参考），这里有一个更彻底的解释：

假设您在索引 ind 处有一个（本地）最大值（如问题的情况）。

第 1 步： 尝试使用两个周围的值插入更精确的最大值位置。这在很多很多地方都有很好的解释，例如 https://www.dsprelated.com/freebooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html 对如何做到这一点有很好的解释，但 TL:DR 版本是：

delta = 0.5*(X[ind-1]-X[ind+1])/(X[ind-1]-2*X[ind]+X[ind+1])
p0 = ind+delta

估计峰值位于p0 （如果你想要更精确的估计，请改用log(X[ind-1])，或者全力以赴使用sinc函数，但对于大多数用途，上面的delta就足够了）

第 2 步：棘手的部分：使用该位置插入相位。 第一个直觉是使用我们刚刚找到的 delta 进行简单的线性插值：

i0 = floor(p0); w = p0-i0; wp = 1-w
ang = wp*angle(X[i0]) + w*angle(X[i0+1])

这将不起作用有多种原因，@hotpaw2 概述了其中的大部分。第一个是这不是你平均角度的方式，因为它们以 2pi 为模，所以 0 和 2pi 应该是相似的。更正确的方法是取归一化复数的平均值：

ang = angle(wp*X[i0]/abs(X[i0]) + w*X[i0+1]/abs(X[i0+1]))

但是，这仍然不正确，因为如果峰值位于 i0 和 i0+1 之间，则相位会在此处翻转 180 度（pi 弧度），从而使平均值非常具有误导性。要修复此“相位翻转”，您必须 (a) 在 fft 之前执行 fftshift（是的，在时域中）或 (b) 翻转 X 的每个奇数索引值的相位（通过将复数相乘来实现带有-1的数字）或（如果您不愿意像我一样触摸FFT），您也可以（c）使用以下代码模拟方法（b）：

i0 = floor(p0); w = p0-i0; wp = 1-w
if (i0 % 2 == 1) { w*=-1; wp*=-1 } # Flip both if i0 odd
ang = angle(wp*X[i0]/abs(X[i0]) - w*X[i0+1]/abs(X[i0+1])) # Note the "-" here!

这将为您提供（大部分）正确的相位，但对于余弦并且位于 fft 窗口的中心。

第 3 步（可选）：如果您需要相位为正弦并且从窗口的开头开始，您需要添加一个校正因子：

ang_beg = ang - (2*pi*p0/N)*N/2 + pi*0.5 = ang - pi*(p0 - 0.5)

（0.5*pi 将 cos 转换为 sin，-p0*pi 转换为窗口的开头）。 这似乎有效，至少在我需要它的 Phase Vocoder 中。希望其他人也会发现这很有用。

顺便说一句，纯正弦波不需要相位插值，如angle(X[i0]) = angle(-X[i0+1])，因此您可以直接使用它。对于实际信号，可能存在一些偏差，因此插值增加了一些鲁棒性，这通常是一个好主意，尽管使用 w 和 wp 和归一化可能有点过分，而 angle(sgn*(X[i0]-X[i0+1)) 通常就足够了。

非常欢迎任何有关这一切的 cmets。我不是 DSP 专家，所以我可能在某些细节上错了，但这似乎确实有效，所以希望其他人也会发现它有用。

Answer 3

您试图从功率谱 (XX) 中获取相位，而您应该从 FFT (X) 中获取相位。变化：

phase_Estimate=angle(XX(ind));

到：

phase_Estimate=angle(X(ind));

Answer 4

可能晚了，但我稍微修改了你的脚本

f=200; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
shift = 30
phase = shift*pi/180; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave
t=0:1/fs:nCyl*1/f; %time base
x=cos(2*pi*f*t+phase); %replace with cos if a cosine wave is desired

NFFT=4096; %NFFT-point DFT
X=fft(x,NFFT); %compute DFT using FFT
XX=2*abs(X(1:NFFT/2+1));
[tt, ind]=max(XX);
phase_Estimate = angle(X(ind)) * 360/(2*pi)

它得出的结果与我的预期非常接近。

我将 x 向量生成更改为余弦，以 phase_Estimate 而非弧度计算度数，并且可以轻松更改输入相移。