python中的fft带通滤波器
声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow
原文地址: http://stackoverflow.com/questions/19122157/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
fft bandpass filter in python
提问by men in black
What I try is to filter my data with fft. I have a noisy signal recorded with 500Hz as a 1d- array. My high-frequency should cut off with 20Hz and my low-frequency with 10Hz. What I have tried is:
我尝试的是用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)
What I get now are complex numbers. So something must be wrong. What? How can I correct my code?
我现在得到的是复数。所以一定是出了什么问题。什么?如何更正我的代码?
采纳答案by Hooked
It's worth noting that the magnitude of the units of your bp
are not necessarily going to be in Hz, but are dependent on the sampling frequency of signal, you should use scipy.fftpack.fftfreq
for the conversion. Also if your signal is real you should be using scipy.fftpack.rfft
. Here is a minimal working example that filters out all frequencies less than a specified amount:
值得注意的是,您的单位的大小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)
We can plot the evolution of the signal in real and fourier space:
我们可以绘制信号在实空间和傅立叶空间中的演变:
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()
回答by Paul R
There's a fundamental flaw in what you are trying to do here - you're applying a rectangular window in the frequency domain which will result in a time domain signal which has been convolved with a sinc function. In other words there will be a large amount of "ringing" in the time domain signal due to the step changes you have introduced in the frequency domain. The proper way to do this kind of frequency domain filtering is to apply a suitable window functionin the frequency domain. Any good introductory DSP book should cover this.
您在这里尝试做的事情存在一个基本缺陷 - 您在频域中应用一个矩形窗口,这将产生一个与 sinc 函数卷积的时域信号。换句话说,由于您在频域中引入的阶跃变化,时域信号中将出现大量“振铃”。进行这种频域滤波的正确方法是在频域中应用合适的窗函数。任何优秀的介绍性 DSP 书籍都应该涵盖这一点。