滤波器实现

卷积和滤波

滤波的数学基础是卷积。对于有限脉冲响应(FIR)滤波器，滤波运算的输出$y(k)$是输入信号$x(k)$与脉冲响应$h(k)$的卷积： $$ y(k)=\sum_{l=-\infty}^{\infty}h(l)x(k-l) $$

如果输入信号也是有限长度的，您可以使用conv/convolve函数来执行滤波运算。例如，要用三阶平均值滤波器对包含五个样本的随机向量进行滤波，可以将$x(k)$存储在向量x中，将$h(k)$存储在向量h中，并求这两个向量的卷积：

import numpy as np
from scipy import signal

x = np.random.randn(5)
h = np.ones(4)/4
y = signal.convolve(x,h)

y的长度比x和h的长度之和小1。

滤波器和传递函数

滤波器的传递函数是其脉冲响应的Z变换。对于FIR滤波器，输出y的Z变换$Y(z)$是传递函数和输入x的Z变换$X(z)$的乘积：

$$Y(z)=H(z)X(z)=\left(h(1)+h(2)z^{-1}+\cdots+h(n+1)z^{-n}\right)X(z)$$

多项式系数$h(1),h(2),...,h(n+1)$对应于第n阶滤波器的脉冲响应的系数。

注意：Matlab中滤波器系数索引从1到(n + 1)，而Python中是从0到n。

FIR 滤波器也称为全零、非递归或移动平均值(MA)滤波器。

对于无限脉冲响应(IIR)滤波器，传递函数不是多项式，而是有理函数。输入和输出信号的Z变换的关系是：

$$Y(z)=H(z)X(z)=\frac{b(1)+b(2)z^{-1}+\cdots+b(n+1)z^{-n}}{a(1)+a(2)z^{-1}+\cdots+a(m+1)z^{-m}}X(z)$$

其中$b(i)$和$a(i)$是滤波器系数。在本例中，滤波器的阶是n和m的最大值。n = 0的IIR滤波器也称为全极点、递归或自回归(AR)滤波器。n和m均大于零的IIR滤波器也称为极点-零点、递归或自回归移动平均值(ARMA)滤波器。缩写AR、MA和ARMA通常应用于与滤波随机过程相关联的滤波器。

使用filter函数进行滤波

对于IIR滤波器，滤波运算不能用简单的卷积来说明，而需要用可从传递函数关系中找到的差分方程来说明。假设$a(1)=1$，将分母移到左侧，并进行逆Z变换，以获得

$$y(k)+a(2)y(k-1)+\cdots+a(m+1)y(k-m)=b(1)x(k)+b(2)x(k-1)+\cdots+b(n+1)x(k-n)$$

对于当前输入、过去的输入以及过去的输出，$y(k)$是

$$y(k)=b(1)x(k)+b(2)x(k-1)+\cdots+b(n+1)x(k-n)-a(2)y(k-1)-\cdots-a(m+1)y(k-m)$$

这是数字滤波器的标准时域表示。从$y(1)$开始，假设一个初始条件为零的因果系统，表示形式等效于

$y(1)=b(1)x(1)$
$y(2)=b(1)x(2)+b(2)x(1)-a(2)y(1)$
$y(3)=b(1)x(3)+b(2)x(2)+b(3)x(1)-a(2)y(2)-a(3)y(1)$
$\quad ; ; ;\vdots$
$y(n)=b(1)x(n)+\cdots+b(n)x(1)-a(2)y(n-1)-\cdots-a(n)y(1)$

要实现此滤波运算，您可以使用MATLAB的filter函数或Python的lfilter函数。filter/lfilter将系数存储在两个行向量中，其中一个为分子，一个为分母。例如，要求解差分方程

$$y(n)-0.9y(n-1)=x(n)\quad \Rightarrow \quad Y(z)=\frac{1}{1-0.9z^{-1}}X(z)=H(z)X(z)$$

b = np.array([1])
a = np.array([1,-0.9])
zi = signal.lfilter_zi(b,a)*0
y,_ = signal.lfilter(b,a,x,zi=zi)

filter/lfilter提供的输出样本的数量与输入样本一样多，即y的长度与x的长度相同。如果a的第一个元素不是1，则filter/lfilter在实现差分方程之前，将系数除以$a(1)$。