数字信号处理实验仿真结果

西安电子科技大学统计与自适应信号处理实验仿真题目基于LMS算法的自适应滤波器的matlab仿真学院电子工程学院专业电路与系统学号学生姓名授课教师撰写日期： 2012年 12 月 20 日在自适应滤波器中，参数可调的数字滤波器一般为FIR 数字滤波器，IIR 数字滤波器或格型滤波器。

图1中，()x n 表示时刻n 的输入信号，()y n 表示时刻n 的输出信号，()d n 表示时刻n 的信号或期望响应信号，()e n 表示时刻n 的误差信号。

误差信号为期望响应信号()d n 与输出信号()y n 之差，记为()()()e n d n y n =-。

自适应滤波器的系统参数受误差信号控制，并根据()e n 的值而自动调整，使之适合下一时刻(1)n +的输入(1)x n +，以使输出信号(1)y n +更加接近期望信号(1)d n +，并使误差信号(1)e n +进一步减小。

当均方误差2[()]E e n 达到最小值时，()y n 最佳地逼近()d n ，系统已经适应了外界环境。

2.2 LMS 算法1）2[()]E e n 与权值W 的关系LMS 自适应滤波器通过算法，当2[()]E e n 最小时，滤波器已经调节出适合外部环境的滤波器权值W 。

我们可以先推导2[()]E e n 与加权系数W 的关系式。

写成矩阵形式：10()()[()][][][()]N T T i i i y j W x j X j W W X j -====∑（1）误差：()()()()[][()]T e j d j y j d j W X j =-=-（2）则：222[()][()[][()]] [()]2[()[()]][][[][()][()][]]T T T T E e n E d j W X j E d j E d j X j W E W X j X j W =-=-+（3）令2[][()[()]],[][[()][()]],[()](0)T dd P E d j X j R E X j X j E d j ϕ===代入式（3），则有：22111[()][()]2[][][][][] (0)(0)2(0)i m i T T NNNdd i m x x i x d i m i E e j E d j P W W R W WW W ϕϕϕ====-+=+-∑∑∑（4）可以从上式看出均方误差2[()]E e n 是加权系数W 的二次函数，它是一个中间上凹的超抛物线曲面，是具有唯一最小值的函数，即2[()]E e n 与W 的关系在几何上是一个“碗形”的多维曲面。

为了简单，设W 是唯一的，则2[()]E e n 与W 的关系成为一个抛物线。

调节加权系数W 使(1)()()W j W j j μ+=-∇（5）其中()j ∇为：2122[]()2()2[]e j Gw j P e ωω⎡⎤∂⎢⎥∂⎢⎥∇==-⎢⎥∂⎢⎥∂⎢⎥⎣⎦（6）为控制收敛速度与稳定性的数量常数，称为收敛因子或自适应常数。

式（5）中第二项前的负号表示当梯度值为正时，则权系数应该小，以使2[()]E e n 下降。

根据式（5）的递推算法，当权系数达到稳定时，一定有()0j ∇= ，即均方误差达到极小，这时权系数一定达到所要求的最佳权系数W Λ。

LMS 算法有两个关键：梯度()j ∇以及收敛因子μ的选择。

按（6）计算()j ∇时，要用到统计量G ，P ，因此有很大困难，故通常用一种粗糙，但却有效的方法，就是()j ∇用ˆ()j ∇代替，即 2112ˆ()2l l e e j e e ωωωω⎡⎤∂∂⎡⎤⎢⎥⎢⎥∂∂⎢⎥⎢⎥⎢⎥⎢⎥∇==⎢⎥⎢⎥∂∂⎢⎥⎢⎥⎢⎥⎢⎥∂∂⎣⎦⎣⎦（7）式（7）的含义是指单个误差样本的平方 作为均方误差 的估计值，从而使计算量大大减少。

从而最终可以推出权系数迭代的LMS 算法为：(1)()2()()W j W j e j X j μ+=-（8）w，根据上式可以逐步递推得到最佳权X j为输入样本向量，只要给定系数迭代的初值(0)()3.1本例通过设计一个二阶加权系数自适应横向FIR滤波器，对一正弦信号加噪声信号进行滤波。

为了实现该功能，先生成一个标准正弦信号()n n，s n和一个随机噪声信号()然后将()x n，再依照由LMS算法推导出来n n相加就得到了加噪后的正弦信号()s n和()的公式（8），设计自适应滤波器算法，对噪声干扰信号进行滤波，最后得到滤波后的信号e n，仿真实验结果如下所示：()图2 μ取0.00026时的滤波效果图当μ取0.00026时得到的效果较好。

前面一段时间较模糊是因为滤波器参数还没有调整到最佳，如图3所示，当t取0.5时，已经找到了最佳权系数。

图3 μ取0.000026时的滤波效果图当μ取0.000026时，滤波结果几乎呈直线，而且线条很粗，说明寻找加权系数的速度太慢了，如图3所示。

图4 μ取0. 26时的滤波效果图当μ取0. 26时，滤波结果也是呈直线状，而且线条很细，有的地方还有毛刺，说明系数参数变化太快，系统还没有调整到最佳加权系数，如图4所示。

图5 μ取1时的滤波效果图当μ取1时，系统输出混乱，如图5所示。

实验结果表明：不同的μ值得到的滤波效果是不同的。

通过实验数据观察得出：μ偏大时，自适应时间越短，自适应过程越快，但它引起的失调也越大，所以导致滤波结果很模糊，输出信号变化较大，当μ大于某个值时，系统输出混乱；μ偏小时，系统比较稳定，输出信号变化小，失调也小，但自适应过程却相应加长的，因此参数μ的选择应从整个系统要求出发，在满足精度要求的前提下，尽量减少自适应时间。

3.2自适应横向滤波器有2个权值，输入随机信号v（n），样本间相互独立，一般用高斯白噪声代替，且功率为0.02，即E((r(n))^2)=0.02，信号周期N＝16个样点，利用Matlab：（1）给出X（n）的图形;（2）编程，分别给出W(0)=[0,0]',u=0.1及W(0)=[4,-10],u=0.05时; a、LMS算法的权值随时间n变换的轨迹：W1(n)～n；W2(n)～n ;b、误差e(n)随时间变化的关系曲线：e(n)～n.实验仿真结果如下所示：图6 实验仿真结果如图6所示，对于输入的随机信号，使用参数可调的数字滤波器对信号进行滤波，此时的滤波器有两个权值，随着迭代次数的增多，滤波器的权值的变化值越来越小，信号的误差也趋于稳定，因此就得到了输出的最佳值。

3.3白噪声经过AR模型的输出作为LMS滤波器的输入，已知：a1=1.558；a2=-0.81；白噪声方差为1.0，均值为0；u＝0.002；利用Matlab实现：（1）给出单次运算和200次运算的LMS算法学习曲线（2）给出单次运算和200次运算的权值随n变换曲线实验结果如下所示：图7 实验仿真结果实验仿真结果如图7所示，随着迭代次数的增多，滤波器的权系数趋于恒定，随机信号的均方误差的变化也越来越小，当均方误差的值最小时，就得到了我们需要的最佳滤波器权系数。

不同的μ值得到的滤波效果是不同的。

通过实验数据观察得出：μ偏大时，自适应时间越短，自适应过程越快，但它引起的失调也越大，所以导致滤波结果很模糊，输出信号变化较大，μ偏小时，系统比较稳定，输出信号变化小，失调也小，但自适应过程却相应加长的，因此参数μ的选择应从整个系统要求出发，在满足精度要求的前提下，尽量减少自适应时间。

四．实验仿真代码4.1 实验仿真结果1：clc;clear all;close all;t=0:1/10000:1-0.0001;s=sin(2*pi*t);n=randn(size(t));x=s+n;w=[0,1];u=0.00026;for i=1:9999y(i+1)=n(i:i+1)*w';e(i+1)=x(i+1)-y(i+1);w=w+2*u*e(i+1)*n(i:i+1);endfigure(1)subplot(4,1,1)plot(t,n);title('Noise signal');xlabel('t');ylabel('n(t)');subplot(4,1,2)plot(t,s);title('Sinusoidal signal');xlabel('t');ylabel('s(t)');subplot(4,1,3)plot(t,x);title('Sinusoidal signal with noise');xlabel('t');ylabel('x(t)'); subplot(4,1,4)plot(t,e);title('the result of filtering');xlabel('t');ylabel('e(t)');4.2实验仿真结果1：clc;clear all;close all;c1=zeros(1,1024);c2=zeros(1,1024);a1 = -0.1950;a2 = 0.95;sigma = 0.0965;u = 0.04;% w=wgn(1,150,sigma);w=sigma*randn(1,1024);for n=1:1024if (n==1)x(n)=w(n);e(n)=x(n);c1(n)=0;c2(n)=0;elseif (n==2)x(n)=w(n)-a1*x(n-1);e(n)=x(n)-c1(n-1)*x(n-1);c1(n)=c1(n-1)+2*u*e(n)*x(n-1); c2(n)=c2(n-1);elsex(n)=w(n)-a1*x(n-1)-a2*x(n-2);e(n)=x(n)-c1(n-1)*x(n-1)-c2(n-2)*x(n-2);c1(n)=c1(n-1)+2*u*e(n)*x(n-1);c2(n)=c2(n-1)+2*u*e(n)*x(n-2);endendfigure(1)；plot(x,'k');figure(2)；plot(c1,'k');hold on;plot(c2,'r--');hold off figure(3)；plot(e);4.3实验仿真结果1：clc;clear allcleara1 = -0.1950;a2 = 0.95;sigma = 0.0965;e2=zeros(151,2048);c11=zeros(151,2048);c21=zeros(151,2048);x=zeros(151,2048);for k=1:150w=sigma*randn(1,2048);c1=zeros(1,2049);c2=zeros(1,2049);k1=0.04;for n=3:2048x(n)=w(n)-a1*x(n-1)-a2*x(n-2);y(n)=c1(n-1)*x(n-1)+c2(n-2)*x(n-2);e(n)=x(n)-y(n);e1(n)=e(n)*e(n);c1(n+1)=c1(n)+2*k1*e(n)*x(n-1);c2(n+1)=c2(n)+2*k1*e(n)*x(n-2);end;e2(k+1,:)=e1(1:2048)+e2(k,:);c11(k+1,:)=c1(1:2048)+c11(k,:);c21(k+1,:)=c2(1:2048)+c21(k,:);end;e31=e2(151,:)/150;c13=c11(151,:)/150;c22=c21(151,:)/150;E=mean(e31);co1=max(c13);co2=min(c22);figure(1)plot(c1,c2,'k');hold on；axis([-1.4 1.6 -2.5 0.5]);xlabel('c1');ylabel('c2');title('Averaged trajectory');figure(2)plot(e1,'r')；hold on；plot(1:2048,E,'g--');hold on;plot(e31,'k');hold off ；xlabel('Length of signal (n)');ylabel('P(n)');title('MSE P(n)learning curve');figure(3)plot(c1,'m');hold on; plot(1:2048,co1,'g--');hold on；plot(c2,'y');hold on;plot(c13,'k');hold on;plot(c22,'r');hold on; plot(1:2048,co2,'g--');hold offxlabel('Length of signal (n)');ylabel('Coefficients');title('c(n)learning curve'); text(200,0.05,'c1(n)'); text(750,0.25,'co1');text(400,-0.6,'c1(n)'); text(750,-0.97,'co1');.。