Scipy库中IIR滤波器的应用

平和少年

已于 2024-04-24 00:09:27 修改

阅读量1.5k

点赞数 40

分类专栏：信号处理文章标签：信号处理

于 2024-04-23 00:40:17 首次发布

本文链接：https://blog.csdn.net/u011572784/article/details/138096712

版权

信号处理专栏收录该内容

43 篇文章 3 订阅

订阅专栏

本文详细讨论了Scipy库中Signal子库如何在IIR滤波器设计中使用，涉及传输函数、ZPK和SOS形式，以及高阶滤波的稳定性问题和分段滤波实践。

摘要由CSDN通过智能技术生成

Scipy函数库是Python中最常用的函数库之一，它里面包含了统计、最优化、插值、积分、线性代数、傅里叶变换、信号/图像处理、稀疏矩阵处理、特殊数学函数等各种模块。在这篇文章中，我们主要总结Scipy函数库中的Signal子库在IIR滤波器设计中的应用。下一篇文章则将对其在FIR滤波器设计中的应用进行说明。

1、IIR 滤波器
（1）IIR滤波基本概念
设输入序列为 $x_n$ ，输出序列为 $y_n$ ，则IIR滤波器可以用差分方程表示为
$a_0 y_n=\sum_{i=0}^{M}b_i x_{n-i}-\sum_{i=1}^{N}a_i y_{n-i} \tag{1}$
显然，为了得到第 $n$ 个输出值，需要用到前 $M$ 个输入值，和前 $N$ 个输出值，即IIR滤波器是一个带反馈的滤波器。
对(1)式两边同时做Z变换，可求得系统函数为
$H(z)=\frac{b_0+b_1 z^{-1}+...+b_M z^{-M}}{a_0+a_1z^{-1}+...+a_Nz^{-N}}=k\frac{\prod_{m=1}^{M}(z-z_m)}{\prod_{n=1}^N(z-z_n)}\tag{2}$
Scipy中生成的滤波器通常 $a_0=1$ ，在上面的表达式中，若 $a_1=a_2=...=a_N=0$ ，则表示FIR滤波器。从(2)中可知，当 $z=z_m$ 时, $H (z) = 0$ ，此时 $z_m$ 称为该滤波器的零点；当 $z=z_n$ 时， $H(z)=\infty$ ，此时 $z_n$ 称为该滤波器的极点； $k$ 称为滤波器的增益。

（2）IIR滤波的表示形式
在Scipy库中，IIR滤波器主要有三种表示形式，下面分别进行说明：

a) 传输函数形式
传输函数是以式(2)中分子和分母的系数来表征的，其中分子系数保存在向量 $b$ 中，分母系数保存在向量 $a$ 中。下面代码生成了4阶butterworth滤波器，其返回值即为式(2)中分子分母的系数向量。由于这种表示形式直接于系统函数相对应，所以最直观也最好理解。但是，这种表示方法的精度通常较低，当滤波器阶数较高时可能产生较大的数值误差。

from scipy.signal import butter
b, a = butter(4, 0.125) 
b
Out[4]: array([0.0009335 , 0.00373399, 0.00560099, 0.00373399, 0.0009335 ])
a
Out[5]: array([ 1.        , -2.97684433,  3.42230953, -1.7861066 ,  0.35557738])

b) 零极点/增益模式(ZPK)
从式(2)中可知，零点和极点分别是使系统函数为0和 $\infty$ 的点。若给定系统零点、极点和增益，则系统函数也唯一确定。Scipy库中，若要使滤波器ZPK模式进行表示，只需令函数的入参output='zpk’即可，下面给出例子：

from scipy.signal import butter
z, p, k = butter(4, 0.125, output='zpk')
z
Out[7]: array([-1., -1., -1., -1.])
p
Out[8]: 
array([0.80586355+0.30839063j, 0.68255862+0.10819419j,
       0.68255862-0.10819419j, 0.80586355-0.30839063j])
k
Out[9]: 0.0009334986129548442

ZPK表示方式的一个缺陷是，Scipy库中并没有可以直接用zpk进行滤波器的库函数，所以首先需要将zpk形式转化为传输函数形式(b,a)或后面会说到的SOS模式。

c) SOS模式
SOS模式的全称是"Second Order Sections"，它将原始滤波器转化成多个二阶滤波器的级联，所以它的返回值是一个 $(n, 6)$ 的矩阵。其中，每一行均保存了一个二阶滤波器的传输函数，每一行的前3个元素表示二阶滤波器传输函数中的分子系数，每一行的后3个元素表示二阶滤波器传输函数中的分母系数。示例如下：

from scipy.signal import butter
sos = butter(4, 0.125, output='sos')
sos
Out[12]: 
array([[ 9.33498613e-04,  1.86699723e-03,  9.33498613e-04,
         1.00000000e+00, -1.36511724e+00,  4.77592250e-01],
       [ 1.00000000e+00,  2.00000000e+00,  1.00000000e+00,
         1.00000000e+00, -1.61172710e+00,  7.44520838e-01]])

根据上面的代码可以看出，一个4阶的滤波器会被拆成两个2阶滤波器的级联形式，所以输出sos为 $(2, 6)$ 的矩阵形式。采用SOS的一个好处是它将高阶滤波器拆成低阶滤波器的级联来实现，所以它的数值稳定性更强，所以一般高阶滤波器的设计时常用SOS表示形式。它的一个缺点是，目前直接用sos进行滤波的函数sosfilt(sos,x)的耗时往往比lfilter(b,a,x)的耗时长。

d) 各种表示形式的相关转换
Scipy库中提供了上述三种表示形式的相互转换的函数：sos2tf, sos2zpk, ss2tf, ss2zpk, tf2sos, tf2zz, tf2zpk, zpk2sos, zpk2ss, zpk2tf。各个函数的具体用法可参考手册。

e) 系数精度有限造成数值计算不稳定的例子说明
在上面介绍传输函数表示形式时，说到它的一个缺点是在高阶滤波器设计时数值稳定性较差。下面我们通过一个例子直观展示该现象。我们生成一个10阶的butterworth带通滤波器，并观察它的阶跃响应，代码如下：

from scipy.signal import butter
from scipy.signal import lfilter
import matplotlib.pyplot as plt
b, a = butter(10, [0.04, 0.16], btype='bandpass')
x = np.ones(125)
y = lfilter(b, a, x)
plt.figure()
plt.plot(y)
plt.show()

在这里插入图片描述

图1.传输函数表示形式的高阶滤波器输出结果

仿真结果如上图所示，显然结果出现极大震荡，是不稳定的。为了探究产生上述现象的原因，我们将传输函数的表达形式转换为ZPK模式(极点都在单位圆内的系统是稳定系统，否则系统是不稳定的)，结果如下，可见有4个极点是在单位圆外的，所以此时系统是不稳定的。

z, p, k = tf2zpk(b, a)
np.abs(p)
Out[21]: 
array([0.95521728, 0.95521728, 1.10894705, 1.10894705, 1.07918597,
       1.07918597, 0.99951079, 0.99951079, 0.8763112 , 0.8763112 ,
       0.85254882, 0.85254882, 0.80590245, 0.80590245, 0.75864141,
       0.75864141, 0.74397519, 0.74397519, 0.76461253, 0.76461253])

我们用SOS表示形式来实现同样的滤波器，代码如下：

from scipy.signal import butter
from scipy.signal import lfilter
from scipy.signal import sosfilt
import matplotlib.pyplot as plt
sos = butter(10, [0.04, 0.16], btype='bandpass', output='sos')
x = np.ones(125)
y = sosfilt(sos, x)
plt.figure()
plt.plot(y)
plt.show()

在这里插入图片描述

图2. SOS表示形式的高阶滤波器输出结果

从上面仿真结果可见，此时波动小很多，系统稳定性得到极大提升，同样地，我们可以观察此时的极点分布情况，可见所有极点都在单位圆内，所以此时系统是稳定的。

z, p, k = sos2zpk(sos)
np.abs(p)
Out[25]: 
array([0.78818242, 0.78818242, 0.79966393, 0.79966393, 0.81809243,
       0.81809243, 0.85431478, 0.85431478, 0.87659171, 0.87659171,
       0.9025228 , 0.9025228 , 0.93645496, 0.93645496, 0.95535599,
       0.95535599, 0.96368222, 0.96368222, 0.98809157, 0.98809157])

（3）IIR滤波器设计示例

a) 低通滤波器
下图中蓝色波形是一个由多个不同频率的正弦信号叠加而成的混合信号。采样率 $f_s=10000Hz$ ，有效信号频率为 $f_0=100Hz$ ，其他干扰信号的频率均大于100Hz，故可以用低通滤波器来滤除干扰信号，提取出有效信号。
在这里插入图片描述

图3. 滤波前后波形对比

如下面的代码所示，首先实现截止频率以Hz为单位的滤波器，然后调用该函数进行滤波器，滤波结果如上图红线所示，有效信号被正确地提取出来了。需要注意的是，我们调用的sosfiltfilt这个滤波函数是一个前/后向滤波函数，即它会分别从前向/后向这两个相反的方向进行滤波，这样可以取得零相位移动的效果。

from scipy.signal import butter, sosfiltfilt
import matplotlib.pyplot as plt

def butter_lowpass(cutoff, fs, order):
    normal_cutoff = cutoff / (0.5 * fs)
    sos = butter(order, normal_cutoff, btype='low', output='sos')
    return sos
    
def butter_lowpass_filtfilt(data, cutoff, fs, order):
    sos = butter_lowpass(cutoff, fs, order)
    y = sosfiltfilt(sos, data)
    return y

if __name__ == "__main__":
    fs = 10000
    f0 = 100
    N = int(10000)
    t = np.arange(0, N, 1) / fs
    sig = np.sin(2 * np.pi * f0 * t)
    fi = np.arange(500, 3000, 500)
    for i in range(len(fi)):
        sig += 0.2 * np.sin(2 * np.pi * fi[i] * t)
    y = butter_lowpass_filtfilt(sig, cutoff=300, fs=fs, order=6)

    plt.figure()
    plt.plot(t, sig, 'b')
    plt.plot(t, y, 'r')
    plt.legend(['original signal', 'filtered signal'])
    plt.show()

滤波过程需要注意，在函数lfilter、sosfilt等函数的滤波过程中，一般均假定初始状态为0，但是如果输入信号的初始状态不为0，则滤波结果中会出现瞬态变化的形态。下面通过一个仿真说明这个问题，仿真代码如下：

from scipy.signal import butter, sosfiltfilt, sosfilt, sosfilt_zi
import matplotlib.pyplot as plt
if __name__ == "__main__":
    n = 101
    t = np.linspace(0, 1, n)
    x = 0.45 + 0.1 * np.random.randn(n)
    sos = butter(8, 0.125, output='sos')

    y = sosfilt(sos, x)
    zi = x[:4].mean() * sosfilt_zi(sos)
    y2, _ = sosfilt(sos, x, zi=zi)

    plt.figure()
    plt.plot(t, x, 'b')
    plt.plot(t, y, 'k')
    plt.plot(t, y2, 'r')
    plt.legend(['original signal', 'filtered signal with intinal 0', 'filtered signal with intinal others'])
    plt.show()

从下面的仿真结果可以看出，若不设置初值，则滤波默认初值为0，此时初始滤波阶段会有一段瞬态响应，将滤波初始设置得更符合实际情况后，该瞬态效应消失。需要注意的是，函数filtfilt和函数sosfiltfilt中已经考虑了滤波初值的问题，所以不用再进行额外的考虑。
在这里插入图片描述

图4. 不同滤波初值的滤波结果对比

高通/带通滤波器的设计可以参考上述低通滤波的设计，这边不再赘述。

b) 对长信号进行分段滤波
在某些场景下，我们希望对一长段信号进行滤波，但是却没有办法一次得到完整的数据，只能一段一段地进行处理。为了和一次长段处理的结果保持一致，需要把每个小段处理后的结果保存下来，作为下一小段处理的初始条件，下面给出一个例子进行说明，仿真代码如下：

from scipy.signal import butter, sosfiltfilt, sosfilt, sosfilt_zi
import matplotlib.pyplot as plt
if __name__ == "__main__":
    fs = 10000
    f0 = 100
    N = int(10000)
    batch_num = 100
    batch_len = int(N / batch_num)
    t = np.arange(0, N, 1) / fs
    sig = np.sin(2 * np.pi * f0 * t)
    fi = np.arange(500, 3000, 500)
    for i in range(len(fi)):
        sig += 0.2 * np.sin(2 * np.pi * fi[i] * t)

    sos = butter(6, 300 / (0.5 * fs), btype='low', output='sos')
    y = sosfilt(sos, sig)

    y2 = np.zeros(len(sig))
    z = np.zeros((sos.shape[0], 2))
    start = 0
    while start < len(sig):
        stop = np.min(np.array([start + batch_len, len(sig)]))
        y2[start:stop], z = sosfilt(sos, sig[start:stop], zi=z)
        start = stop

    plt.figure()
    plt.plot(t, sig, 'b')
    plt.plot(t, y, 'ko-')
    plt.plot(t, y2, 'r.-')
    plt.legend(['original signal', 'long filtered signal', 'batch filtered signal'])
    plt.show()

从下图可见，按照上面的处理方式，分段滤波结果和整段滤波结果是一致的。
在这里插入图片描述

图5. 分段滤波结果和整段滤波结果对比

暂时就总结这些，后面有些新的重要的应用再不断往这里补充完善。

平和少年

关注

40
点赞
踩
6

收藏

觉得还不错? 一键收藏
0
评论
Scipy库中IIR滤波器的应用

为了探究产生上述现象的原因，我们将传输函数的表达形式转换为ZPK模式(极点都在单位圆内的系统是稳定系统，否则系统是不稳定的)，结果如下，可见有4个极点是在单位圆外的，所以此时系统是不稳定的。其中，每一行均保存了一个二阶滤波器的传输函数，每一行的前3个元素表示二阶滤波器传输函数中的分子系数，每一行的后3个元素表示二阶滤波器传输函数中的分母系数。从上面仿真结果可见，此时波动小很多，系统稳定性得到极大提升，同样地，我们可以观察此时的极点分布情况，可见所有极点都在单位圆内，所以此时系统是稳定的。
复制链接

扫一扫