在NumPy中使用FFT时的频率单位

Question

我在NumPy中使用FFT函数进行一些信号处理.我有一个数组signal ,每小时有一个数据点,共有576个数据点.我使用以下代码signal来查看其傅里叶变换.

t = len(signal)
ft = fft(signal,n=t)
mgft=abs(ft)
plot(mgft[0:t/2+1])

我看到两个峰值,但我不确定x轴的单位是什么,即它们如何映射到小时？任何帮助,将不胜感激.

Answer 1

给定采样率FSample和变换块大小N,您可以使用以下关系计算频率分辨率deltaF,采样间隔deltaT和总捕获时间capT:

deltaT = 1/FSample = capT/N
deltaF = 1/capT = FSample/N

还要记住的是,FFT从返回值0至FSample,或等价-FSample/2于FSample/2.在你的情节,你已经跌落-FSample/2到0分手了.NumPy包含一个辅助函数来为您计算所有这些:fftfreq.

为了您的价值观deltaT = 1 hour和N = 576,你deltaF = 0.001736 cycles/hour = 0.04167 cycles/day,从-0.5 cycles/hour到0.5 cycles/hour.因此,如果您在bin 48(和bin 528)处具有幅度峰值,则对应于频率分量at48*deltaF = 0.0833 cycles/hour = 2 cycles/day.

通常,在计算FFT之前,应该对时域数据应用窗口函数,以减少频谱泄漏.汉恩窗口几乎从来都不是一个糟糕的选择.您还可以使用该rfft函数跳过-FSample/2, 0输出的一部分.那么,你的代码将是:

ft = np.fft.rfft(signal*np.hanning(len(signal)))
mgft = abs(ft)
xVals = np.fft.fftfreq(len(signal), d=1.0) # in hours, or d=1.0/24 in days
plot(xVals[:len(mgft)], mgft)