1

我有一个简单的正弦函数作为 sin(2*pi f t+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。

4

4 回答 4

3

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

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

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

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

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

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

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

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

于 2015-04-15T15:50:39.867 回答
3

这实际上比最初看起来更难回答。@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功能，但对于大多数目的，上面的增量就足够了）

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

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)fftshift在 fft 之前执行（是的，在时域中）或 (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])因此您可以直接使用它。对于实际信号，可能会有一些偏差，因此插值增加了一些鲁棒性，这通常是一个好主意，尽管使用wandwp和归一化可能会过大，angle(sgn*(X[i0]-X[i0+1))通常就足够了。

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

于 2019-05-13T21:18:01.550 回答
1

XX当您应该从 FFT ( )中获取相位时，您正试图从功率谱 ( ) 中获取相位X。改变：

phase_Estimate=angle(XX(ind));

至：

phase_Estimate=angle(X(ind));
于 2015-04-15T10:23:08.050 回答
0

可能晚了，但我稍微改变了你的剧本

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 计算度数而不是弧度，并且可以轻松更改输入相移。

于 2022-02-04T21:03:29.833 回答