python中的fft带通滤波器

Question

我尝试用fft过滤我的数据.我有一个嘈杂的信号记录500Hz作为1d阵列.我的高频应切断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)

我现在得到的是复数.所以一定是错的.什么？我该如何更正我的代码？

Answer 1

值得注意的是,你的单位的大小bp不一定是以Hz为单位,而是取决于信号的采样频率,你应该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()

Answer 2

你在这里尝试做的事情有一个根本的缺陷 - 你在频域中应用一个矩形窗口,这将产生一个时域信号,该信号已经与sinc函数进行了卷积.换句话说,由于您在频域中引入的步骤变化,时域信号中将存在大量"振铃".进行这种频域滤波的正确方法是在频域中应用合适的窗函数.任何好的入门DSP书都应该涵盖这一点.