错误结果绘图窗口 FFT

Posted 2023-02-25

【中文标题】错误结果绘图窗口 FFT

【英文标题】：Bad result plotting windowing FFT

【发布时间】：2016-08-05 04:40:05

【问题描述】：

我正在使用 python 和 scipy 来理解窗口化，我做了一个图来查看窗口在 FFT 下的行为，但结果不是我所期望的。 情节是：

中间的情节是纯 FFT 情节，这是我得到奇怪东西的地方。 然后我改变了三角。获取泄漏的函数，将 1 直接用于数组的 300 个第一项，结果：

代码：

sign_freq=80
sample_freq=3000
num=np.linspace(0,1,num=sample_freq)
i=0
#wave data:
sin=np.sin(2*pi*num*sign_freq)+np.sin(2*pi*num*sign_freq*2)
while i<1000:
    sin[i]=1
    i=i+1


#wave fft:
fft_sin=np.fft.fft(sin) 
fft_freq_axis=np.fft.fftfreq(len(num),d=1/sample_freq)

#wave Linear Spectrum (Rms)
lin_spec=sqrt(2)*np.abs(np.fft.rfft(sin))/len(num)
lin_spec_freq_axis=np.fft.rfftfreq(len(num),d=1/sample_freq)

#window data:
hann=np.hanning(len(num))

#window fft:
fft_hann=np.fft.fft(hann) 

#window fft Linear Spectrum:
wlin_spec=sqrt(2)*np.abs(np.fft.rfft(hann))/len(num)

#window + sin
wsin=hann*sin

#window + sin fft:
wsin_spec=sqrt(2)*np.abs(np.fft.rfft(wsin))/len(num)
wsin_spec_freq_axis=np.fft.rfftfreq(len(num),d=1/sample_freq)

fig=plt.figure()
ax1 = fig.add_subplot(431)
ax2 = fig.add_subplot(432)
ax3 = fig.add_subplot(433)
ax4 = fig.add_subplot(434)
ax5 = fig.add_subplot(435)
ax6 = fig.add_subplot(436)
ax7 = fig.add_subplot(413)
ax8 = fig.add_subplot(414)

ax1.plot(num,sin,'r')
ax2.plot(fft_freq_axis,abs(fft_sin),'r')
ax3.plot(lin_spec_freq_axis,lin_spec,'r')
ax4.plot(num,hann,'b')
ax5.plot(fft_freq_axis,fft_hann)
ax6.plot(lin_spec_freq_axis,wlin_spec)
ax7.plot(num,wsin,'c')
ax8.plot(wsin_spec_freq_axis,wsin_spec)

plt.show()

编辑： 如 cmets 中所要求的，我以 dB 标度绘制函数，获得更清晰的图。非常感谢@SleuthEye！

【问题讨论】：

求大神们，请解释一下剧情中每个框是什么，“纯FFT”是什么意思，具体是什么问题。

与您的预期有何不同？第二幅图显示了插入 1 的直流分量在 0 处的峰值。不直接明显的是远离峰值的泄漏水平。如果您以对数刻度（分贝）显示频谱幅度，您会看得更清楚。

@SleuthEye：来自正弦波和窗口+正弦的线性光谱似乎工作正常，问题是汉宁函数的问题，没有显示任何精确，我预计是这样的： link.问候！ :)

错误的问题标题。你得到了一个很好的结果，但期望值很差。 Windows 降低了噪声（在光谱的平坦部分），而不是使峰更精确。

如果您关心ax5 上的图，我看到的主要问题是您正在绘制复数值。尝试绘制abs(fft_hann)。对于ax6，鉴于您拥有的汉宁窗口的带宽（~1Hz），您可能应该更改x轴的比例以查看任何内容（例如ax6.set_xlim([0, 10])）。

【参考方案1】：

看来有问题的情节是由以下人生成的情节：

ax5.plot(fft_freq_axis,fft_hann)

生成图表：

而不是预期的graph from Wikipedia。

情节的构建方式存在许多问题。首先是这个命令本质上是试图绘制一个复值数组（fft_hann）。结果，您实际上可能会收到警告ComplexWarning: Casting complex values to real discards the imaginary part。要生成一个看起来像 Wikipedia 中的图表，您必须使用以下方法获取幅度（而不是实部）：

ax5.plot(fft_freq_axis,abs(fft_hann))

然后我们注意到仍然有一条线贯穿我们的情节。看着np.fft.fft's documentation：

结果中的值遵循所谓的“标准”顺序：如果A = fft(a, n)，则A[0] 包含零频率项（信号的总和），对于实际输入，它始终是纯实数。然后A[1:n/2] 包含正频率项，A[n/2+1:] 包含负频率项，按负频率递减的顺序排列。 [...] 例程np.fft.fftfreq(n) 返回一个数组，给出输出中相应元素的频率。

确实，如果我们打印fft_freq_axis，我们可以看到结果是：

[ 0.  1.  2. ..., -3. -2. -1.]

要解决这个问题，我们只需将数组的上下部分替换为np.fft.fftshift：

ax5.plot(np.fft.fftshift(fft_freq_axis),np.fft.fftshift(abs(fft_hann)))

那么您应该注意，***上的图表实际上以decibels 中的幅度显示。然后你需要做同样的事情：

ax5.plot(np.fft.fftshift(fft_freq_axis),np.fft.fftshift(20*np.log10(abs(fft_hann))))

我们应该越来越近了，但结果并不完全相同，如下图所示：

这是因为 Wikipedia 上的图实际上具有更高的频率分辨率并在其振荡时捕获频谱的值，而您的图在较少的点处对频谱进行采样，并且其中许多点接近于零幅度。为了解决这个问题，我们需要得到窗口在更多频率点的频谱。 这可以通过对 FFT 的输入进行零填充，或者更简单地将参数 n（输出的所需长度）设置为远大于输入大小的值来完成：

N = 8*len(num)
fft_freq_axis=np.fft.fftfreq(N,d=1/sample_freq)
fft_hann=np.fft.fft(hann, N) 
ax5.plot(np.fft.fftshift(fft_freq_axis),np.fft.fftshift(20*np.log10(abs(fft_hann))))
ax5.set_xlim([-40, 40])
ax5.set_ylim([-50, 80])

【讨论】：

你无法想象我有多感激！，我欠你一块巧克力！这就是我无法理解的，非常感谢！