>>> from env_helper import info; info()
页面更新时间： 2023-04-15 21:29:19
运行环境：
Linux发行版本: Debian GNU/Linux 12 (bookworm)
操作系统内核: Linux-6.1.0-7-amd64-x86_64-with-glibc2.36
Python版本: 3.11.2
2.7. 滤波器设计¶
scipy.signal库提供了许多信号处理方面的函数。在这一节,让我们来看看如何利用signal库设计滤波 器,查看滤波器的频率响应,以及如何使用滤波器对信号进行滤波。 假设如下导入signal库:
>>> import scipy.signal as signal
>>> import scipy.optimize as opt
>>> import numpy as np
下面的程序设计一个带通IIR滤波器:
>>> b, a = signal.iirdesign([0.2, 0.5], [0.1, 0.6], 2, 40)
这个滤波器的通带为0.2f0到0.5f0,阻带为小于0.1f0和大于0.6f0,其中f0为1/2的信号取样频率, 如果取样频率为8kHz的话,那么这个带通滤波器的通带为800Hz到2kHz。通带的最大增益衰减为 2dB,阻带的最小增益衰减为40dB,即通带的增益浮动在2dB之内,阻带至少有40dB的衰减。 iirdesgin返回的两个数组b和a, 它们分别是IIR滤波器的分子和分母部分的系数。其中a[0]恒等于1。 下面通过调用freqz计算所得到的滤波器的频率响应:
>>> w, h = signal.freqz(b, a)
freqz返回两个数组w和h,其中w是圆频率数组,通过w/pi*f0可以计算出其对应的实际频率。h是w中 的对应频率点的响应,它是一个复数数组,其幅值为滤波器的增益,相角为滤波器的相位特性。 下面计算h的增益特性,并转换为dB度量。由于h中存在幅值几乎为0的值,因此先用clip函数对其裁剪 之后,再调用对数函数,避免计算出错。
>>> # 本程序用各种fmin函数求卷积的逆运算
>>> 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:
>>> pl.plot(w/np.pi*4000, power)
[<matplotlib.lines.Line2D at 0x7f9e1d9ce5d0>]
在实际运用中为了测量未知系统的频率特性,经常将频率扫描波输入到系统中,观察系统的输出,从 而计算其频率特性。下面让我们来模拟这一过程。 为了调用chirp函数以产生频率扫描波形的数据,首先需要产生一个等差数组代表取样时间,下面的语 句产生2秒钟取样频率为8kHz的取样时间数组:
>>> t = np.arange(0, 2, 1/8000.0)
然后调用chirp得到2秒钟的频率扫描波形的数据:
>>> sweep = signal.chirp(t, f0=0, t1 = 2, f1=4000.0)
频率扫描波的开始频率f0为0Hz,结束频率f1为4kHz,到达4kHz的时间为2秒,使用数组t作为取样时 间点。 下面通过调用lfilter函数计算sweep波形经过带通滤波器之后的结果:
>>> 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] 为了和系统的增益特性图进行比较,需要获取输出波形的包络,因此下面先将输出波形数据转换为能 量值:
>>> out = 20*np.log10(np.abs(out))
为了计算包络,找到所有能量大于前后两个取样点(局部最大点)的下标:
>>> index = np.where(np.logical_and(out[1:-1] > out[:-2], out[1:-1] > out[2:]))[0] + 1
最后将时间转换为对应的频率,绘制所有局部最大点的能量值:
>>> pl.plot(t[index]/2.0*4000, out[index] )
>>> pl.show()
