- 微信
- 微博
  
  分享文章到微博
- 复制链接
  
  复制链接到剪贴板

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

aqhs 发表于 2022/07/10 17:02:21 2022/07/10

【摘要】滤波器是一个系统，它对输入系统的信号起到滤波的作用。滤波器设计是建立一个信号系统，使得信号通过该系统处理后可以得到所期望的输出。滤波的目的不同导致信号滤波方法多样。自适应滤波器是指根据环境的变化，使用自适应算法来改变滤波器的参数和结构的滤波器，LMS和RLS自适应滤波算法就是根据估计误差实时调节线性滤波器参数的算法。简述了FIR和IIR数字滤波器设计、周期信号增强或抵消、自定义模板匹配滤波等。

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

滤波器是一个系统，它对输入系统的信号起到滤波的作用。滤波器设计是建立一个信号系统，使得信号通过该系统处理后可以得到所期望的输出。滤波的目的不同导致信号滤波处理的方法各种各样。自适应滤波器是指根据环境的变化，使用自适应算法来改变滤波器的参数和结构的滤波器，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)余弦变换滤波

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