matlab - 带汉宁窗的傅立叶滤波器后恢复信号

Question

我有一个复杂的信号，可以在傅立叶空间中看到，并且想要过滤一些我不想要的频率。我在网上阅读了我在进行傅立叶变换之前应该应用Hanning窗口的方法，以避免泄漏。

因此，如下面的代码所示，我正在做的是将汉宁窗口应用于我的数据，然后对其进行傅立叶变换。作为测试，我想看看是否不过滤任何东西，是否可以恢复原始信号。但是，信号在边缘变为零。

现在，我了解到这是因为汉宁窗口过滤器的末尾也变为零。在那种情况下，如何应用汉宁窗口，进入频域并恢复信号时回到我的时域？如果我的信号在末尾变为零，则当我尝试过滤所需的频率时，时域中的结果在边缘仍将变为零。

我的方法缺少什么/做错了什么？感谢您提供的任何帮助!

这是我在做什么的示例代码:

import sys
import matplotlib
def fourier(time,array):

    fft = np.fft.fft(array*np.hanning(len(array)))

    Npts = len(array)
    spacing_array = time[::-1][:-1][::-1] - time[:-1]

    if np.mean(spacing_array) - spacing_array[0] > 1.e-16:
            print "time axis not equally separated. cannot compute fft"
            sys.exit()
    spacing = spacing_array[0]

    freq = np.fft.fftfreq(Npts, spacing)

    return freq,fft

if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt

    # generate a sample signal
    sample_rate = 100.0
    nsamples = int(2.e3)
    t = np.arange(nsamples) / sample_rate
    x = np.cos(2*np.pi*0.5*t) + 0.2*np.sin(2*np.pi*2.5*t+0.1) + \
            0.2*np.sin(2*np.pi*15.3*t) + 0.1*np.sin(2*np.pi*16.7*t + 0.1) + \
                0.1*np.sin(2*np.pi*23.45*t+.8)
    y = np.cos(7*np.pi*1.1*t) + 3*np.sin(0.2*np.pi*2.8*t+1) + \
            0.2*np.sin(8*np.pi*2.7*t) + 1*np.sin(2*np.pi*t + 2.1) + \
                0.1*np.sin(0.2*np.pi*0.45*t+1.4)
    z = x*np.exp(y*1j)

    z_freq,z_fft = fourier(t,z)

    plt.clf()
    plt.figure(figsize=(8,12))

    plt.subplot(4,1,1) # original signal
    plt.plot(t,np.absolute(z))

    plt.subplot(4,1,2) # fourier transform
    plt.semilogy(sorted(z_freq),[b for (a,b) in sorted(zip(z_freq,np.absolute(z_fft)/nsamples))])

    # filtering
    plt.subplot(4,1,3)
    idx = np.where(np.abs(z_freq)>2.0)
    z_fft[idx]=0
    z_filter = np.fft.ifft(z_fft)
    plt.plot(t,np.real(z_filter))

    z_freq,z_fft = fourier(t,z_filter)
    plt.subplot(4,1,4)
    plt.semilogy(sorted(z_freq),[b for (a,b) in sorted(zip(z_freq,np.absolute(z_fft)/nsamples))])
    plt.show()

从中得出的图像如下:
Real space and Frequency space of signal before and after filtering with Hanning window

Answer 1

刺杀您的几个问题。

在那种情况下，如何应用汉宁窗口，进入频域并恢复信号时回到我的时域？

达到numpy的fft实现的数值精度后，先进行傅立叶变换，然后再进行傅立叶逆变换，即可返回原始信号。但是，如果将信号乘以fft之前的汉宁窗口，则fft'ing和ifft'ing的结果将是信号乘以hanning窗口。因此，您提供的第三个情节是有意义的。它是原始信号，乘以汉宁窗口(并通过消除高于2.0的频率进行平滑)。

我在网上阅读了我在进行傅立叶变换之前应该应用Hanning窗口的方法，以避免泄漏。

hanning窗口用于尽可能减少频谱泄漏。在理想情况下，您可以将有限的数据片段用作切片，以在整个空间中重新创建整个函数。当有限的数据片段无法完全满足此目的时，您将招致一定数量的光谱泄漏。

汉宁窗口会衰减信号的开始和结束，以帮助满足“重复图块”约束。话虽如此，我认为频谱泄漏不是您所提供数据的问题，但也许是您正在分析的另一个数据集中的问题。

有关频谱泄漏的更多信息，请转到why does spectral leakage arise in a fft。

如果我的信号在末尾变为零，则当我尝试过滤所需的频率时，时域中的结果在边缘仍将变为零。

不一定是这种情况。您可能在频域中消除的频率实际上是使功能在时域为零的模式。

第四情节

您的第四个情节很有趣。在其中，您将绘制高频消除时间序列的傅立叶变换结果。该功能会以+/- 2的速度快速衰减，因为您已经消除了这些频率。 fft不能精确到0的事实是因为离散傅立叶变换不精确。

快速编码注释

您的代码的一行有趣地使用了数组切片。我花了一秒钟，但我想评论一下你所做的事情。

time[::-1][:-1][::-1] # Reverse array, disregard "new" last element, ...
                      # ... then reverse again
time[1:]  # Don't consider the first element.
time[1:] == time[::-1][:-1][::-1] # These are exactly the same.
# Result: True