3.6.2 滤波器设计

signal模块提供了许多滤波器设计的函数。在下面的实例中，我们设计一个IIR带通滤波器，并查看其频率响应，最后使用它对频率扫描信号进行滤波计算。

    sampling_rate = 8000.0
    
    # 设计一个带通滤波器：
    # 通带为0.2*
4000 - 0.5*
4000
    # 阻带为<0.1*
4000, >0.6*
4000
    # 通带增益的最大衰减值为2dB
    # 阻带的最小衰减值为40dB
    b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40) ❶
    
    # 使用freq计算滤波器的频率响应
    w, h = signal.freqz(b, a) ❷
    
    # 计算增益
    power = 20*
np.log10(np.clip(np.abs(h), 1e-8, 1e100)) ❸
    freq = w / np.pi *
 sampling_rate / 2

❶首先用iirdesign( )设计一个IIR带通滤波器。这个滤波器的通带为0.2f0到0.5f0，阻带为小于0.1f0和大于0.6f0，其中f0为信号取样频率的一半。如果取样频率为8kHz，那么这个带通滤波器的通带为800Hz到2kHz。通带的最大增益衰减为2dB，阻带的最小增益衰减为40dB，即通带的增益浮动在2dB之内，阻带至少有40dB的衰减。

iirdesgin( )返回两个数组b和a，它们分别是IIR滤波器的分子和分母部分的系数。其中a[0]恒等于1。❷调用freqz( )计算所得到的滤波器的频率响应。freqz( )返回两个数组w和h，其中w是圆频率数组，通过ωf0/π可以计算出与其对应的实际频率。h是w中对应频率点的响应，它是一个复数数组，其幅值表示滤波器的增益特性，相角表示滤波器的相位特性。

❸计算h的增益特性，并使用dB进行度量。由于h中存在幅值几乎为0的值，因此先用clip( )对其裁剪之后，再调用对数函数，避免计算出错。

在实际运用中为了测量未知系统的频率特性，经常将频率扫描波输入到系统中，观察系统的输出，从而计算其频率特性。下面让我们模拟这一过程：

    # 产生两秒钟的取样频率为sampling_rate Hz的频率扫描信号
    # 开始频率为0，结束频率为sampling_rate/2
    t = np.arange(0, 2, 1/sampling_rate) ❶
    sweep = signal.chirp(t, f0=0, t1=2, f1=sampling_rate/2) ❷
    # 对频率扫描信号进行滤波
    out = signal.lfilter(b, a, sweep) ❸
    # 将波形转换为能量
    out = 20*
np.log10(np.abs(out)) ❹
    # 找到所有局部最大值的下标
    index = signal.argrelmax(out, order=3)  ❺
    # 绘制滤波之后的波形的增益
    pl.figure(figsize=(8, 2.5))
    pl.plot(freq, power, label=u"带通IIR滤波器的频率响应") 
    pl.plot(t[index]/2.0*
4000, out[index], label=u"频率扫描波测量的频谱", alpha=0.6) ❻
    pl.legend(loc="best")

❶为了调用chirp( )产生频率扫描波形的数据，首先需要产生一个表示取样时间的等差数组，这里产生两秒的取样频率为8kHz的取样时间数组。❷然后调用chirp( )得到两秒的频率扫描波形的数据。频率扫描波的开始频率f0为0Hz，结束频率f1为4kHz，到达4kHz的时间为两秒，使用数组t作为取样时间点。❸最后调用lfilter( )计算频率扫描波形经过带通滤波器之后的结果。

❹为了和系统的增益特性图进行比较，需要获取输出波形的包络，因此先将输出波形数据转换为能量值。❺为了计算包络，调用argrelmax( )找到out数组中所有局部最大值的下标，order参数指定局域最大值的范围，这里的值为3表示所有的局域最大值都是连续7个元素（前后各三个元素）中的最大值。❻最后将时间转换为对应的频率，绘制所有局部最大点的能量值。

图3-32显示了freqz( )计算的频谱和频率扫描波得到的频率特性，可以看到结果是一致的。

图3-32 用频率扫描波测量的频率响应