为什么使用 FFT 在信号中舍入频率值？

Question

所以，我试图弄清楚如何在实践中使用 DFT 来检测信号中的普遍频率。我一直试图围绕傅立叶变换是什么以及 DFT 算法如何工作，但显然我还有很长的路要走。我编写了一些代码来生成信号（因为目的是处理音乐，所以我生成了一个主要的 C 和弦，因此是奇怪的频率值），然后尝试回到频率数。这是我的代码sr = 44100 # sample ratex = np.linspace(0, 1, sr) # one second of signaltpi = 2 * np.pidata = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)freqs = np.fft.fftfreq(sr)fft = np.fft.fft(data)idx = np.argsort(np.abs(fft))fft = fft[idx]freqs = freqs[idx]print(freqs[-6:] * sr)这给了我[-262.  262. -330.  330. -392.  392.]不同于我编码的频率（261.63、329.63 和 392.0）。我做错了什么，我该如何解决？

Answer 1

实际上，如果帧持续T数秒，则 DFT 的频率是k/THz，其中 k 是整数。因此，过采样不会提高估计频率的准确性，只要这些频率被识别为 DFT 幅度的最大值。相反，考虑到持续 100 秒的较长帧会导致 DFT 频率之间的间距为 0.01Hz，这可能足以产生预期的频率。通过将峰值频率估计为其相对于功率密度的平均频率，可能会更好。

图 1：即使在应用了 Tuckey 窗口之后，加窗信号的 DFT 也不是 Dirac 的总和：在峰值底部仍然存在一些频谱泄漏。在估计频率时必须考虑到该功率。

另一个问题是帧的长度不是信号周期的倍数，它可能不是周期性的。尽管如此，DFT 的计算方式好像信号是周期性的，但在帧的边缘是不连续的。它会引起被描述为频谱泄漏的杂散频率。加窗是处理此类问题和缓解与人为不连续相关问题的参考方法。实际上，窗口的值在靠近帧的边缘处不断减小到零。 有一个窗口函数列表，在scipy.signal 中有很多窗口函数可用。一个窗口被应用为：

tuckey_window=signal.tukey(len(data),0.5,True)

data=data*tuckey_window

在这一点上，出现最大幅度的频率仍然是 262、330 和 392。应用窗口只会使峰值更加明显：加窗信号的 DFT 具有三个不同的峰值，每个都具有中心瓣和旁瓣，具体取决于窗口的 DFT。这些窗口的波瓣是对称的：中心频率因此可以计算为峰值的平均频率，与功率密度有关。

import numpy as np

from scipy import signal

import scipy

sr = 44100 # sample rate

x = np.linspace(0, 1, sr) # one second of signal

tpi = 2 * np.pi

data = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)

#a window...

tuckey_window=signal.tukey(len(data),0.5,True)

data=data*tuckey_window

data -= np.mean(data)

fft = np.fft.rfft(data, norm="ortho")

def abs2(x):

        return x.real**2 + x.imag**2

fftmag=abs2(fft)[:1000]

peaks, _= signal.find_peaks(fftmag, height=np.max(fftmag)*0.1)

print "potential frequencies ", peaks

#compute the mean frequency of the peak with respect to power density

powerpeak=np.zeros(len(peaks))

powerpeaktimefrequency=np.zeros(len(peaks))

for i in range(1000):

    dist=1000

    jnear=0

    for j in range(len(peaks)):

        if dist>np.abs(i-peaks[j]):

             dist=np.abs(i-peaks[j])

             jnear=j

    powerpeak[jnear]+=fftmag[i]

    powerpeaktimefrequency[jnear]+=fftmag[i]*i

powerpeaktimefrequency=np.divide(powerpeaktimefrequency,powerpeak)

print 'corrected frequencies', powerpeaktimefrequency

由此产生的估计频率为 261.6359 Hz、329.637Hz 和 392.0088 Hz：它比 262、330 和 392Hz 好得多，并且满足这种纯无噪声输入信号所需的 0.01Hz 精度。

Answer 2

DFT 结果箱在频率上由 Fs/N 分隔，其中 N 是 FFT 的长度。因此，DFT 窗口的持续时间限制了 DFT 结果箱频率中心间距方面的分辨率。

但是，对于低噪声（高 S/N）中分离良好的频率峰值，您可以通过在 DFT 结果箱之间插入 DFT 结果来估计频率峰值位置，以达到更高的分辨率，而不是增加数据的持续时间。您可以尝试对粗略频率峰值位置估计进行抛物线插值，但加窗 Sinc 插值（本质上是香农-惠特克重建）将提供更好的频率估计精度和分辨率（假设感兴趣的频率峰值周围有足够低的本底噪声，例如，在您的人工波形情况下没有附近的正弦曲线）。

Answer 3

由于您希望获得 0.01 Hz 的分辨率，因此您需要对至少 100 秒的数据进行采样。您将能够解析高达约 22.05 kHz 的频率。