In [2]:
import scipy.signal as signal
import scipy.optimize as opt
import numpy as np
下面的程序设计一个带通IIR滤波器：
In [3]:
b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40)
In [4]:
w, h = signal.freqz(b, a)
freqz 返回两个数组w和h，其中w是圆频率数组，通过 w/pi*f0 可以计算出其对应的实际频率。
h是w中的对应频率点的响应，它是一个复数数组，其幅值为滤波器的增益，相角为滤波器的相位特性。
下面计算h的增益特性，并转换为dB度量。由于h中存在幅值几乎为0的值，
因此先用 clip 函数对其裁剪之后，再调用对数函数，避免计算出错。
本程序用各种 fmin 函数求卷积的逆运算
In [5]:
import scipy.signal as signal
import numpy as np
import scipy as sp
import pylab as pl
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))
通过下面的语句可以绘制出滤波器的增益特性图，这里假设取样频率为8kHz：
In [6]:
pl.plot(w/np.pi*4000, power)
Out[6]:
[<matplotlib.lines.Line2D at 0x7f3486a856a0>]
在实际运用中为了测量未知系统的频率特性，经常将频率扫描波输入到系统中，观察系统的输出， 从而计算其频率特性。下面让我们来模拟这一过程。
为了调用 chirp 函数以产生频率扫描波形的数据，首先需要产生一个等差数组代表取样时间，
下面的语句产生2秒钟取样频率为8kHz的取样时间数组:
In [7]:
t = np.arange(0, 2, 1/8000.0)
然后调用 chirp 得到2秒钟的频率扫描波形的数据:
In [8]:
sweep = signal.chirp(t, f0=0, t1 = 2, f1=4000.0)
In [9]:
out = signal.lfilter(b, a, sweep)
lfilter 内部通过如下算式计算IIR滤波器的输出:
通过如下算式可以计算输入为x时的滤波器的输出，其中数组x代表输入信号，y代表输出信号:
y[n] = b[0]x[n] + b[1]x[n − 1] + · · · + b[P ]x[n − P ] − a[1]y[n − 1] − a[2]y[n − 2] − · · · − a[Q]y[n − Q]
为了和系统的增益特性图进行比较，需要获取输出波形的包络，因此下面先将输出波形数据转换为能量值:
In [10]:
out = 20*np.log10(np.abs(out))
为了计算包络,找到所有能量大于前后两个取样点(局部最大点)的下标:
In [11]:
index = np.where(np.logical_and(out[1:-1] > out[:-2], out[1:-1] > out[2:]))[0] + 1
最后将时间转换为对应的频率,绘制所有局部最大点的能量值:
In [12]:
pl.plot(t[index]/2.0*4000, out[index] )
pl.show()
