python中的fft带通滤波器

Posted 2023-03-10

【中文标题】python中的fft带通滤波器

【英文标题】：fft bandpass filter in python

【发布时间】：2013-10-07 23:10:09

【问题描述】：

我尝试使用 fft 过滤我的数据。我有一个以 500Hz 作为一维数组记录的噪声信号。我的高频应该用 20Hz 截止，我的低频应该用 10Hz 截止。 我试过的是：

fft=scipy.fft(signal) 
bp=fft[:]  
for i in range(len(bp)): 
    if not 10<i<20:
        bp[i]=0

ibp=scipy.ifft(bp)

我现在得到的是复数。所以一定有什么地方不对劲。什么？如何更正我的代码？

【参考方案1】：

您在这里尝试执行的操作存在一个根本缺陷 - 您在频域中应用了一个矩形窗口，这将产生一个与 sinc 函数卷积的时域信号。换句话说，由于您在频域中引入的阶跃变化，时域信号中会出现大量“振铃”。进行这种频域过滤的正确方法是在频域中应用合适的window function。任何优秀的 DSP 入门书籍都应该涵盖这一点。

【讨论】：

+1。相关问题：dsp.stackexchange.com/questions/6220/…

【参考方案2】：

值得注意的是，bp 的单位幅度不一定是赫兹，而是取决于信号的采样频率，你应该使用scipy.fftpack.fftfreq 进行转换。此外，如果您的信号是真实的，您应该使用scipy.fftpack.rfft。这是一个最小的工作示例，可以过滤掉所有小于指定数量的频率：

import numpy as np
from scipy.fftpack import rfft, irfft, fftfreq

time   = np.linspace(0,10,2000)
signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time)

W = fftfreq(signal.size, d=time[1]-time[0])
f_signal = rfft(signal)

# If our original signal time was in seconds, this is now in Hz    
cut_f_signal = f_signal.copy()
cut_f_signal[(W<6)] = 0

cut_signal = irfft(cut_f_signal)

我们可以绘制信号在实数和傅立叶空间中的演变：

import pylab as plt
plt.subplot(221)
plt.plot(time,signal)
plt.subplot(222)
plt.plot(W,f_signal)
plt.xlim(0,10)
plt.subplot(223)
plt.plot(W,cut_f_signal)
plt.xlim(0,10)
plt.subplot(224)
plt.plot(time,cut_signal)
plt.show()

【讨论】：

谢谢！但是如果我想要一个带通滤波器，我必须做两次吗？一次低通，一次高通？

@meninblack 当然。正如您在示例中看到的那样，我只用矩形截止过滤了低频（有关此问题的注释和改进，请参见 Paul R 的回答）。我试图做的是给你一个简单的，最小的工作答案，让你看到哪里出了问题（错误的单位，等等......）。