用汉宁窗傅里叶滤波后恢复信号

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 的结果将是 信号乘以汉宁窗。因此，您提供的第三个情节是有道理的。它是原始信号，乘以汉宁窗（并通过消除高于 2.0 的频率进行平滑）。

我在网上看到我应该在进行傅里叶变换之前应用汉宁窗以避免泄漏。

汉宁窗用于尽可能减少频谱泄漏。在理想情况下，您可以将有限的数据片段用作图块，在所有空间中重新创建整个函数。当您有限的数据片段不能完全满足该目的时，您将产生一定量的频谱泄漏。

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

有关频谱泄漏的更多信息，请转至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