滤波器设计、自适应滤波、匹配滤波方法及应用

滤波器设计、自适应滤波、匹配滤波方法及应用

滤波器是一个系统，它对输入系统的信号起到滤波的作用。滤波器设计是建立一个信号系统，使得信号通过该系统处理后可以得到所期望的输出。滤波的目的不同导致信号滤波处理的方法各种各样。自适应滤波器是指根据环境的变化，使用自适应算法来改变滤波器的参数和结构的滤波器，LMS(Least Mean Square，最小均方)和RLS(Recursive Least Squares，递归最小二乘) 自适应滤波算法就是根据估计误差实时调节线性滤波器参数的算法。对于给定幅频特性的线性滤波器，描述了FIR(有限冲击响应)和IIR(无限冲击响应)数字滤波器、梳状滤波器、平滑微分器、Savitzky-Golay平滑器。简述了三种特殊滤波器，即：周期信号增强或抵消、自定义模板匹配滤波、余弦变换滤波。最后用实例说明了滤波器的使用及效果。

(1)自适应滤波原理

设滤波器的输入随机序列为信号与噪声的叠加：x(t) = s(t) + w(t)，s(t)表示信号的真值，w(t)表示噪声，滤波器输出y(t)是信号s(t)估值，自适应滤波器工作过程可用下图表示。

(2)LMS自适应滤波算法

自适应滤波器的权系数随着输入随机序列而不断地得到调整，LMS(Least Mean Square，最小均方)算法就是其中的一种简单的算法，它可以描述如下：

一般我们选择0<μ<1/EX,，EX是输入信号的总功率，当μ越大，收敛速度越快；当μ越小，则收敛速度越慢。

(3)RLS自适应滤波算法

设t时刻的滤波器权系数向量W(t)为：

输入观测数据向量X(t)为：

设P(t)是p×p维矩阵，初始值为P(0)=α I，α是很大的标量，I是单位矩阵，那么，

RLS(Recursive Least Squares，递归最小二乘)自适应滤波算法为：

其中0<λ≤ 1是遗忘因子，λ越小收敛速度越快，λ越大收敛速度越慢，当λ=1时在经过一段时间后相当于恒权滤波。

(4)自适应滤波应用

(4.1)自适应AR参数估计

从上述自适应LMS和RLS算法中我们可以看出, 如果将滤波器输入改为x(t-1)，期望输出s(t)改成x(t), 则上述算法得到的权系数就相当于AR(p)模型的参数估计值, 即a(i)=-wi，i=1,2 ,…, p。

(4.2)周期信号的抵消或增强

假设随机信号x(t)是周期信号s(t)(如工频及其谐波干扰)和宽带信号w(t)之和，x(t) = s(t) + w(t)，由于宽带信号在延时足够长之后，信号的前后可视为不相关；而周期信号则依然相关。为了分离这两种信号，可把x(t)延迟M点之后s(t)=x(t-M)作为参考信号，这时输出y(t)为周期信号，则e(t)=x(t-M)-y(t)为宽带随机信号。

(4.3)自适应逆滤波

设我们所观察到的信号x(n)是由原始信号s(n)通过冲激响应为h(n)的系统产生的，即x(n)=s(n)*h(n)或X(z)=S(z)H(z)。我们希望在不知道系统冲激响应h(n)或传递函数H(z)的情况下，获取原始信号s(n)。

自适应逆滤波的原理如图所示，如果图中的d(i)不能得到，则可将已知的示标信号p(i)与原输入信号加在一起作为未知系统的输入，并同时将p(i)作为d(i)。输出y(i)就为原信号的估计值。自适应均衡和自适应解卷积问题都可归结为上述自适应逆滤波。

(5)FIR和IIR滤波器

数字滤波器(DF)是对数字信号实现滤波的线性时不变系统。数字滤波器分为低通、高通、带通和带阻滤波器，它们的理想频率响应曲线如图所示。数字滤波器分为有限冲激响应(FIR)滤波器和无限冲激响应(IIR)滤波器。

滤波器理想幅频响应特性，(a)低通 (b)高通 (c)带通 (d)带阻

(6)FIR滤波器设计

有限冲击响应(FIR)数字滤波器设计包括低通、带通、带阻、高通滤波器设计，这里采用窗函数方法，也称为傅里叶级数法。

(6.1)低通滤波器设计

假设低通滤波器的归一化截止频率为fc(fc=fa/fs，fa为实际截止频率，fs为信号的采样率)，滤波器的阶数为2M+1，则低通滤波器的系数B(i)为：

其中w(i)是窗函数，例如我们采用Hamming窗，则：

w(i)=0.54-0.46cos(π(M-i)/M)

(6.2)带通滤波器设计

设带通滤波器的通带为f1<f<f2 ，f为归一化频率(f=1对应采样率的一半)，则：

(6.3)高通滤波器设计

(6.4)带阻滤波器设计

(7)IIR滤波器设计

无限冲击响应(IIR)数字滤波器设计有多种方法，最常用的一种就是双线性变换法。首先设计出模拟低通滤波器，然后将其变换成低通、高通、带通和带阻滤波器。

(7.1)确定模拟低通滤波器的阶数

(7.2)模拟低通滤波器设计

模拟低通巴特沃兹(Butterworth)滤波器以Butterworth近似函数作为滤波器的系统函数，它的幅度平方函数表示为：

式中N是正数,表示滤波器的阶数，Ωc为通带的截止频率，从幅度平方函数中可确定出系统的函数Ha(s)：

其极点为：

这2N个极点在s平面上以象限对称，等间隔地分布在半径为Ωc的圆上，极点间的角度为π/N。当N为奇数时，在s平面实轴0与π之间处有极点；而当N为偶数时，则第一个极点必然出现在π/(2N)处。极点决不会出现在虚轴上，这样系统才有可能是稳定的。取左半平面上的极点作为Ha(s)的极点sk ，系统函数为：

(7.3)双线性变换

利用双线性变换从低通原型模拟滤波器到实际数字滤波器的变换(Ωp为实际截止频率，T为采样间隔)：

(8)梳状滤波器

梳状滤波器(Comb filter)的传递函数为：

式中参数：0<r<1，0<u<1，当r>u时为陷波器，当r<u时为增强器，滤波器阶数n与频率点相关。

如果陷波(增强)频率点为f的整数倍，采样率为Fs，那么n=Fs/f。

(9)Savitzky-Golay平滑

一种适合数据平滑的特殊低通滤波器，称之为萨维兹凯-戈雷(Savitzky-Golay)滤波器，或者最小二乘、数据平滑多项式滤波器。Savitzky-Golay滤波器更能保留信号相对极大值、极小值和宽度等特征。设滑动窗口的长度为n=2*m+1，对应的信号数据为x(-m), x(-m+1), …, x(0), …, x(m-1), x(m)，那么建立以下标i为自变量的k-1次多项式：

估计多项式系数aj，j=0,1, …, k-1，使得x(i)与y(i)之间的均方误差最小。

以y(0)=a0作为x(0)的平滑输出值。滑动窗继续向前滑动一点，再用新窗口内数据估计多项式系数，用同样的方法求下一点的平滑值，对于长数据该平滑器计算量偏大。

Savitzky-Golay滤波器参考文献：[1][美]W.H.Press, S.A.Teukolsky, W.T.Vetterling, B.P.Flannery 著. 傅祖芸，赵梅娜，丁岩 等译. C语言数值算法程序大全(第二版)[M]. 北京 电子工业出版社，1995年10月第1版. 第550-555页.

Numerical Recipes in C, The Art of Scientific Computing, Second Edition, Cambridge University Press, 1988, 1992. pp550-555.

[2]宁智 编写. Microsoft C 科学与工程工具库[M]. 北京 学苑出版社，1993年12越第1版，第44-46页.

[3]《数学手册》编写组. 数学手册[M]. 北京 高等教育出版社，1979年5月第1版，第903-911页.

参考文献

[1] 韩曾晋 编. 自适应控制系统[M]. 北京 机械工业出版社，1983年6月第1版.

[2] 何振亚 编. 数字信号处理的理论与应用（上、下）[M]. 北京 人民邮电出版社，1983年第1版.

[3] 张树京 编著. 统计信号处理[M]. 北京 人民邮电出版社，1986年第1版.

[4] 李素芝，万建伟 编著. 时域离散信号处理[M]. 北京 国防科技大学出版社，1994年8月第1版.

[5] 聂能，尧德中，谢正祥 主编. 生物医学信号数字处理技术及应用[M]. 北京 科学出版社，2005年8月第1版.

[6] 韩崇昭，王月娟，万百五 编著. 随机系统理论[M]. 西安交通大学出版社，1987年3月第1版.

(10)FIR、IIR滤波器设计举例

以对当前通道进行滤波为条件设计指定频带的滤波器，并且在滤波器幅频特性显示中以当前段信号的对数功率谱为背景。

(10.1)FIR滤波器设计与滤波。菜单操作：《处理》→ 《FIR滤波器设计》。出现如下对话框，输入FIR滤波器阶数，滤波器频带范围，可以设计滤波器并对当前通道信号进行滤波。还可以将设计好的滤波器系数保存到Fir_Coefs.txt中，便于之后使用，例如《处理》→《FIR外部定义滤波器》，将系数读入后对信号滤波。

(10.2)多通道等采样率下FIR滤波器设计与滤波。当通道数目大于1且各通道信号采样率相同时，以如下对话框来替代上述对话框，此时可用设计好的滤波器对所有通道进行滤波。

(10.3)IIR滤波器设计与滤波。菜单操作：《处理》→ 《IIR滤波器滤波》，出现如下对话框。IIR滤波器设计有可能出现系统不稳定情况。

(10.4)梳状滤波。菜单操作：《处理》→《梳状滤波》。以下是在特定频率点上的陷波器(左边)和增强器(右边)。对于陷波器来说有0<u<r<1，对于增强器来说有0<r<u<1。

(10.5)平滑与微分器。菜单操作：《处理》→《平滑与微分》。出现对话框，用户可以选择平滑或者微分器，观察滤波器的幅频特性，并对当前通道或者全部通道(在采样率相同的情况下)进行滤波。

点击“特性”按钮显示滤波器幅度频率响应曲线以及对于的FIR滤波器系数，在多通道同采样率的情况下勾上“全部通道处理”对所有通道全程实施滤波。

(10.6) Savitzky-Golay平滑。菜单操作：《处理》→《Savitzky-Golay平滑》，出现如下对话框输入参数：

输入滑动窗口长度之一半为5(窗长为2*5+1=11)，多项式阶数为4，对信号进行平滑。并且用11点平滑结果进行比较。如下图所示，黑色曲线为原始数据，红色为Savitzky-Golay平滑结果，蓝色为11点滑动平均。可以看出Savitzky-Golay平滑能够保持好极值的幅度。

(11)自适应实时滤波

(11.1)自适应信号增强或抵消

float RLSCalc(float xactual, float xref, float *w, int nw, float update, int init, short index)

RLSCalc函数采用递推最小二乘(RLS)算法做自适应实时信号抵消或增强，也可以自适应地估计自回归(AR)模型系数。

float LMSCalc(float xactual, float xref, float *w, int nw, float update, int init, short index)

LMSCalc函数采用递推最小均方误差(LMS)算法做信号抵消或增强，也可以自适应地估计AR模型系数。

函数说明：

假设t时刻观测到的数据x(t)是信号和噪声的叠加：x(t) = s(t) + e(t)，x(t)是观测数据，s(t)是信号，e(t)是噪声，与信号s(t)对应的参考输入信号为sr(t)；通过加权系数wi(t)，i=1,2, …, p，对t时刻的未知信号s(t)进行估计，得到估计值sg(t) = w1(t)*x(t) + w2(t)x(t-1) + … + wp(t)x(t-p+1)，此时信号的估计误差为es(t ) = sr(t) - sg (t)，采用RLS或者LMS算法用es(t)去不断修正权重系数wi(t)。此时函数调用对应的参数为观测数据输入xactual=x(t)，参考信号输入xref=sr(t)，此时函数返回值sg(t)是信号增强输出，eg(t)=x(t)-sg(t)是信号抵消输出。用于自适应AR(p)模型参数估计时，函数输入参数为xactual=x(t-1)，xref=x(t)，t时刻的AR模型系数ai=-wi(t)，AR模型为x(t) + a1x(t-1) + a2x(t-2) + … + apx(t-p) = w(t)。

参数：

float xactual：观测数据输入值；

float xref：参考信号输入值；

float *w：输入输出参数，加权向量，长度为p，w[i], i=0,1, …, p-1；

int p：权向量w[]的长度；

float update：更新速度，此值越大更新速度越快，越小更新越慢，0≤ update≤1；

init：操作选择器，=1对信号增强器进行初始化；=0执行信号增强操作；=2结束操作；

信号增强或者抵消的操作分三步：

①先对信号增强器做初始化(仅调用一次, init=1)；

②再对每一点输入数据求出增强结果(多次调用, init=0)；

③最后当不再对信号进行处理时关闭信号增强器(仅调用一次, init=2)。

short index：自适应滤波器选择，index=0-7，可同时使用8个信号自适应滤波器。

返回：函数返回滤波后信号增强的结果sg(t)(当init=0时)，eg(t)=x(t)- sg (t)是信号抵消结果。

k = 0, 1, …, 7。

(11.2) 信号通过线性系统

float SignalToLinearSys(float *A, int n, float *B, int m, float x, int mode, short index)

线性滤波系统的传递函数为：

如果x(t)为输入信号，y(t)为输出信号，滤波计算公式为：

参数：

float *A：传递函数分母系数，A[i], i = 0…n-1，A[0]=1；

int n：系数A[]的长度n，n大于或者1且小于127；

float *B：传递函数分子系数，B[i], i = 0…m-1；

int m：系数B[]的长度为m，m大于或者1且小于127；

float x：输入的当前时刻信号；

int mode：mode=0表示开始的几个数，初始化， mode=1输入的数据产生滤波后输出。

int index：0=< index<8，表示可以同时进行8个滤波。

返回：滤波器输出。

(12)周期信号增强或消减

设观测数据x(t)是周期信号s(t)和高斯白噪声w(t)之叠加：x(t)=s(t)+w(t)，s(t)=s(t+T)，这里T是周期。

通过叠加平均来消除噪声w(t)的影响得到s(t)的估计sa(t)，如果叠加的周期数目为m，则：

这里假设w(t)的均值为零。消减周期信号后得到的信号为：wa(t)=x(t)-sa(t)。在进行叠加消除干扰之前需要已知周期信号的周期T。

菜单操作：《处理》→ 《周期信号增强》，出现如下对话框，输入参数对当前通道信号进行处理。

信号中存在周期为0.2秒的方波，原始信号的采样率为100Hz，因此输入参数周期间隔点数为20，叠加次数为m=10个周期，得到如下增强后的周期信号和消减周期信号后的信号序列。

以下是一段信号，采样率为100Hz，从功率谱估计看信号中存在两个周期信号，一个是7Hz，周期为0.143秒，一个是10Hz，周期为0.1秒。叠加周期为m=10，取得如下结果以及对应的功率谱。

原始信号，增强后周期信号(周期为0.14s、0.1s)所对应的功率谱

一段实采信号增强其中10Hz成分，采样率为200Hz，叠加周期为10。

对应的功率谱

(13)自定义模板匹配滤波

在信号{x(t), t = 0, 1, …, n-1}中挑选一小段{s(i), i = 0, 1, …, m-1}作为模板，对信号进行模板匹配滤波。首先对模板进行归一化，设a是s(i)的平均值，归一化模板为{h(i), i = 0, i, …, m-1}，则

在进行滤波之前，先对信号x(t)去基线，去基线的滑动平均窗口长度为m。设匹配点为I0，则模板匹配滤波输出y(t)为：

匹配点I0通常是指要与信号的对准点，如h(I0)是模板的顶点，下图所示的模板波形。

模板匹配滤波的菜单操作：《处理》→ 《定义模板滤波》，此时菜单项《定义模板滤波》已经打上勾，表示此时在波形中拖动鼠标选择模板波形，选择模板方法是找到需要滤波通道中合适的模板波形，展开波形以观察细节，在要选取的波形起点按下鼠标左键不松开，往右拖动鼠标，抵达波形终点松开鼠标，即完成从信号中挑选模板。

松开鼠标后会弹出对话框，显示模板波形及参数(心电信号模板)：

在对话框中可以进一步调整模板参数，点击“滤波”按钮即开始对当前通道实施模板匹配滤波。以下是一段滤波结果。

再举一个例子，信号采样率为100Hz，选择模板波形及匹配滤波结果如下：

模板匹配滤波，选择的模板很重要，应该是需要提取的特征波形。

(14)余弦变换滤波

余弦变换信号滤波的方法是将一段信号进行余弦变换(采用快速算法)，将处于通带外的系数设置成零，之后再实施反余弦变换，得到滤波后的信号。因此可以根据需要在频域定义多个通带。可以使用余弦变换对信号全程或者当前信号进行滤波。当多通道信号的采样率相同时，可以选择对信号全部通道进行滤波处理。菜单操作：《处理》→ 《余弦变换滤波》，出现如下对话框。

在对话框中显示了当前段信号的余弦变换系数曲线，在系数显示图中可以拖放鼠标后弹出对话框来定义通带或者阻带频率范围。

一个例子：原信号的采样率为1000Hz，带通滤波频带为5-15Hz，滤波结果如下。余弦变换滤波是分段进行的，在输出信号段与段的连接处可能出现不连续突跳。

---

联系作者：chengbowork@163.com

相关文章

多通道信号分析软件系统
生物医学信号处理与分析软件系统设计
曲线拟合软件
离散小波变换用于信号滤波
多维特征参数机器学习算法
多维特征参数机器学习软件
Kohonen自组织特征映射神经网络(环形和球面形网络)
矩阵的三维图形显示软件
图片浏览软件工具
主成分分析(K-L变换)与信号的分解与合成(滤波)
信号的样本熵序列计算
信号的双谱分析
信号的经验模态分解(EMD)
希尔伯特(Hilbert)变换信号瞬时频率计算
给定概率分布的随机变量仿真
信号的特征参数计算
数据分布点纹图
多通道信号数据压缩存储
图像处理基本算法及软件简介
信号功率谱估计
特殊函数计算器