现代信号分析
以步长为0.005为例）：
load('ecgca102_edfm.mat'); E=val(1,:)./10000; E1=val(3,:)./10000; thorax=E(1:3000); %胸部信号，参考信号 dn=thorax.'; abdomen=E1(1:3000);% 腹部信号，输入信号 xn=abdomen.'; M=20; %阶数 itr=3000;% 迭代次数 en=zeros(itr,1); %误差序列 W=zeros(M,itr);% 权系数 mu=0.5;%步长因子 s=zeros(itr,1);%均方误差 %迭代计算 for k=M:itr x=xn(k:-1:k-M+1); %滤波器M个抽头的输入 y=W(:,k-1).'*x; %滤波器输出 en(k)=dn(k)-y; %第k次迭代的误差 W(:,k)=W(:,k-1)+2*mu*en(k)*x; %权系数更新 s(k)=s(k-1)+abs(en(k)).^2; c(k)=sqrt(s(k)/k);%均方误差 end %求最优时滤波器的输出序列 yn=inf*ones(size(xn)); for k=M:length(xn) x=xn(k:-1:k-M+1); yn(k)=W(:,end).'*x; end figure; subplot(311),plot((1:3000),xn-yn,'r',(1:3000),dn,'b'); legend('胎儿心电图','母亲心电图'); xlabel('采样点数');ylabel('幅度'); subplot(312);plot(c); xlabel('迭代次数');ylabel('均方误差'); subplot(313);plot(W(1,:),'r');hold on; plot(W(2,:));plot(W(3,:),'k'); xlabel('迭代次数');ylabel('滤波器系数'); legend('权系数1','权系数2','权系数3'); 修改步长参数，即mu，分别取0.005，0.05，0.5。运行程序，得到的仿真图形分别如图4-6、图4-7、图4-8所示。
图4-6 步长为0.005时的LMS算法 图4-7 步长为0.05时的LMS算法
- 21 -
图4-8 步长为0.5时的LMS算法
比较上面3幅图，可以很明显的发现，随着步长因子的增加，LMS算法的收敛速度，失调量也随着变化，当步长因子?增大，系统的收敛速度变快，但是相应的失调量也越来越大（如图4.7所示），当步长因子?大到某个数值时，系统不再收敛，也就是系统不再稳定（如图4.8所示）；但是这并不意味着?越小越好，因为步长因子?越小，则系统收敛速度变慢，这在实际应用中也是不容许的。因此参数?的选择应从整个系统要求出发,在满足精度要求的前提下,尽量减少自适应时间。
3.2 用RLS算法得到胎儿心电图
递归最小二乘（RLS）算法，其实现步骤有点类似卡尔曼滤波器，也是通过几个递归公式来循环。采样点数3000，滤波器阶数设置为80，初始化条件设置为：R?1(0)??I,?是小的正数，取0.001。
w(0)?wI?0。选择改变的参数为遗忘因子?，分别取0.97、0.7、0.5，观察滤波效果。同样，由
于滤波器系数过多，这里只观测前三个滤波器系数与迭代次数的曲线。
RLS算法得到胎儿心电图的MATLAB源代码如下（此处以遗忘因子?=0.97为例）： %迭代计算 for k=M:itr x=xn(k:-1:k-M+1);%滤波器M个抽头的输入 K(:,k)=(Rxx*x)./(la+x.'*Rxx*x);%更新滤波器增益 e(k)=dn(k)-W(:,k-1).'*x;%误差序列 W(:,k)=W(:,k-1)+K(:,k)*e(k);%更新权系数 Rxx=(Rxx-K(:,k)*x.'*Rxx)/la;%更新自相关矩阵逆更新 y(k)=W(:,k).'*x;%滤波器输出 e1(k)=dn(k)-y(k); s(k)=s(k-1)+abs(e1(k)).^2; c(k)=sqrt(s(k)/k);%均方误差 end figure; subplot(311), plot((1:3000),xn-y,'r',(1:3000),dn,'b'); legend('胎儿心电图','母亲心电图'); xlabel('采样点数');ylabel('幅度'); subplot(312);plot(c); xlabel('迭代次数');ylabel('均方误差'); subplot(313);plot(W(1,:),'r');hold on; plot(W(2,:));plot(W(3,:),'k'); xlabel('迭代次数');ylabel('滤波器系数'); legend('权系数1','权系数2','权系数3');
- 22 -
现代信号分析
仿真波形分别如图4-9、图4-10、图4-11所示。
图4-9 遗忘因子为0.97时的RLS算法 图4-10 遗忘因子为0.7时的RLS算法
图4-11 遗忘因子为0.5时的RLS算法
比较上面3幅图，根据均方误差与迭代次数的关系曲线图可以发现当?=0.9时，系统在1800点处收敛；而当?=0.5时，系统在500点就收敛了，所以得出结论：遗忘因子?越小，系统的跟踪能力越强，但同时可以发现其对噪声也越敏感；相反，遗忘因子?值越大，系统跟踪能力减弱，但对噪声不敏感，收敛时估计误差也越小。所以，在用普通递推RLS算法时，一定要对?有个准确的取值，才能保证系统性能的最佳状态，所以一般文献要求?的取值范围为0。 .95???0.995
五、维纳滤波器
题目：观测信号为y(k)?s(k)?e(k),其中有信号是s(k)为恒量平稳序列，其统计特征已求得为
?E[s(k)]?0?E[s2(k)]??2s?
噪声e(k)是零均值白噪声，且与有用信号不相关，即
- 23 -
?E[e(k)]?02 ??E[e(i)e(j)]??e?(i?j)
??E[e(k)s(k)]?0求维纳滤波器？（10分）
维纳滤波器实际上就是在最小均方误差条件下探确定滤波器的冲击响应h(n)或系统函数H(z)，也就是求解维纳-霍夫方程的问题。因噪声e(k)是零均值白噪声，且与有用信号不相关,根据已知条件
?E[s(k)]?0?E[s2(k)]??2
s??E[e(k)]?0?2 及 ?E[e(i)e(j)]??e?(i?j)
??E[e(k)s(k)]?0 可求得
rsy(n)?E?s(k)y?(k)?
?E?s(k)s?(k?n)??E?s(k)e?(k?n)??rs(n)??s2ry(n)?E?y(k?n)y?(k)?（s(k)为恒量平稳信号）
?E?s(k?n)?e?(k?n)??E?s(k)?e?(k)??rs(n)?re(n)??s2??e2?(i?j) 则维纳霍夫方程为?Rs?Re?h?rd。
以一阶维纳滤波器为例，则维纳霍夫方程可写为
??s2??e2?s2??h(0)???s2???2? ??222???s??e??h(1)???s???s 求解该方程得
?h(0)?1??h(1)?2?2??2??se 维纳滤波器为：
??s2??2? ??s?H(z)? 均方误差
?s22???2s2e(1?z?1)
- 24 -
现代信号分析
?(k) ??Es(k)?s?2??s2?e2 ?rs(k)??h(k)r(k)?222?s??ek?01?sy
六、FIR维纳滤波器
题目：假定下图所需的信号的
d(n)是一个正弦序列
d(n)=
sin(n?0??),?0在0~?之间任取一值，噪声序列?1?n?和?2?n?都是AR(1)过程，分别由如下
的一阶差分方程产生：
v1(n)??v1(n?1)?e(n)v2(n)??0.6v1(n?1)?e(n)
式中0???1,e(n)是零均值，单位方差的白噪声，与d(n)不相关，取不同的?值，如0.1，0.9，试采,2个不同阶数的FIR维纳滤波器对信号x(n)?d(n)??1(n)进行滤波。并分析其结果。(12分)
维纳滤波器（Wiener filter）是由数学家维纳（Rorbert Wiener）提出的一种以最小平方为最优准则的线性滤波器。其数字系统描述如下图所示，它是要从噪声观测x(n)=d(n)+v(n)中恢复信号d(n)。
d(n)δ(n)W(z)^d(n)+-+e(n)
图6-1 维纳滤波器的一般框图
题目中所示系统的Wiener-Hopf方程可以推到如下：
将v2(n)作为维纳滤波器的输入，以估计噪声v1(n)，则Wiener-Hopf方程为：其中Rv2Rv2??rv1v2，是v2(n)的自相关矩阵，rv1v2是期望信号v1(n)和滤波器输入v2(n)之间的互相关矢量。对互相关量，可推导
- 25 -