python实现 快速傅里叶变换(fft)、滤波器设计(butter/chebv)

频谱分析

fft的意义在于频谱分析，基于python实现的fft、filter代码如下。其中滤波器分别用了巴特沃斯和契比雪夫二个设计，参数见源码的解释。比较复杂的是契比雪夫的第二个参数的意义。官方原文的解释为The maximum ripple allowed below unity gain in the passband. Specified in decibels, as a positive number.

---

from scipy import signal
import numpy as np
import math
import matplotlib.pyplot as pl

def fft(yt, sampling_rate, fft_size=None):  
    if fft_size is None:  
        fft_size = len(yt)  
    yt = yt[:fft_size]  
    yf = abs(np.fft.rfft(yt)/fft_size)  
    freqs = np.linspace(0, 1.0*sampling_rate/2, 1.0*fft_size/2+1)  
    return freqs, yf  

fs = 500 # sample rate
nyq = fs * 0.5 # nyqust rate, the maximun rate that can be rebuild
x = np.linspace(0,1,fs) # 1 second timeline
yt = np.array([math.cos(100.0*2*math.pi*i) + 2.0*math.sin(120.0*2*math.pi*i) for i in x])

pl.subplot(221)
pl.plot(x,yt)
pl.title("original")

freqs, yf = fft(yt, fs)
pl.subplot(222)
pl.plot(freqs, yf)
pl.title("fft")

b, a = signal.butter(4, 110/nyq, "lowpass") # lowpass
# 阶数；最大纹波允许低于通频带中的单位增益。以分贝表示，以正数表示；频率(Hz)/奈奎斯特频率（采样率*0.5）
b, a = signal.cheby1(4, 5, 110/nyq, "lowpass") # lowpass
yt  = signal.filtfilt(b, a, yt)
pl.subplot(223)
pl.plot(x,yt)
pl.title("lp")

freqs, yf = fft(yt, fs)
pl.subplot(224)
pl.plot(freqs, yf)
pl.title("fft")

pl.show()