第5章
CHAPTER 5


随机信号的滤波





一般来说，确定性信号的滤波可以认为是对信号不同频谱分量的取舍，例如低通滤波器保留信号的低频分量，丢弃信号的高频分量。对于随机信号，滤波可以认为是对其功率谱不同分量的取舍，在时间域对应的操作就是对自相关函数的修改。随机信号滤波的本质是对其统计特征和概率特征的修改，而不是对波形的修改。因此，滤波器的设计准则也将基于相关统计特征，这和确定性信号滤波有显著的不同。本章在介绍滤波器基本概念的基础上，重点分析基于平方误差准则的滤波器及预测器设计，即
Wiener滤波、卡尔曼滤波、最小二乘滤波和线性预测； 然后介绍基于输出信噪比准则的滤波器
(匹配滤波器)设计； 最后介绍自适应滤波。本章只考虑数字随机信号的滤波。
5.1数字滤波器的基本概念
一个典型的数字滤波器框图如图51所示。



图51数字滤波器框图


设输入信号为x(n)，输出信号为y(n)，这里只考虑线性滤波，滤波过程可用以下差分方程来表示： 

y(n)=-∑M-1i=1aiy(n-i)+∑N-1i=0bix(n-i)(51)


由式(51)可知该数字滤波器的传递函数为

H(z)=∑N-1i=0biz-i1+∑M-1i=1aiz-i(52)

其单位冲激响应函数为

h(n)=z-1(H(z))(53)

y(n)=h(n)x(n)=∑∞i=0h(i)x(n-i)(54)

式中，n<0时，有h(n)=0，这样的滤波器系统称为因果系统。
如果冲激响应函数是有限长的，即

h(n)=hn,0≤n≤N

0,其他(55)

则称此滤波器为有限冲激响应(Finite Impulse Response，FIR)滤波器，否则，称为无限冲激响应(Infinite Impulse Response，IIR)滤波器。
如果h(n)满足如下条件： 

h(n)=0,n<0

∑∞n=0|h(n)|<C(56)

则称此滤波器是因果的，并且是稳定的。在这一章讨论稳定、因果的FIR或IIR滤波器。
随机信号滤波的一般问题模型如图52所示，其中x(n)是观测信号，d(n)是期望信号。目标是通过滤波器参数的合理设置，使得滤波器输出y(n)与期望信号尽可能接近，即使得误差信号e(n)=d(n)-y(n)尽量小。

图53给出了产生观测信号的两种模型，其中，图53(a)表示观测数据是由目标信号s(n)和噪声干扰w(n)组成的，图53(b)表示信号不仅受到加性噪声的干扰，还存在未知系统引起的畸变。假设s(n)与w(n)是广义平稳的随机信号，希望通过数据处理后，使噪声受到抑制，同时增强并恢复目标信号。根据期望信号的不同，一般可分为如下三种情况： 
(1) d(n)=s(n)，此时希望从x(n)中估计出s(n)，这是一个滤波问题； 
(2) d(n)=s(n+k)，k>0，此时希望估计x(n)将来的值，这是一个预测问题； 
(3) d(n)=s(n-k)，k>0，此时希望估计x(n)以前的值，这是一个平滑问题。



图52随机信号滤波的一般问题模型




图53产生观测信号的模型



5.2Wiener滤波
Wiener（维纳）滤波器的设计准则是： 考虑线性时不变滤波器，以使误差信号达到最小均方误差
，它是目前应用最广泛的滤波器。
5.2.1最小均方误差准则与正交性原理
设线性滤波器的冲激响应函数为h(n),n≥0，即

y(n)=∑∞k=0h(k)x(n-k)(57)

误差信号为

e(n)=d(n)-y(n)=d(n)-∑kh(k)x(n-k)(58)


最小均方误差(Minimum MeanSquared Error，MMSE)准则是使均方误差E［|e(n)|2］最小
，以此来设计滤波器。
构造如下目标函数： 

J=E［|e(n)|2］(59)

要使目标函数最小，应有Jh(k)=0,k=0,1,2,…。
由于

Jh(k)=E［|e(n)|2］h(k)

=Ee(n)［e(n)］h(k)

=Ee(n)d(n)-∑∞k=0h(k)x(n-k)h(k)

=-E［e(n)x(n-k)］(510)

Jh(k)=0，即

E［e(n)x(n-k)］=0，k=0,1,2,…(511)

这说明在MMSE准则下，误差e(n)与每个输入样本x(n-k)都是正交的，这就是所谓的正交性原理。式(511)和第3章中提到的正交性原理是一致的。
当滤波器的选择满足式(511)的正交性准则时，滤波器最优输出为

y(n)=∑∞k=0h(k)x(n-k)(512)

不难证明： E［e*(n)y(n)］=0，即输出信号与误差信号是正交的。
满足正交性准则的Wiener滤波器的均方误差为

Emin=E(|e(n)|2)

=Ed(n)-∑+∞k=0h(k)x(n-k)e(n)

=Ed(n)e(n)

=Ed(n)2-∑+∞k=0h(k)x(n-k)d(n)

=σ2d-∑+∞k=0h(k)Rdx(k)(513)


5.2.2WienerHopf正则方程
由正交性原理E［e(n)x(n-m)］=0可知

Ed(n)-∑∞k=0h(k)x(n-k)x(n-m)=0(514)

经简单计算，可得

Rdx(m)-∑∞k=0h(k)Rx(m-k)=0(515)

上式等价于

∑∞k=0h(k)Rx(m-k)=Rdx(m)(516)

此方程称为WienerHopf方程，又称为Wiener滤波器的正则方程。
对于FIR Wiener滤波器，有

∑Nk=0h(k)Rx(m-k)=Rdx(m)(517)

取m=0,1,2,…,N，可得如下方程组： 

Rx(0)Rx(-1)…Rx(-N)
Rx(1)Rx(0)…Rx(-N+1)
︙︙︙
Rx(N)Rx(N-1)…Rx(0)h(0)
h(1)
︙
h(N)=Rdx(0)
Rdx(1)
︙
Rdx(N)(518)

根据式(248)自相关矩阵定义，上式可简写为

RTxh= rdx


如果通过观测数据可得到自相关函数R^x(m)及互相关函数R^dx(m)的估值，则根据方程组(518)可求解出滤波器系数h=R-Txrdx。
下面考虑一种特殊情况，假设信号传输过程受到噪声干扰，观察信号为x(n)=s(n)+w(n)，滤波的目的是降低噪声的影响，考虑两种特殊情况： 
(1) 若d(n)=s(n)，且s(n)与w(n)互不相关，则

Rdx(m)=E［s(n)x(n-m)］

=Es(n)(s(n-m)+w(n-m))

=Rs(m)(519)

同时有

Rx(m)=Rs(m)+Ru(m)=Rs(m)+δ(m)σ2(520)

此时WienerHopf方程为

∑Nk=0h(k)Rx(m-k)=Rs(m)，m=0,1,2,…,N(521)


(2) 若d(n)=s(n+α),α>0，则Rdx(m)=Rs(m+α)。
此时WienerHopf方程为

∑Nk=0h(k)Rx(m-k)=Rs(m+α)，m=0,1,2,…,N(522)


例51已知x(n)=s(n)+w(n)，其中信号s(n)是AR(1)过程，s(n)=0.6s(n-1)+u(n)，u(n)是N(0,0.64)的白噪声，w(n)是N(0,1)标准高斯过程。试设计一个长度为2的
Wiener滤波器估计s(n)。
解s(n)的功率谱为

S(ω)=0.641-0.6e-jω2=0.641.36-1.2cosω(523)


对功率谱做傅里叶反变换，得到s(n)的相关函数Rs(m)=0.6m。根据信号模型可得

Rx(m)=Rs(m)+1·δ(m)= 0.6m+δ(m)(524)

由WienerHopf方程可得

2h(0)+0.6h(1)=1

0.6h(0)+2h(1)=0.6(525)

解得滤波器的系数为

h(0)=0.451

h(1)=0.165(526)


图54给出了例51中各信号的功率谱，可以看到式(526)所给出的滤波器能显著降低噪声的影响。


图54例51中滤波器的滤波效果



5.2.3Wiener滤波器的求解
求解Wiener滤波器的方法大致有如下3种： 
(1) 对于FIR滤波器，可直接用WienerHopf方程组求解。
(2) 对IIR滤波器，由WienerHopf方程可知

∑+∞k=0h(k)Rx(m-k)=Rdx(m)(527)


对式（527）进行z变换可得

H(z)Px(z)=Pdx(z)(528)

H(z)=Pdx(z)Px(z)(529)

因为Rx(m)，Rdx(m)是双边序列，这里需采用双边z变换。对式(529)采用双边逆z变换可求出滤波器系数h(k)，但由于采用双边变换，h(k)不是因果的，因此实用性较差，下面给出一种IIR滤波器的因果实现。


图55Wiener滤波器


(3) 将Wiener滤波器看成两个滤波器的级联，其中第一个滤波器为因果IIR系统，第二个滤波器为因果FIR系统，如图55所示。

图55中，第一个系统称为白化滤波器，输出i(n)是功率谱为σ2的白噪声，传递函数为G-1(z)； 第二个系统的传递函数是Q(z)。
图55中所示的Wiener滤波器的传递函数为H(z)=Q(z)G(z)，所以只需求出G-1(z)和Q(z)即可。
设x(n)的功率谱密度函数为Sx(z)|z=ejω，G(z)是可逆的最小相位系统，由X(z)=G(z)I(z)，可知

Sx(z)=σ2G(z)·G(z-1)(530)

假设Sx(z)能分解为Sx(z)=S+x(z)·S-x(z-1)，其中一部分零极点在单位圆内
，即|z|<1［记为S+x(z)］，而另一部分在单位圆外，即|z|>1［记为S-x(z)］，则

G(z)=S+x(z)(531)

实际上由于Sx(z)是有理谱密度，总可以分解为Sx(z)=S+x(z)·S-x(z-1)的形式，这就是所谓的谱因式分解定理。
考查图55中第二个滤波器，有

y(n)=∑∞k=0q(k)i(n-k)(532)

式(532)中，利用正交性原理，第二个滤波器的WienerHopf方程(参阅式（517）)为

∑∞k=0q(k)Ri(m-k)=Rdi(m),m=0,1,2,…(533)

由于i(n)是白噪声，所以Ri(m-k)=δ(m-k)σ2，代入式(532)，可得

q(m)σ2=Rdi(m)(534)

q(m)=1σ2Rdi(m),m=0,1,2,…(535)


设Rdi(m)的z变换为Γ(z)=∑∞k=-∞Rdi(k)z-k，由于Q(z)=∑∞k=0q(k)z-k，所以

Q(z)=1σ2Γ(z)+(536)

式中，Γ(z)+=∑∞k=0Rdi(k)z-k是Rdi(m)的单边z变换。
求解Rdi(m)过程如下： 

Rdi(m)=Ed(n)i(n-m)

=E∑∞k=0v(k)x(n-m-k)·d(n)

=∑∞k=0v(k)Rdx(m+k)(537)

其中v(k)是第一个滤波器的冲激响应函数，即

1G(z)=V(z)=∑∞k=0v(k)z-k(538)

对式(537)进行z变换

Γdi(z)=∑∞m=-∞∑∞k=0v(k)Rdx(m+k)·z-m

=∑∞k=0v(k)∑∞m=-∞Rdx(m+k)z-(m+k)zk

=∑∞k=0v(k)zk∑∞m=-∞Rdx(m+k)z-(m+k)

=∑∞k=0v(k)zk·Γdx(z)

=v((z)-1)Γdx(z)

=1G((z)-1)Γdx(z)(539)

所以

Q(z)=［Γdi(z)］+σ2=1σ2Γdx(z)G((z)-1)+(540)

所以最终的IIR Wiener滤波器的传递函数为

H(z)=Q(z)G(z)=1σ2G(z)Γdx(z)G((z)-1)+(541)

下面举例说明该方法的求解过程。
例52已知x(n)=s(n)+w(n)，其中信号s(n)是AR(1)过程，s(n)=0.8s(n-1)+u(n)，u(n)是N(0,0.36)的白噪声，w(n)是N(0,1)标准高斯过程。试设计一个因果IIR Wiener滤波器估计s(n)。
解s(n)是AR(1)过程，其自相关函数对应的z变换为

Ss(z)=0.36(1-0.8z-1)(1-0.8z)(542)

由题意可知

Sx(z)=Ss(z)+1=1.6(1-0.5z-1)(1-0.5z)(1-0.8z-1)(1-0.8z)(543)

可知σ2=1.6，且

G(z)=1-0.5z-11-0.8z-1(544)

互相关函数Rdx(m)=Rds(m)=Rss(m)，所以，Γdx(z)=Ss(z)。

Γdx(z)G(z-1)+=0.36(1-0.8z-1)(1-0.8z)1-0.5z1-0.8z+

=0.36(1-0.8z-1)(1-0.5z)+

=0.61-0.8z-1(545)

所以得到Wiener滤波器的传递函数为

H(z)=11.61G(z)Γxd(z)G(z-1)+

=11.61-0.8z-11-0.5z-10.61-0.8z-1

=3811-0.5z-1(546)

求z逆变换，得到滤波器的系数(单位冲激响应函数)为

h(n)=38·12n，n≥0(547)


图56给出了例52中滤波器的降噪滤波效果，可以看到式(547)所给出的一阶IIR滤波器性能优于图54中的FIR滤波器。


图56例52中滤波器的滤波效果



5.3线性预测
Wiener滤波是一个有广泛通用性的滤波理论。本节将用Wiener滤波理论来分析随机信号的可预测性问题，这里只考虑线性预测。线性预测的基本含义是用随机信号过去时刻的线性组合预测当前时刻的取值，线性预测是考察随机信号时间相关特性的重要手段，广泛应用于多个领域。

设x(n)在n时刻之前p个数据为x(n-p),x(n-p+1),…,x(n-1)，利用这p个数据
预测n时刻的值x(n)(如图57(a)所示)，称为“前向预测”，线性预测表示为

x^(n)=-∑pk=1akx(n-k)(548)



图57前向预测示意图


式(548)中，x^(n)为预测值，ak为预测系数。与5.1节所提到的预测问题相比，上述预测问题忽略了观测噪声的影响，可更集中考虑信号本身的特性。后续的推导也可拓展到噪声下的预测问题。式(548)称为无延时预测，当存在延时D时，线性预测表示为(如图57(b)所示)

x^(n)=-∑pk=1akx(n-D-k)(549)




有延时预测与无延时预测的预测器求解采用相同的方法，这里只讨论无延时预测。定义预测误差为

e(n)=x(n)- x^(n)(550)

下面基于最小均方误差准则推导预测器系数求解过程，预测均方误差为

ρ=E|e(n)|2=Ex(n)+∑pk=1akx(n-k)2(551)

欲使ρ最小，令ρak=0，可得

E［x(n-k)e(n)］=0，k=1,2，…,p(552)

式（552）和式(511)一致，是正交性原理在线性预测中的体现。要得到最佳预测器的求解，先对
式（552）求共轭，并用m代替k，可得

Ex(n-m)x(n)+∑Pk=1akx(n-k)=0，m=1,2，…,p(553)

可以得到如下包含预测器系数的方程： 

Rx(m)=-∑pk=1akRx(m-k)，m=1,2，…,p(554)

或以下矩阵形式： 

Rx(0)Rx(-1)…Rx(-p+1)
Rx(1)Rx(0)…Rx(-p+2)
︙︙︙
Rx(p-1)Rx(p-2)…Rx(0)a1
a2
︙
ap=-Rx(1)
Rx(2)
︙
Rx(p)(555)

此时最小的ρ值(最小预测均方误差)为

ρmin=E［x(n)x(n)-x^(n)］=Rx(0)+∑pk=1a*kRx(k)(556)


方程(554)与方程(555)称为线性预测WienerHopf方程，解此方程可以求解出线性预测模型的预测系数ak及最小预测均方误差ρ。
前面叙述的模型使用了p个n时刻之前的数据预测x(n)为“前向预测”。对平稳信号，时间相关性与时间起点无关。因此，亦可用n,n-1,…,n-p+1的数据向后预测x(n-p)，这称为“后向预测”。
为了表示清楚，把前向预测(forward prediction)模型记为

x^f(n)=-∑pk=1afkx(n-k)

ef(n)=x(n)- x^f(n)

ρf=E［|ef(n)|2］(557)


利用x(n),x(n-1),…,x(n-p+1)预测x(n-p)，称为后向预测(backward prediction)，如图58所示，其数学模型为

x^b(n-p)=-∑pk=1abkx(n-p+k)(558)




图58后向预测示意图


误差记为： eb(n)=x(n-p)-x^b(n-p)。
均方误差记为：  ρb=E［|eb(n)|2］。
利用正交性原理，令ρb最小，同样可得后向预测的WienerHopf方程

E［x(n-p+m)e(n)］=0，m=1,2，…,p(559)

Ex(n-p+m)x(n-p)+∑pk=1abkx(n-p+k)=0，m=1,2，…,p(560)

Rx(-m)=-∑pk=1abkRx(k-m)，m=1,2,…,p(561)

式（561）的矩阵表达形式为

Rx(0)Rx(1)…Rx(p-1)
Rx(-1)Rx(0)…Rx(p-2)
︙︙︙
Rx(-p+1)Rx(-p+2)…Rx(0)ab1
ab2
︙
abp=-Rx(-1)
Rx(-2)
︙
Rx(-p)(562)


式（562）和式(555)形式类似，利用自相关函数的共轭对称特性，对式(562)两边取共轭，对照式(555)，得到afk=［abk］，即前向预测系数和后向预测系数互为共轭。利用和前向预测类似的方法，可求得预测误差为

ρbmin=Rx(0)+∑pk=1［ab(k)］Rx(-k)(563)


注意到预测误差肯定是实数，对式（563）取共轭并利用预测系数关系和自相关函数特性，不难得出式（563）和式(556)是等价的，即当预测器长度相同时，前向预测和后向预测的误差相等。前向预测和后向预测的关系总结如下： 
(1) ρbmin=ρfmin； 
(2) 若afk、abk是实数，则afk=abk。
若afk、abk是复数，则afk=［abk］ (共轭关系)。

前向预测、 后向预测在变长度预测器的迭代求解中有重要应用，将在后面结合参数化功率谱估计深入分析。接下来考虑当随机信号可以用第4章线性模型描述时，其最佳线性预测的结果。首先考虑AR模型，一个p阶AR模型如下： 

x(n)=-∑pk=1φkx(n-k)+u(n)(564)

其中，u(n)是零均值白噪声，φk是模型系数。利用模型参数构造如下预测模型： 

x^(n)=-∑pk=1φkx(n-k)(565)

这时预测误差刚好为u(n)，因为u(n)和过去时刻的模型输出不相关，不难验证式（565）满足正交性原理，因此是最佳线性预测。
对比式(564)和式(565)，可看到p阶AR模型和p阶线性预测器是一致的，定义

A(z)=1+∑pk=1φkz-k(566)

AR模型和线性预测的关系如图59所示，由于滤波器A(z)能将x(n)变成一个白噪声u(n)
［预测误差信号e(n)］，所以称A(z)为白化滤波器，或反滤波器。AR模型、白化滤波器及线性预测器分别示于
图59(a)、图59(b)和图59(c)。


图59AR模型、白化滤波器和线性预测器



5.4卡尔曼滤波
在应用中，Wiener滤波理论必须把用到的全部数据存储起来，估计接收信号的自相关矩阵，然后再进行计算。按照这种滤波方法设置的专用计算机的存储量与计算量比较大，并且很难进行实时处理。在解决非平稳过程的滤波问题时，难以得到高效方法。到20世纪50年代中期，随着空间技术的发展，这种方法越来越不能满足实际应用的需要，面临着新的挑战。1960年和1961年，卡尔曼(R.E.Kalman)和布西(R.S.Bucy)提出了递推滤波算法，成功地将状态变量法引入滤波理论中，用消息与干扰的状态空间模型代替了通常用来描述它们的协方差函数，将状态空间描述与离散时间更新联系起来，适于计算机直接进行运算，而不用去寻求滤波器冲激响应的明确公式。这种方法得出的是表征状态估计值及其均方误差的微分方程，给出的是递推算法。这就是著名的卡尔曼理论，或称卡尔曼布西滤波。
卡尔曼滤波不要求保存过去的测量数据，当新的数据到来时，根据新的数据和前一时刻诸量的估值，借助于系统本身的状态转移方程，按照一套递推公式，即可算出新诸量的估值。与Wiener滤波器不同，卡尔曼滤波器能够利用先前的运算结果，再从当前数据提供的最新信息，即可得到当前的估值。卡尔曼递推算法大大减少了滤波装置的存储量和计算量，并且突破了平稳随机过程的限制，使卡尔曼滤波器适用于对时变信号的实时处理。
卡尔曼滤波基于信号的状态空间模型，如下式所示： 

xn+1=Fxn+Gun+ wn(567)

yn+1=Hxn+1+�瘙
經n+1(568)

其中，xn为p×1维系统状态向量，F为p×p维系统矩阵，un为k×1维输入控制向量，wn为系统动态噪声向量，G为系统控制矩阵，大小为p×k。式(567)表示系统特性可由若干参数形成的状态向量描述。当前的状态与前一个时刻的状态、当前输入和系统噪声有关，系统可能存在多路输入。
yn+1为观测向量，大小为r×1。H为观测矩阵，大小为r×p，描述系统状态对输出的影响。�瘙
經n+1为观测噪声向量，大小为r×1。
式(567)中系统状态随时间变化，因此式(567)称为状态方程，式(568)称为观测方程。
对于给定的系统，状态空间模型为已知，状态空间模型具有良好的一般性，对许多实际问题可用状态空间建模，实际上第4章的线性模型可以认为是状态空间模型的特例，例如例52中所描述的噪声中的信号AR(p)模型： y(n)=x(n)+v(n)。若定义系统的状态为xn=［xn,xn-1,…,xn-p+1］T(为描述方便，这里将x(n)简记为xn)，可用状态空间模型表示为

xn+1
xn
︙
xn-p+2=-φ1-φ2…-φp
10
0
10xn
xn-1
︙
xn-p+1+1
0
︙
0un+1(569)

yn+1
yn
︙
yn-p+2=1
1

1xn+1
xn
︙
xn-p+2+vn+1
vn
︙
vn-p+2(570)


下面讨论卡尔曼滤波的具体流程。卡尔曼滤波关注在系统矩阵、控制矩阵、观测矩阵已知的情况下根据观测数据对状态空间的估计。假设在时刻n，基于n时刻以前所获得的全部知识，对状态变量xn做出一个预测估计，记为x^n|n-1，则预测估计的误差为

en|n-1= xn- x^n|n-1(571)

称为预测误差。预测误差是零均值的，其协方差矩阵为

Cn|n-1=E［en|n-1eHn|n-1］=E［(xn- x^n|n-1)(xn- x^n|n-1)H］(572)

在预测估计x^n|n-1的基础上，利用n时刻所获取的新观测数据yn进一步改善对xn的估计，记为x^n，称为更新估计，更新估计通过下式完成

x^n= x^n|n-1+ Kn(yn-Hx^n|n-1)(573)

其中，Kn为待定的增益矩阵，称为卡尔曼增益，(yn-Hx^n|n-1)称为新息，代表由新观测数据yn获得的新信息。
更新估计的误差记作en，则有

en= xn- x^n= xn-［x^n|n-1+ Kn(yn-Hx^n|n-1)］(574)


卡尔曼滤波的实质是寻找适当的增益矩阵Kn，使更新估计的均方误差达到最小。更新估计的均方误差为

E［|en|2］=E［(xn- x^n)H(xn- x^n)］(575)

而更新估计的协方差矩阵为

Cn=E［eneHn］=E［(xn- x^n)(xn- x^n)H］(576)


比较以上两式可知，更新估计的均方误差E［|en|2］为协方差矩阵Cn的对角元素之和，因此，要使均方误差最小，就等效于使协方差矩阵对角元素之和最小。
将式(573)代入式(574)，可得

en=xn- x^n

=xn-［x^n|n-1+ Kn(yn-Hx^n|n-1)］

=xn- x^n|n-1- Kn(Hxn+ �瘙
經n-Hx^n|n-1)

=en|n-1- KnHen|n-1- Kn�瘙
經n

=(I- KnH)en|n-1- Kn�瘙
經n(577)

由式(576)得到

Cn=(I- KnH)Cn|n-1(I- KnH)H- KnRvKHn(578)

其中，Rv是观测噪声的自相关矩阵。注意式
（578）用到了信号与噪声的不相关特性。令

DnDHn=HCn|n-1HH- Rv(579)

将式(578)展开，并将DnDHn代入，得

Cn= Cn|n-1- KnHCn|n-1- Cn|n-1HHKHn+ KnDnDHnKHn(580)

可见，Cn中含有Kn的二次项和一次项，对其进行配方，得到

Cn= Cn|n-1+(KnDn-A)(KnDn-A)H-AAH(581)

其中，A为未知矩阵，式（581）展开为

Cn=Cn|n-1+ KnDnDHnKHn- KnDnAH-ADHnKHn+AAH-AAH

=Cn|n-1+ KnDnDHnKHn- KnDnAH-ADHnKHn(582)

比较式(580)和式(582)，可得

A= Cn|n-1HH(DHn)-1(583)


观察式(581)，其第一项和第三项均与增益矩阵Kn无关，只有第二项包含Kn。考虑到Cn及第二项为对称矩阵，故其主对角线元素是非负的，若要使Cn的主对角线元素之和达到最小，须要求第二项的主对角线元素之和达到最小，即为零。因此，有

KnDn-A=0(584)

则

Kn=AD-1n= Cn|n-1HH(DHn)-1D-1n

=Cn|n-1HH(HCn|n-1HH- Rv)-1(585)

Kn即为使更新估计均方误差达到最小的卡尔曼增益矩阵。
将式(585)代入式(578)，可得更新估计的协方差矩阵

Cn=(I- KnH)Cn|n-1(586)


由上述分析可知，从预测估计x^n|n-1和预测误差协方差矩阵Cn|n-1出发，即可求得卡尔曼增益矩阵Kn，然后由新的观测向量yn计算出新的信息(yn-Hx^n|n-1)，再加上x^n|n-1，就完成了对状态变量xn的更新估计x^n。
为了使上述估计过程能递推地进行，还需要由n时刻的估计更新x^n及其协方差矩阵Cn，来计算n+1时刻的预测估计x^n+1|n及其协方差矩阵Cn+1|n。有了x^n+1|n和Cn+1|n，就可以重复上述估计过程，进而得到x^n+1和Cn+1，完成n+1时刻的更新估计。
由状态方程

xn+1=Fxn+Gun+ wn(587)


若动态噪声wn是零均值的，并且不同n时刻的wn是不相关的，作为对xn+1的预测估计x^n+1|n，暂时忽略wn的作用，于是取

xn+1|n=Fx^n+Gun(588)

其中，输入控制向量un为已知，那么n+1时刻的预测误差为

en+1|n= xn+1- x^n+1|n=Fxn+Gun+wn-(Fx^n+Gun)=Fen+ wn(589)

注意到wn与en是互不相关的，那么n+1时刻预测误差协方差矩阵为

Cn+1|n=E［en+1|neHn+1|n］

=E［(Fen+ wn)(Fen+ wn)H］

=FCnFH+ Rw(590)

其中，Rw是wn的自相关矩阵。式
（590）可进一步用于n+1时刻卡尔曼增益矩阵的计算，从而形成迭代计算。根据最小均方误差准则，可以确定卡尔曼滤波器的初始值

x^-1=E［x-1］

C-1=E［(x-1- x^-1)(x-1- x^-1)H］(591)

式（591）需要已知模型状态的统计量信息。实际应用中当无法得到这些先验信息时，卡尔曼滤波可以采用更灵活的初始化手段，后续的迭代更新会优化估计结果。例如x^-1可以用全零或全1向量初始化，C-1可以用单位矩阵初始化。
至此，完成了向量卡尔曼滤波算法的推导，总结其递推过程如下： 
步骤1： 建立状态空间模型，如式(567)和式(568)所示。
步骤2： 设置初始化条件，即n=-1时，采用如式(591)所示的初始化，即

x^n=E［xn］


Cn=E［(xn- x^n)(xn- x^n)H］


步骤3： 预测，如式(588)所示

xn+1|n=Fx^n+Gun


步骤4： 计算预测误差的协方差矩阵，如式(590)所示

Cn+1|n=FCnFH+ Rw


步骤5： 计算卡尔曼增益如下［参见式(585)］：

Kn+1= Cn+1|nHH(HCn+1|nHH- Rv)-1


步骤6： 更新对状态的估计如下：

x^n+1= x^n+1|n+ Kn+1(yn+1-Hx^n+1|n)


步骤7： 估计误差的协方差矩阵如下［可参见式(586)］：

Cn+1=(I- Kn+1H)Cn+1|n


步骤8： 令n=n+1，重复步骤3~8，直到估计结束。
当待估计的系统状态向量恰好为一维时，即当待估计的状态变量只有一个时，向量卡尔曼滤波器便退化成为单变量卡尔曼滤波器，二者的递推过程完全一致。


图510卡尔曼滤波算法框图


根据上述递推过程，画出向量卡尔曼滤波算法的框图，如图510所示。

卡尔曼滤波的核心是对状态的预测和估计。当将卡尔曼滤波用于解决实际问题时，需要建立适当的状态空间模型，使得卡尔曼滤波模型中的状态为待估计的参数。例如，例51所示的AR模型滤波问题，采用式(569)和式(570)所描述的状态空间模型，简写为

xn+1=Fxn+Gun(592)

yn+1=Hxn+1+�瘙
經n+1(593)


基于卡尔曼滤波算法可得到x(n)的估计，卡尔曼滤波实现的功能是降噪滤波，并且可以证明，卡尔曼滤波的理论性能等价于因果IIR Wiener滤波器的性能。
在上述问题中，如果感兴趣的不是信号，而是模型的参数，则需修改状态空间模型。例如在无线通信的时变信道估计中，发送端发送已知的训练序列，接收端根据接收信号的训练序列估计信道参数，下面建立用于时变信道估计的状态空间模型。假设n时刻的信道参数表示为hn(l)，0≤l≤L-1，其中L表示信道响应的长度。进一步假设信道参数响应随时间的变换可以用一个AR(1)模型描述，即

hn(l)=alhn-1(l)+un,l，0≤l≤L-1(594)


定义系统状态向量为hn=［hn(0),hn(1),…,hn(L-1)］T，可得到如下状态方程： 

hn+1(0)
hn+1(1)
︙
hn+1(L-1)=a0
a1

aL-1hn(0)
hn(1)
︙
hn(L-1)+un，0
un，1
︙
un,L-1(595)

上式可推广到更一般的时变信道参数模型和时不变信道，请读者自行思考。接下来考虑观测方程，接收信号是信道参数和输入信号的卷积，即

yn=∑L-1l=0hn(l)xn-l+wn(596)

其中，wn表示信道加性噪声。定义输出向量为yn=［yn,yn-1,…,yn-p+1］T，可得到如下观测方程： 

yn+1
yn
︙
yn-p+2=xn+1xn…xn-L+2
xnxn-1…xn-L+1
︙︙︙
xn-p+2xn-p+1…xn-p-L+3hn+1(0)
hn+1(1)
︙
hn+1(L-1)+wn+1
wn
︙
wn-p+2(597)




与前述状态空间模型的最大不同是式(568)中的观测矩阵是时变的，但在所考虑的信道估计问题中，训练序列，即观测矩阵是已知的，其时变性不会影响卡尔曼滤波的处理流程。
卡尔曼滤波和Wiener滤波都基于最小均方误差准则，但两者有很大不同，主要体现如下。
(1) 卡尔曼滤波和Wiener滤波需要不同的统计量信息。由式(518)、式(529)所确定的Wiener滤波器需要计算目标信号和观测信号的互相关函数及观测信号的自相关函数，卡尔曼滤波只需系统噪声和观测噪声的统计量信息。卡尔曼滤波还需要初始化状态空间模型的系统状态，但如前面所提到，系统状态初始化可采用灵活的初始化策略。
(2) 卡尔曼滤波更适用于时变系统，这主要是因为Wiener滤波求解一个固定的滤波器并应用于所有情况，而卡尔曼滤波则采用每个符号更新的迭代计算。例如由式(595)、式(597)所描述的时变信道参数估计问题，难以用Wiener滤波直接求解。
5.5最小二乘滤波
前述Wiener滤波器需要接收信号相关二阶矩信息，卡尔曼滤波需预先知道输入控制信号和噪声的二阶矩。实际应用中这些统计信息往往无法预先得到，仅有部分观测信号及其对应的期望信号。一种解决方案是从可用的数据估计出需要的二阶矩，再采用前述的MMSE滤波器； 另一种是采用新的最小化性能标准，直接基于可用数据设计出最佳滤波器。本节将介绍第二种方法。
考虑式(57)所示的滤波问题，这里只考虑FIR滤波器，设线性滤波器的冲激响应函数为h(k),0≤k≤K-1，则滤波器输出信号可表示为

y(n)=∑K-1k=0h(k)x(n-k)(598)

令d(n)表示期望信号，误差信号为

e(n)=d(n)-y(n)=d(n)-∑K-1k=0h(k)x(n-k)(599)


采用基于3.7节的最小二乘准则求解滤波器系数，其目标函数是最小化每个时隙点误差模的平方和，如下所示： 

minh(k)∑N-1n=0|e(n)|2(5100)

其中，N代表观测信号及期望信号的长度。下面推导滤波器参数的求解算法。定义观测信号向量为

y= y(N-1),…,y(1),y(0)T(5101)


观测信号整体可用以下矩阵相乘的方式描述


y=y(N-1)
y(N-2)
︙
y(1)
y(0)=x(N-1)x(N-2)…x(N-K)
x(N-2)x(N-3)…x(N-K-1)
︙︙︙
x(1)x(0)…x(-K+2)
x(0)x(-1)…x(-K+1)h(0)
h(1)
︙
h(K-2)
h(K-1)(5102)



上式简记为

y=Xh(5103)

定义如下误差信号向量

e= e(N-1),…,e(1),e(0)T(5104)

误差信号能量可写为

E=eHe=(dH- hHXH)(d-Xh)

=dHd- hHXHd- dHXh+ hHXHXh(5105)

根据第3章的求导公式，可得

(dHd)h=0，(hHXHd)h=0

(dHXh)h=(XHd)，(hHXHXh)h= (XHX)Th

由E/h=0，可得

(XHX)h= XHd(5106)

h^= (XHX)-1XHd(5107)

最小二乘解和式(518)中的Wiener解有相同的形式。定义如下基于时间平均的自相关函数估计和基于时间平均的互相关估计

R^=1N∑N-1n=0x(n)xH(n)(5108)

c^=1N∑N-1n=0x*(n)d(n)(5109)

易知R^=(XHX)*/N，c^=XHd/N，利用共轭对称矩阵的特性，得到最小二乘解的另一种形式

h= R^-Tc^(5110)

因此，可以认为最小二乘解是利用时间平均代替统计平均的Wiener解。将式(5107)代入式(5105)可以得到最小二乘解的最优误差为

Emin=Ed- c^HR^-1c^(5111)

可以证明，这和Wiener解的形式也是相同的。当采用加权最小二乘法时，误差重写为

E= eHe=(dH- hHXH)W(d-Xh)(5112)

其中，W是加权矩阵，其作用是体现不同误差分离重要性的差异性。加权最小二乘解的推导可采用相同的方法。
5.6匹配滤波器
信噪比是雷达、声呐探测系统的基本指标，在雷达发展的初期，一般采用信噪比作为衡量雷达接收机抗干扰性能的指标。对于含有噪声的观测信号，接收滤波器的首要任务是降低噪声功率，即提高信噪比。1943年，North从最大信噪比准则出发，建立了匹配滤波器理论。简单而言，匹配滤波器就是最大化滤波器输出的信号噪声功率比，并且这里的信号功率一般指信号的瞬时功率。


图511连续时间线性系统


首先以连续时间系统阐述匹配滤波器的含义。这里只考虑实信号，如图511所示用于降噪的线性系统，输入为

x(t)=s(t)+w(t)(5113)


其中，s(t)为确定性目标信号，w(t)为零均值高斯白噪声，方差记为σ2。滤波器是线性系统，它满足迭加原理。因而，可得到滤波器的输出为

y(t)=s0(t)+w0(t)(5114)




注意到s0(t)是确定性信号，根据线性时不变系统特性，可表示为

s0(t)=12π∫+∞-∞S(Ω)H(Ω)ejΩtdΩ(5115)

w0(t)是随机信号，根据第4章所推导的线性时不变系统特性，w0(t)的功率谱为

Sw0(Ω)=Sw(Ω)|H(Ω)|2=σ2|H(Ω)|2(5116)

得到噪声的平均功率为

E［w20(t)］=12π∫+∞-∞σ2|H(Ω)|2dΩ(5117)

定义在某个时刻滤波器输出端信号的瞬时功率与噪声的平均功率比值为

Dt0=s20(t0)E［w20(t)］(5118)

把式(5115)、式(5117)代入，得到

Dt0=12π∫+∞-∞S(Ω)H(Ω)ejΩt0dΩ2∫+∞-∞σ2|H(Ω)|2dΩ(5119)

要求Dt0的最大值，相当于求上确界，可利用如下许瓦兹不等式

∫+∞-∞a(Ω)b(Ω)dΩ2≤∫+∞-∞a(Ω)2dΩ∫+∞-∞b(Ω)2dΩ(5120)

其中，a(Ω)，b(Ω)是积分有限的任意函数，上式中等号成立的充分必要条件是

a(Ω)=cb(Ω)(5121)

其中，c为常数。在式(5120)中，令

a(Ω)=σejΩt0H(Ω)，b(Ω)=S(Ω)/σ(5122)

则

∫+∞-∞S(Ω)H(Ω)ejΩt0dΩ2≤∫+∞-∞σ2H(Ω)2dΩ∫+∞-∞(S(Ω)2/σ2)dΩ(5123)

代入式(5119)，得到

Dt0≤12πσ2∫+∞-∞S(Ω)2dΩ(5124)

要使上式成立，只需令

H(Ω)=ce-jΩt0S*(Ω)/σ2(5125)




图512离散时间线性系统示意图



可以看到，这时滤波器是由目标信号完全确定的，因此h(t)称为匹配滤波器。接下来考虑离散时间线性系统，见图512。以从时间域推导匹配滤波器表达式的方式进行讨论。假设输入信号为有限长度序列s(n) (n=0,1,…,N-1)，离散时间线性滤波器的单位样值响应为h(n)，h(n)在区间s(n) (n=0,1,…,N-1)上取值，在该区间之外为零，噪声w(n)为零均值高斯白噪声，方差为σ2。

根据离散时间线性系统理论，输出信号s0(n)为

s0(n)=∑N-1k=0s(k)h(n-k)(5126)

输出噪声信号的平均功率为(参见习题4.1)

E(w0［n］2)=σ2∑N-1k=0h(k)2(5127)

任意取一时刻，不失一般性，取时刻N-1，瞬时功率为s0(N-1)2。定义离散时间线性滤波器的输出信噪比为

d0=s0(N-1)2E［w20(n)］(5128)

将式(5126)和式(5127)代入式(5128)，得

d0=∑N-1k=0s(k)h(N-1-k)2σ2∑N-1k=0h2(k)(5129)

令s=［s(0)s(1)…s(N-1)］T，h=［h(N-1)h(N-2)…h(0)］T，则

d0=1σ2hTs2hHh=1σ2(h)Hs2hHh(5130)

根据许瓦兹不等式可得

(h)Hs2≤((h)H(h))(sHs)=(hHh)(sHs)(5131)

当且仅当h=cs时，式(5131)等号成立，其中c为常数，因此

d0≤1σ2(hHh)(sHs)hHh=1σ2(sHs)(5132)

令E= sHs=∑N-1n=0s［n］2代表信号的能量，所以

d0≤1σ2(sHs)=Eσ2(5133)

当且仅当h=cs时式(5133)等号成立，c只影响滤波器的放大倍数，可令c=1，所以，当

h(N-1-n)=s(n)，n=0,1,…,N-1(5134)

时，输出信噪比达到最大，或等价于

h(n)=s(N-1-n)，n=0,1,…,N-1(5135)

输出的最大信噪比为

d0max=E/σ2(5136)


例53假定输入信号为

s(n)=1，0≤n≤3

0，其他(5137)

求匹配滤波器的单位样值响应和匹配滤波器的输出信号。
解信号长度N=4，则匹配滤波器的单位样值响应为

h(n)=s(N-1-n)=1，0≤n≤3

0，其他(5138)

输出信号为

s0(n)=∑N-1k=0s(k)h(n-k)=n+1，0≤n≤3
7-n，4≤n≤6
0，其他(5139)

输入信号、匹配滤波器的单位样值响应和输出信号的波形如图513所示。


图513匹配滤波器输出





5.7自适应滤波
5.7.1自适应滤波器的基本概念

Wiener滤波器、卡尔曼滤波器都是以已知信号和噪声的统计特征为基础，具有固定的滤波器系数。因此，仅当实际输入信号的统计特征与设计滤波器所依据的先验信息一致时，这类滤波器才是最佳的。否则，这类滤波器不能提供最佳性能。在实际中，往往难以预知这些统计特性，故实现不了真正的最佳滤波。在没有任何关于信号和噪声的先验知识的条件下，自适应滤波器利用前一时刻已获得的滤波器参数自动调节现时刻的滤波器参数，以适应信号和噪声未知或随机变化的统计特性，从而实现最优滤波。
自适应信号处理的研究工作始于20世纪中叶。美国通用电气公司研究出了简单的自适应滤波器，用以消除混杂在有用信号中的噪声和干扰。Widrow B.等于1967年提出的自适应滤波理论，可使自适应滤波系统的参数自动地调整达到最佳状况，而且在设计时，只需要很少的或根本不需要任何关于信号与噪声的先验统计知识，这种滤波器的实现差不多像Wiener滤波器那样简单，而滤波性能几乎与卡尔曼滤波器一样好。20世纪70年代中期，Widrow B.等提出自适应滤波器及其算法，发展了最佳滤波设计理论。
自适应信号处理的应用领域包括通信、雷达、声呐、地震学、导航系统、生物医学和工业控制等。自适应滤波器是相对固定滤波器而言的，固定滤波器属于经典滤波器，它滤波的频率是固定的,自适应滤波器滤波的频率则是自动适应输入信号而变化的，所以其适用范围更广。自适应滤波器的一般结构框图如图514所示，其中d(n)为滤波器输出期望信号，x(n)为滤波器输入信号，e(n)为滤波器输出和期望信号之间的误差信号。自适应滤波器以当前时刻的滤波器输出为基础，以误差最小化为准则，更新当前时刻的滤波器系数。
自适应滤波器根据应用环境的不同有不同的结构。下面结合实际应用介绍几种自适应滤波器的实现结构。图515给出了自适应噪声消除的原理图。S(n)为噪声干扰下的目标信号，n0,n1为互不相关的白噪声。自适应噪声消除的目的是将信号S(n)从噪声中分离出来，其基本原理是将参考噪声n1经过滤波器处理后接近实际噪声n0。



图514自适应滤波器的结构示意图




图515自适应噪声消除原理图



自适应滤波在生物医学信号处理中的一个典型应用
是“胎儿心电图自适应干扰清除”，如图516所示。其中目标信号——胎儿心电信号受到母亲心电信号的干扰，消除这种干扰的一种有效处理方式是将母亲胸部探头得到的心电信号(基本不包含胎儿心电信号)作为参考信号，经过自适应滤波器的滤波后接近胎儿心电信号中的干扰信号。


图516胎儿心电图自适应干扰消除


下面两小节介绍两种经典的自适应滤波算法，为了简化讨论，只考虑实信号。
5.7.2LMS自适应滤波器
本节介绍LMS（Least Mean Square，最小均方）自适应滤波器。
考虑阶数为N的FIR自适应滤波器

y(n)=∑Nk=0h(k)x(n-k)(5140)

为了下面推导方便，引入信号的向量描述如下： 

h=［h(0),h(1),…,h(N)］T

x(n)=［x(n),x(n-1),…,x(k-N)］T(5141)

误差信号为

e(n)=d(n)-y(n)

=d(n)- hTx(n)(5142)

均方误差为

ξ=E(|e(k)|2)

=E［d2(k)-2d(k)xT(k)h+ hTx(k)xT(k)h］

=E［d2(k)］+ hTE［x(k)xT(k)］hT-2E［d(k)xT(k)］h(5143)

记Rx=E(x(k)xT(k))为滤波器输入信号的自相关矩阵，rdx=E［d(k)xT(k)］为输入信号和期望信号的互相关向量，得到均方误差的简洁表达式为

ξ=E［|e(k)|2］= hTRxh-2rdxTh+Dd(5144)


当输入信号是广义平稳的随机信号时，由上式可见ξ是滤波器系数h的二次函数。我们称ξ(h)为误差性能函数，因为Rx是半正定矩阵，误差性能函数的曲面是一个向下凹的抛物面,典型的二维均方误差性能函数如图517所示。

从图517可以看到误差性能函数有唯一最小点，通过对滤波器系数求偏导数并令偏导数为零可求出滤波器系数，这也是Wiener滤波所采用的方法。在自适应滤波中，我们希望得到一种迭代求解的形式，使得滤波器系数可实时跟踪信号统计特性的变化。由于误差性能函数是二次函数，连续可导并且有唯一全局最小点，因此可以用降低的梯度下降法进行迭代。梯度下降法的示意图如图518所示，其基本思想是:要达到最优解，不管初始权值如何选择，只要在调整过程中，权值的调整沿ξ(h)的负梯度方向进行，就可以保证最终收敛到最优解。



图517二维均方误差性能函数




图518梯度下降法示意图



可见，梯度法可表示为如下的迭代公式： 

hn+1= hn+μ(-ξ(h)|h= hn)(5145)

这里引入时间n，hn表示时刻n的滤波器系数(第n次迭代求解的结果)，μ为迭代步长，一般是小于1的常数。式(5145)的梯度为

ξ(h)|h= hn=2Rxhn-2rdx(5146)


式(5146)需计算自相关矩阵，这种梯度计算方法不能应用于非平稳信号，和自适应滤波的初衷不一致。Widrow等提出了一种近似的计算梯度的有效方法，其原理是用当前时刻的误差|e(n)|2代替均方误差E(|e(n)|2)。此时有

ξ^(h)= e(n)2(5147)

e(n)=d(n)- hTx(n)(5148)

ξ^(h)|h= hn=-2e(n)x(n)(5149)

容易证明ξ^(h)|h=hn是ξ(h)|h=hn的无偏估计。此时梯度法迭代公式为

hn+1=hn+2μe(n)x(n)(5150)

此算法称为WidrowHopf  LMS算法，计算非常简单，并且不涉及信号的统计量，每个时刻的更新迭代只需当前接收信号及期望信号，符合自适应滤波实时更新信号特征的初衷。
LMS算法的实现步骤如下： 
步骤1： 初始化，选择合适的步长μ，初始化h0的值(一般为零或随机初始化)，设n=0； 
步骤2： 迭代更新，如式(5148)和式(5150)所示： 

e(n)=d(n)- hTnx(n)

hn+1= hn+2μe(n)x(n)


步骤3： 判断是否收敛，如果不收敛，令n=n+1，返回步骤2，一般如果滤波器在连续两次的更新小于某个门限值则认为算法收敛。
LMS算法简单，易于软硬件实现，应用非常广泛。从算法流程可以预见μ值的选择对算法有重要影响。μ越大，收敛越快，但容易振荡，μ值小，收敛慢，但比较平稳。
下面基于式(5146)分析μ值的选择对算法稳定性的影响。根据式(5146)把更新算法重写为

hn+1= hn-2μ(Rxhn- rdx)(5151)

注意到最优滤波器系数(记为hopt)满足WienerHolf方程Rxhopt=rdx，式
(5151)可改写为

hn+1=(I-2μRx)hn+2μRxhopt(5152)

其中，I是大小和Rx相等的单位矩阵。令�瘙
經n=hn-hopt，式(5152)可写为

�瘙
經n+1=(I-2μRx)�瘙
經n(5153)

记Rx的特征向量矩阵为Q，得到的特征值分解为

Rx=QΛQ-1(5154)

其中，Λ是包含Rx特征值的对角矩阵，由于实随机信号的自相关矩阵是对称矩阵，即Rx=QΛQ-1=RTx=Q-TΛQT，所以

Q-1= QT(5155)

式(5153)可重写为

�瘙
經n+1=(QQ-1-2μQΛQ-1)�瘙
經n

=Q(I-2μΛ)Q-1�瘙
經n(5156)

式(5156)左乘Q-1，并记v-n=Q-1�瘙
經n，得到

v-n+1=(I-2μΛ)v-n(5157)

递归应用式(5157)可得

v-n+1=(I-2μΛ)2v-n-1

=(1-2μλ1)n
(1-2μλ2)n

(1-2μλN-1)nv-0(5158)

要使算法稳定，应有

limn→∞v-n=0(5159)

这时有limn→∞Q-1�瘙
經n=0，即

limn→∞Q-1(hn- hopt)=0(5160)

即

limn→∞hn= hopt(5161)

其收敛到最优解。
要使式(5159)成立，必须有

limn→∞I-2μλin=0,i=0,1,2,…,N-1(5162)

即

|1-2μλmax|<1(5163)

所以

0<μ<1λmax(5164)

式中，λmax是Rx最大的特征值，这就是梯度法收敛的充分必要条件。一般而言，LMS算法收敛速度取决于μ和λmin，μ越大，收敛速度越快，λmin越大，收敛速度越快。通常定义d=λmax/λmin为谱动态范围，d越大，收敛时间越长。

5.7.3RLS自适应滤波
RLS(Recursive Least Squares，递推最小二乘法)自适应滤波器的设计准则是从滤波器开始运行到当前时间，滤波器系数的更新是令总的平方误差达到最小，令n表示当前时间，总的平方误差一般定义为

ξ(n)=∑nj=0λn-je(j)2=∑nj=0λn-jd(j)- hTnx(j)2(5165)

其中，e(j)是瞬时误差，λ(0<λ≤1)称为遗忘因子，用于降低“旧”数据的影响，提高滤波器的跟踪能力。
最优滤波器系数将使得ξ(n)最小，即满足以下方程

ξ(n)hn=0(5166)

基于式(5165)可得到包含最优滤波器系数的方程

R(n)hn=c(n)(5167)

其中

R(n)=∑nj=0λn-jx(j)xT(j)(5168)

c(n)=∑nj=0λn-jx(j)d(j)(5169)

如果R(n)是满秩矩阵，求解方程(5167)可得到滤波器系数，相应的解称为最小二乘解。但是直接求解上述方程需要求矩阵的逆矩阵。下面推导求解的迭代形式，首先R(n)的更新可以表示为

R(n)=λR(n-1)+x(n)xT(n)(5170)

c(n)=λc(n-1)+x(n)d(n)(5171)

将上述两式代入式(5167)，得到

［R(n)-x(n)xT(n)］hn-1=c(n)-x(n)d(n)(5172)

经过简单计算，得到

R(n)hn-1+x(n)［d(n)- hTn-1x(n)］=c(n)(5173)

注意d(n)-hTn-1x(n)刚好是应用前一时刻滤波器所造成的误差，记为d(n)-hTn-1x(n)=e(n-1)并在式(5173)两边乘以R-1(n)，得到

hn-1+ R-1(n)x(n)e(n-1)= R-1(n)c(n)= hn(5174)

定义如下自适应增益向量

g(n)= R-1(n)x(n)(5175)

最终得到

hn= hn-1+g(n)e(n-1)(5176)

式(5176)的更新依然涉及矩阵求逆，但其中的求逆可用如下迭代公式求解： 

A-xxT-1= A-1-A-1xxTA-T1+ xTA-1x(5177)

将R(n)=λR(n-1)+x(n)xT(n)代入式(5177)得到

［λR(n-1)-x(n)xT(n)］-1=λR-1(n-1)-λ-1R-1(n-1)x(n)xT(n)R-T(n-1)λ-11+λ-1xT(n)R-1(n-1)x(n)(5178)

定义

P(n)= R-1(n)(5179)

可得到如下迭代更新

P(n)=λP(n-1)-λ-1P(n-1)x(n)xT(n)PT(n-1)λ-11+λ-1xT(n)P(n-1)x(n)(5180)


最终得到最小二乘解的递归实现，称为递归最小二乘(Recursive Least Square，RLS)算法，具体步骤如下： 
步骤1： 初始化，设R(-1)=I，P(-1)=I，选择合适的遗忘因子λ初始化的值(一般为小于1但接近1的实数)，初始化滤波器系数，记为h-1，设n=0； 
步骤2： 迭代更新，如下式所示： 

e(n-1)=d(n)- hTn-1x(n)

P(n)=λP(n-1)-λ-1P(n-1)x(n)xT(n)PT(n-1)λ-11+λ-1xT(n)P(n-1)x(n)

g(n)=P(n)x(n)

hn= hn-1+g(n)e(n-1)


步骤3： 判断是否收敛，如果不收敛，令n=n+1，返回步骤2，收敛判断可以用和LMS算法相同的准则。
总之，RLS算法的计算复杂度高于LMS算法，但收敛速度较快。
本章习题
1. 设噪声中存在L个具有随机相位的复正弦信号

xn=∑Li=1Aiej(ωin+i)+vn

式中，i,i=1,2,…,L为均匀分布的随机相位，它们是互相独立的； vk为零均值与σ2v方差的白噪声，且与i互相独立。
(1) 证明E［ejie-ji］=δik,k=1,2,…,L； 
(2) 证明xn的自相关函数

rxx(k)=E［xnxn-k］=σ2vδ(k)+∑Li=1|Ai|2ejωik； 


(3) 将xn通过一个传递函数为A(z)=a0+a1z-1+…+aMz-M的滤波器滤波，滤波后的输出为yn，证明输出功率为

Py=E［y*nyn］=aHRxxa=σ2vaHa+∑Li=1Ai|2|A(ωi)|2

式中

a=［a0,a1,…,aM］T

Rxx=［rxx(i,j)］=［rxx(i-j)］

A(ωi)=∑Mm=0ame-jωim


2. 接收信号为x(t)=s(t)+v(t)，s(t)、v(t)均为零均值平稳信号，目标信号自相关函数为Rss(τ)=e-|τ|，噪声的自相关函数为Rvv(τ)=e-2|τ|，信号与噪声互不相关，请根据5.2.2节的思路导出因果连续IIR Wiener滤波器。
3. 设系统模型为xn=0.6xn-1+ωn，观测方程为zn=xn+vn，其中ωn为方差σ2ω=0.82的白噪声，vn为方差σ2n=1的白噪声，vn与xn不相关。试求其阶数为2的离散FIR Wiener滤波器并计算其均方误差。

4. 利用第3题的已知条件，计算因果IIR Wiener滤波器和非因果IIR Wiener滤波器，计算相应的均方误差并和第3题的FIR滤波结果进行比较。
5. 非因果FIR线性滤波器h(n)的输入/输出关系为

y(n)=∑Lk=-Lh(k)x(n-k)


设目标信号为d(n)，误差信号为e(n)=d(n)-y(n),试推导该非因果
FIR线性滤波器的WienerHopf正则方程。
6. 在以下AR(p)模型中

x(n)=-∑pk=1akx(n-k)+u(n)


u(n)是零均值，方差为σ2u的高斯白噪声。证明以下线性预测是线性最优预测并求
出预测误差。

x^(n)=-∑pk=1akx(n-k)


7. 证明阶数相同的线性前向预测系数和后向预测系数互为共轭。
8. 一维平稳随机信号x(n)的系统状态方程和观测方程分别为

s(n)=0.5s(n-1)+w(n)

x(n)=s(n)+v(n)

其中，w(n)、v(n)为白噪声，方差均为1。s(n)、w(n)、v(n)两两互不相关。试导出相应的卡尔曼迭代滤波公式，分析当k→∞时，预测参数的收敛值，并和最佳线性预测器进行比较。
9. x(n)的标量系统状态方程和观测方程定义为

s(n)=as(n-1)+w(n)

x(n)=s(n)+v(n)


写出从x(n)恢复出s(n)的卡尔曼滤波迭代计算过程，求出等效的系统传输函数（由x(n)到s(n)的传输函数）。
10. 对第9题中的标量状态空间模型，试分析其预测误差随时间增长的变化情况，是否存在随时间增长，预测误差逐渐下降的趋势？
11. 如下时变状态空间模型

xn+1= Fnxn+ wn

yn+1= Hn+1xn+1+ �瘙
經n+1


证明存在无观测噪声的等价状态空间模型，如下式所示

x′n+1=F′nx′n+w′n

yn+1=H′n+1x′n+1


12. 基于式(5105)的加权最小二乘滤波误差为

E= eHe=(dH- hHXH)W(d-Xh)


试分别针对实数模型和复数模型，推导最佳滤波系数及其对应的误差。
13. 线性滤波器的输入为x(t)=s(t)+v(t)，v(t)为高斯白噪声，功率谱密度为N0，s(t)为方波信号，定义如下： 

s(t)=A，0≤t≤A

0，其他


设计匹配滤波器并求最大输出信噪比。当采用如下滤波器进行降噪时，输出最大信噪比是多大？要使信噪比尽量大，如何选取参数a？

h(t)=e-at，0≤t≤A

0，其他


14. 考虑随机信号的匹配滤波问题，假设接收信号模型为

x(n)=s(n)+v(n)


其中，s(n)为平稳随机信号，w(n)为平稳噪声，s(n)与w(n)相互独立且其自相关矩阵已知，分别记为Rss和Rvv。考虑线性因果FIR滤波器h(n)，0≤n≤L。证明滤波器输出可表示为

y(n)= hHx(n)


其中，h=［h(0),…,h(L)］T，x(n)=［x(n),…,h(n-L)］T，证明滤波器输出信号的信噪比可表示为

SNR=hHRsshhHRvvh


15. 利用第14题的结论，证明当噪声是白噪声时，使输出信噪比达到最大的匹配滤波器为Rss最大特征值对应的特征向量。以上结论是否可拓展到有色噪声？
16. 考虑如图519所示单权自适应线性组合器。设开关S是断开的，E［x2(n)］=1，E［x(n)x(n-1)］=0.5，E［d2(n)］=4，E［d(n)x(n)］=-1，E［d(n)x(n-1)］=1，试导出性能函数表达式，并给出性能函数的图形。

17. 如图519所示，设开关S闭合，按第16题的要求再做一次。当只有一个权因子时，
性能函数表达式表示什么类型的曲线？
18. 考虑如图520所示自适应噪声对消系统，试：



图519单权自适应线性组合器




图520自适应噪声对消系统



(1) 写出性能表面表达式； 
(2) 确定自适应增益的范围； 
(3) 写出这种情况下的LMS算法。

19. 给定如图521所示系统辨识结构，试给出递归LMS算法。

20. 对于如图522所示逆模拟系统，设对所有i≠n，有s(n)与s(i)相互独立，且Rss(0)=1，同时设n0(n)与s(n)是相互独立的白噪声，且Rnn(0)=Pn，试推导以下功率谱的表达式： Sss(z)、Sxx(z)、Sdd(z)、Sdx(z)。



图521系统辨识结构




图522逆模拟系统



21. 信号系统结构
如图523所示，有一部分信号泄漏到参考信道，设n0(n)是总功率为N的白噪声(独立)。试用输入功率谱Sss(z)表示最佳Wn(z)。


图523信号系统结构