用scipy.signal设计IIR(无线脉冲响应滤波器)

import scipy.signal as signal
import pylab as pl
import numpy as np

b,a=signal.iirdesign([0.2,0.5],[0.1,0.6],2,40)
w,h=signal.freqz(b,a)
power=20*np.log10(np.clip(np.abs(h),1e-8,1e100)) #对幅值几乎为0的h裁剪，并使用dB来度量
pl.plot(w/np.pi*4000,power)

t=np.arange(0,2,1/8000.0) #2s钟、取样时间为8KHz的取样时间数组
sweep=signal.chirp(t,f0=0,t1=2,f1=4000.0) #扫描波形
pl.figure(2)
pl.plot(t,20*np.log10(sweep))

out=signal.lfilter(b,a,sweep) #扫描频率从带通滤波器出来后的结果
out=20*np.log10(np.abs(out))  #转化为dB
#获取输出波形的包络
#找到所有能量大于前后两个取样点(局部最大点)的小标,out必须为数组
index=np.where(np.logical_and(out[1:-1]>out[:-2],out[1:-1]>out[2:]))[0]+1
pl.figure(3)
pl.plot(t[index]/2.0*4000,out[index])#将时间转换为对应的频率，绘制所有局部最大点的能量值

'''
滤波器带通为0.2f0到0.5f0，阻带为小于0.1f0以及大于0.6f0，带通的最大增益衰减为2dB，阻带的
最小增益衰减为40dB
f0为信号采样频率的一半
b为iir滤波器的分子系数
a为iir滤波器的分母系数,a[0]恒为1

freqz(b,a) 计算所得到的滤波器的频率响应
w是原频率数组，其对应的实际频率为w*f0/pi,f0=fs/2
h是w对应频率点的响应，h是一个复数，其幅值表示增益特性，相角表示滤波器的相位特性

clip函数，第二个参数为最小值，第三个参数为最大值
数组中的元素小于最小值则被改为最小值，大于最大值则被改为最大值。

chirp函数，第一个参数为取样时间，第二个参数为频率扫描波的开始频率f0，第四个参数为
结束频率f1，第三个参数为f0到f1的时间。
注：需要学习chirp信号

lfilter函数计算波形通过滤波器后的结果
'''