使用NumPy FFT提取相位信息

Question

我正在尝试使用快速傅里叶变换来提取单个正弦函数的相位移。我知道在纸上，如果我们将函数的变换表示为T，则有以下关系：

然而，我发现虽然我能够准确捕捉余弦波的频率，但除非我以极高的速率进行采样，否则相位是不准确的。例如：

import numpy as np
import pylab as pl

num_t = 100000
t = np.linspace(0,1,num_t)
dt = 1.0/num_t
w = 2.0*np.pi*30.0
phase = np.pi/2.0

amp = np.fft.rfft(np.cos(w*t+phase))
freqs = np.fft.rfftfreq(t.shape[-1],dt)

print (np.arctan2(amp.imag,amp.real))[30]

pl.subplot(211)
pl.plot(freqs[:60],np.sqrt(amp.real**2+amp.imag**2)[:60])
pl.subplot(212)
pl.plot(freqs[:60],(np.arctan2(amp.imag,amp.real))[:60])
pl.show()

使用num = 100000个点，我得到了1.57173880459的相位。

使用num = 10000个点，我得到了1.58022110476的相位。

使用num = 1000个点，我得到了1.6650441064的相位。

出了什么问题？即使有1000个点，每个周期有33个点，这应该足以解决。也许有一种方法可以增加计算频率点的数量吗？是否有任何方法可以在“低”点数下完成此操作？

编辑：通过进一步实验，似乎我需要每个周期约1000个点才能准确提取相位。为什么？！

编辑2：进一步实验表明，精度与每个周期的点数有关，而不是绝对数字。增加每个周期采样点的数量会使相位更准确，但如果信号频率和采样点数都按相同因素增加，则精度保持不变。

- KBriggs

---

如果您增加点的数量会发生什么？1,000,000个？10,000,000个？相位是否渐近地接近pi/2？ - wallyk

我不是专家，但随着采样率的降低，FFT频率和相位结果中不可避免地会出现越来越多的噪声。较低的采样率会导致频率分量幅度的扩散和主分量相位的偏移。你还会看到结果中主频率幅度的减小，因为其他分量增加了。 - DisappointedByUnaccountableMod

是的，随着点数增加，相位会越来越准确。我应该补充一下，重要的不是绝对点数，而是每个周期的点数，这使得它更准确。在1kHz下采样30Hz信号10秒钟与在1kHz下采样1kHz信号1秒钟具有相同的准确性，但在1kHz下采样3Hz波形1秒钟与在10kHz下采样30Hz波形1秒钟具有相同的准确性。 - KBriggs

当您将量化添加到例如8或12位的样本中时会发生什么？ - DisappointedByUnaccountableMod

我不确定我理解了。你所说的“量化...位”是什么意思？ - KBriggs

Answer 1

你的点并没有均匀地分布在区间上，你在末尾重复了一个点：0 等同于 1。当你取得的点越多时，这个问题就不那么重要了，但仍然会产生一些误差。你可以完全避免这种情况，linspace 有一个标志可以解决这个问题。此外，它还有一个标志可以直接返回数组中的 dt。

做

t, dt = np.linspace(0, 1, num_t, endpoint=False, retstep=True)

替代

t = np.linspace(0,1,num_t)
dt = 1.0/num_t

然后它可以正常工作 :)

- Ilja

---

啊，所以我的dt错了，差了1/N的因素？看起来好像是可以工作的...这有点尴尬 - 显然需要阅读一下linspace的文档。谢谢！ - KBriggs

不，你的dt并不是问题，这只是为了方便而已。问题在于，在结束那个时期和开始下一个时期时，你考虑了同一个点两次。 - Ilja

但是我如何获得准确的真实世界数据呢？如果我使用实时时间序列，我不能保证最后一个点不对应第一个点，对吧？ - KBriggs

Answer 2

在未旋转的FFT结果中，结果区间内的相位值仅在输入信号在FFT长度内完全为整数周期时才正确。你的测试信号不是这样的，因此FFT测量到的是测试正弦波端点间信号不连续性的部分相关内容。更高的采样率将创建与正弦波略有不同的最后一个端点，从而可能产生较小的不连续性。

如果要减少此FFT相位测量误差，请将测试信号创建为参考测试向量的确切中心（样本N/2），而不是第一个样本的测试相位，并进行fftshift操作（旋转N/2），以便在长度为N的结果FFT输入向量中的第一个和最后一个点之间没有信号不连续性。

- hotpaw2

---

你能给一个简短的例子吗？在我有实时数据而且不能控制周期性的真实世界例子中，我如何提取有意义的相位值？ - KBriggs

对于任意信号，只需将FFT相位结果参考到FFT向量的中心，而不是任何端点。您可以通过使用fftshift预处理步骤或后处理来实现此目的：取决于bin编号是奇数还是偶数，您可能需要翻转相位的符号（这是另一个域中的fftshift）。 - hotpaw2

在“现实世界”中，您无法控制FFT孔径的周期性，因此通常会对数据进行窗口处理，并且绝对希望在隆起（窗口函数的中心）附近测量相位。由于常见的窗口函数通常在那里接近零，因此没有信号可以测量端点处的相位。 - hotpaw2

您能否推荐一本相关的参考书籍来讨论这个问题？我不是很理解如何在窗口中间测量相位，或者如何将相位参考到FFT向量的中心？ - KBriggs

好的，没问题。等我有时间把它恰当地表达出来后，我会稍后发布的。谢谢。 - KBriggs

Answer 3

这段代码可能会有所帮助：

def reconstruct_ifft(data):
"""
In this function, we take in a signal, find its fft, retain the dominant modes and reconstruct the signal from that

Parameters
----------
data : Signal to do the fft, ifft

Returns 
-------
reconstructed_signal : the reconstructed signal

"""
N = data.size
yf = rfft(data)

amp_yf = np.abs(yf) #amplitude

yf = yf*(amp_yf>(THRESHOLD*np.amax(amp_yf)))

reconstructed_signal = irfft(yf)

return reconstructed_signal

0.01是您想要保留的fft振幅的阈值。将THRESHOLD设置得更大（大于1没有任何意义）会导致较少的模式和更高的均方根误差，但确保更高的频率选择性。 （请调整Python代码的TAB）

- RAHUL MALLIK