现代信号处理大作业

邱辉1010010218现代信号处理大型作业一．试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：F＝1.7KHz,F＝2.3KHz,Fprs＝8KHz,A≥70dB。rmin与通常的滤波器相比，互补滤波器具有优良的结构特性和结构特性，具有较低的噪声能量和系数敏感性，其定义如下：一组滤波器H(ZH(ZH(Z)如果满足下式：12nnH(e)0<w<2jwKk1则称这组滤波器为幅度互补滤波器;如果满足下式：2nH(e)1,0<w<2jwkk1则称这组滤波器为功率互补滤波器，同时互补滤波器还应该满足：nH(z)(z)其中A(z)为全通函数，适当的选择全通函数，可以kk1使两带函数具有所需要的低通和高通特性。(1)对Fp、Fr进行预畸tg('ptg('r(2)计算*，判断是否等于1，即该互补滤波器是否为互补'c'p'r'c镜像滤波器(3)计算相关系数第页1邱辉1010010218110.1Apmin0.1Armin－];12'p'k;rk'1k2;11k''q0;21kqq2q15q150q13;509000Nlg(1/16)/lgq;2;(N)iu1i;(N)2(2m12q(qsin(u)mm(m2N2m;'m0i12(q2cos(u)mmNm1N4N;2N2NN;12vk'2'2/k);iiiv;iN2i11'22i1i12v;iN2i1i'222i(4)互补镜像滤波器的数字实现2A;;ii2i22BiiiAZ2H(Z)i;N1i1AZ12ii第页2邱辉1010010218BZ2H(Z)Z;iN1i1BZ222ii1H(Z)[H(Z)H(Z)];2L121．二带滤波器组（1）源程序：clear;clf;Fp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;forjj=1:Ma=0;form=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%Nisodd,u=jendab=0;form=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endbW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));endfori=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfori=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfori=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfori=1:N2第页3邱辉1010010218b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LP=zeros(size(w));HP=zeros(size(w));forn=1:length(w)z=exp(j*w(n)*2*pi);H1=1;fori=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2));endH2=1/z;fori=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endLP(n)=abs((H1+H2)/2);HP(n)=abs((H1-H2)/2);endplot(w,LP,'b',w,HP,'r');holdon;xlabel('digitalfrequency');ylabel('amptitude');（2）运行结果：见图1图1二带数字滤波器组2．四带滤波器组（1）源程序：第页4邱辉1010010218clear;clf;Fp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;forjj=1:Ma=0;form=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%Nisodd,u=jendb=0;form=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));Endfori=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfori=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfori=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfori=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LLP=zeros(size(w));LHP=zeros(size(w));HLP=zeros(size(w));HHP=zeros(size(w));forn=1:length(w)z=exp(j*w(n)*2*pi);H1=1;fori=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2));endH21=1;fori=1:N1H21=H21*(a(i)+z^(-4))/(1+a(i)*z^(-4));第页5邱辉1010010218endH2=1/z;fori=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endH22=1/(z^2);fori=1:N2H22=H22*(b(i)+z^(-4))/(1+b(i)*z^(-4));endLP=((H1+H2)/2);HP=((H1-H2)/2);LLP(n)=abs((H21+H22)/2*LP);LHP(n)=abs((H21-H22)/2*LP);HHP(n)=abs((H21+H22)/2*HP);HLP(n)=abs((H21-H22)/2*HP);endplot(w,LLP,'b',w,LHP,'r',w,HLP,'k',w,HHP,'m')holdonxlabel('digitalfrequency');ylabel('amptitude');（2）运行结果：见图2图2四带数字滤波器组二、根据《现代数字信号处理》第四章提供的数据，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1）Levison算法2）Burg算法第页6邱辉10100102183）ARMA模型法4）MUSIC算法1Levinson算法Levinson算法用于求解Yule-Walker方程，是一种按阶次进行递推的算法，即首先以AR（0）和AR（1）模型参数作为初始条件，计算AR（2）模型参数；然后根据这些参数计算AR（3）参数，等等，一直到计算出AR（p）模型参数为止，需要的运算量数量级为2，其中p为AR模型的阶数。p它利用了方程系数矩阵的对称性和Toeplitz性质，是一种高效的算法。Levinson算法的优点是计算简单，步步监视均方误差，其缺点是需要由观测数据计算自相关值，当观测数据短时，误差较大；当观测数据长时，计算量大；并且会产生谱峰漂移和谱线分裂。算法步骤如下：（1）由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：1ˆNmR(m)xxxnxmn()(),1mN*Nn0（2）根据估计所得到的自相关函数，用下面的迭代公式估算数：AR模型参kR(0)aR(i)2*kk,ii1kDaRk(1),ia0k,0kk,ii0Dk2kk1aa22k12k1kaa*,i,kkik,ik1k,kik1,k1k1（3）对于AR（p）模型，按以上述递推公式迭代计算直到k1p时为止。将算出的模型参数代入下式即可得到功率谱估计值：2ˆS(e)jpxx2p1aejwip,ii1第页7邱辉1010010218Levinson算法的MATLAB源程序如下：X为输入序列，p为ARLevinson（X，pfunctionS=Levinson(X,p)N=length(X);form=0:N-1R(m+1)=sum(conj(X([1:N-m])).*X([m+1:N]))/N;enda=-R(2)/R(1);sigma2=(1-abs(a)^2)*R(1);gamma=-a;fork=2:psigma2(k)=R(1)+sum(a.*conj(R([2:k])));D=R(k+1)+sum(a.*R([k:-1:2]));gamma(k)=D/sigma2(k);sigma2(k)=(1-abs(gamma(k))^2)*sigma2(k);a=[a-gamma(k)*conj(fliplr(a)),-gamma(k)];endsigma2=real(sigma2);f=linspace(-0.5,0.5,512);fork=1:512S(k)=10*log10(sigma2(end)/(abs(1+sum(a.*exp(-j*2*pi*f(k)*[1:p])))^2));endplot(f,S);title(['Levinson:',int2str(p),'阶']);xlabel('归一化频率'),ylabel('相对谱/dB');分别对应于10阶、25阶、40阶、55阶的Levinson算法的源程序和执行结果如下所示：clfp=[10254055];fork=1:4subplot(2,2,k);Levinson(x,p(k));end第页8邱辉1010010218图3lenxison算法2Burg算法Burg算法一方面希望利用已知数据段两端以外的未知数据（它对这些未知数据不作主观臆测）小相位的。Burg算法与自相关法和协方差法不同，它不直接计算AR参数，而是先估计反射系数，然后利用Levinson递推算法由反射系数求得AR参数，估计反射系数时所使用的准则是前向和后向预测误差功率估计的平均值最小准则，在这里，预测误差功率估计仍然用时间平均来代替集合平均。Burg法估计反射系数的准则表示为：N1[e(n)e(n)]min22ppnp其优点是实现简单，在一定程度上可以克服Levinson算法的谱峰漂移和谱线分裂的缺点。算法步骤如下：（1）设输入数据序列为(01，对前后向预测误差之和求偏xnnN导，得反射系数N12e(n)e(nk1*k1nkk1N(e(n)2e(n2)k1k1nk前后向预测误差递推公式：第页9邱辉1010010218()1()enenkkk11e(n)(*enkkk1aaa,i1,2,3,...,k,a1k,ik1,ikk1,kik1,0（2k=pAR模型参数计算功率谱，计算功率谱的公式同上面算法。Burg算法的MATLAB源程序如下：其中，参数X为输入序列，p为AR模型的阶数，函数调用形式为：Burg（X，pfunctionS=Burg1(X,p)N=length(X);ef=X;eb=X;sigma2=sum(X*X')/N;a=[];fork=1:pefp=ef(k+1:end);ebp=eb(k:end-1);gamma(k)=2*efp*ebp'/(efp*efp'+ebp*ebp');sigma2(k+1)=(1-abs(gamma(k))^2)*sigma2(k);ef(k+1:end)=efp-gamma(k)*ebp;eb(k+1:end)=ebp-gamma(k)'*efp;a=[a-gamma(k)*conj(fliplr(a)),-gamma(k)];endf=linspace(-0.5,0.5,512);fork=1:512S(k)=10*log10(sigma2(end)/(abs(1+sum(a.*exp(-j*2*pi*f(k)*[1:p])))^2));endplot(f,S);title(['Burg:',int2str(p),'阶']);xlabel('归一化频率'),ylabel('相对谱/dB');分别对应于10阶、25阶、40阶、55阶的Burg算法的源程序和执行结果如下所示：clfp=[10254055];fork=1:4subplot(2,2,k);Burg1(x,p(k));第页10邱辉1010010218end图4Burg算法3ARMA模型法当采用AR模型谱估计方法，特别是采用Burg法时，能得到可靠的高分辨率估计。但当噪声污染了数据时，只有采用ARMA模型才能获得良好的谱估计。采用ARMA模型，以较少的模型参数就能改善AR谱估计的性能。ARMA模型法的第一步是模型定阶。定阶方法基本上分为两类，一类是然后再逐一比较，计算量很大。另一类是线性代数法，这类方法不用预先估计模型参数，计算简单。对于本题，最终确定阶数p=17,q=6。ARMA模型法的第二步是计算模型参数。ARMA模型的参数估计有两种方法，一种是基于最大似然估计的最佳估计法，这种方法计算复杂。另一种是次最佳的方法，它先估计模型的AR参数，再估计MA参数，而如果仅仅是计算功率谱，则计算MA参数的步骤可以省略。本题采用修正的Yule-Walker法估算AR参数，它需要求解如下的（M-q）×p阶超定方程组R(q)R(qR(qpR(qa1a(RqpRqp(2)Rq()()2Rqpa(MR(M2)R(Mp)R(M)Rp第页11邱辉1010010218其中M是自相关函数的最大延迟。若R(m)用自相关函数的无偏估计1ˆNm()x()(1nxmnmNR'*Nmn0代替，则有其中aRrTaaaa12p为方程组的最小二乘解。R可通过SVDMATLAB函数pinv求解。求出AR模型参数后，将其构成A(z)并与原来的ARMA模型串联，这样求解B(z)的问题就可以转化为求解MA模型参数的问题。算法步骤如下：（1）用x(n)通过A(z)，得到y(n)。（2）用一无穷阶的AR模型近似MA模型。用Burg算法可得到此近似AR模型的参数以及激励白噪声的功率。一般此AR模型的阶数应大于MA模型阶数的两倍，以获得较好的近似效果。（3）可以证明，将上一步求出的近似AR模型参数视为时间序列，则MA模型就可视为一线性预测滤波器，按最小均方误差准则就可以求出型参数，方法同样可用Burg算法。MA模这样，ARMA模型的参数就全部估计出来了，根据以下公式即可算出功率谱：22q1bejwiˆiS(e)i1j2xxpp1aejwiii1ARMA算法的MATLAB源程序如下：其中，参数X为输入序列，p为AR模型的阶数，q为MA模型的阶数，函数调用形式为：ARMA（X，p，qfunctionS=ARMA(X,p,q)N=length(X);M=N-1;r=xcorr(X,'unbiased');fork=1:pR(:,k)=(r([N+q-k+1:N+M-k])).';end第页12邱辉1010010218a=(-pinv(R)*(r([N+q+1:N+M]).')).';Y=filter([1,a],[1],X);pp=5*q;[S,A,K]=Burg(Y,pp);P=K(end);[S,A,K]=Burg(A,q);b=A;f=linspace(-0.5,0.5,512);fork=1:512S(k)=10*log10(P*(abs(1+b*(exp(-j*2*pi*f(k)*[1:q]).'))/abs(1+a*(exp(-j*2*pi*f(k)*[1:p]).')))^2);endplot(f,S);title(['ARMA:(',int2str(p),',',int2str(q),')阶']);xlabel('归一化频率'),ylabel('相对谱/dB');分别对应于ARMA（30，10ARMA（30，1ARMA（40，10ARMA（50，10）的ARMA算法的源程序和执行结果如下所示：clfsubplot(2,2,1);ARMA(x,30,10);subplot(2,2,2);ARMA(x,30,1);subplot(2,2,3);ARMA(x,40,10);subplot(2,2,4);ARMA(x,50,10);第页13邱辉1010010218图5ARMA算法4MUSIC算法MUSIC算法是根据自相关矩阵的特征分解，把噪声和信号区分开来。这种方法抗干扰能力强、稳定性好、分辨率高。MUSIC算法的步骤：（1）构造自相关阵R(0)R(R((N)R(0)R((N2))RRR(NR(N2)R(0)自相关函数可用有偏估计代替。（2）计算R的特征向量ˆ，=1,2,…,。iNi（3）估计R的维数M，方法有AIC和MDL法。（4）根据剩余特征向量计算MUSIC谱。1ˆP(f)NHˆev2iiM1其中，ee,e,,e)。TjN(jj第页14邱辉1010010218对于向量相乘，可用FFT来减少运算量。MUSIC算法的MATLAB源程序如下：其中，参数X为输入序列，p为AR模型的阶数，函数调用形式为：MUSIC（X，pfunctionS=MUSIC(X,p)N=length(X);r=xcorr(X,'biased');clearRfork=1:NR(:,k)=(r([N-(k-1):2*N-k])).';end[V,D]=eig(R);f=linspace(-0.5,0.5,512);S0=fft(V(:,p+1:end),512);fork=1:512S(k)=10*log10(1/(S0(k,:)*S0(k,:)'));endS=[S(257:512)S(1:256)];plot(f,S);title(['MUSIC:',int2str(p),'维']);xlabel('归一化频率'),ylabel('相对谱/dB');分别对应于10维、25维、40维、55维的MUSIC算法的源程序和执行结果如下所示：clfp=[10254055];fork=1:4subplot(2,2,k);MUSIC(x,p(k));end第页15邱辉1010010218图6MUSIC算法第页16邱辉10100102181为均衡带限信号所引起失真的横向或格型自适应均衡器(其中横向FIR系统长M=11),系统输入是取值为±1的随机序列()，其均值为零；xn参考信号()(7)；信道具有脉冲响应：dnxncos(2(n2)/)]/，nh(n)0其它式中w用来控制信道的幅度失真(w=2～4,如取w=2.9,3.1,3.3,3.5等)，且信道受到均值为零、方差(相当于信噪比为30dB)的高斯白噪2v声v(n)的干扰。试比较基于下列五种算法的自适应均衡器在不同信道失(对应于每一种情况,在同一坐标下画出其学习曲线):1）横向/格-梯型结构LMS算法;2）横向/格-梯型结构RLS算法;并分析其结果。延时++信道+图1所谓自适应滤波器，就是利用前一时刻已经获得的滤波器参数等结果，自动的调整现在时刻的滤波器参数，以适应信号和噪声未知的或随时间变化的统计特性，从而实现最优滤波。性能优于以前所研究的维纳滤波或者卡尔曼滤波器。第页17邱辉1010010218㈠基于横向结构的LMS算法1.设计思路：用横向自适应滤波器的输出与期望响应信号的误差调节滤波器抽头系数，使得均方误差最小。LMS算法实际采用最速下降法得出抽头系数的迭代公式。2.设计参数：ekdkxT(k)WWk1Wk(k)(k)(k)(k)3.源程序：clear;clf;w=[2.9,3.1,3.3,3.5];M=11;sgm2=0.001;%噪声方差T=7;N=400;%训练次数iteration=500;L=iteration+T+M-1;%单个输入信号长度u=0.025;%迭代步长value=zeros(length(w),L-M+1-T);%-------forww=1:length(w)h=ones(1,3);h(1)=0.5*(1+cos(2*pi/w(ww)));h(3)=h(1);e2=zeros(N,L-M+1-T);%-------forn=1:Nrand('seed',n*N);X=sign(2*rand(1,L)-1);X=randsrc(1,L,[+1,-1],1)D=zeros(size(X));;%产生一个随机信号序列formm=T+1:LD(mm)=X(mm-T);end%对应的参考信号U=conv(X,h);randn('seed',n*N);V=randn(size(U)).*sqrt(sgm2);%产生高斯噪声R=U+V;%自适应滤波器输入信号%-------W=zeros(M,1);%滤波器参数的初始值为0form=1:iterationr=R(T+m:T+m+M-1)';y=r'*W;e=D(T+m+M-1)-y;%误差信号e2(n,m)=e.^2;W=W+2*u*e*r;%滤波器参数迭代第页18邱辉1010010218endendvalue(ww,:)=mean(e2);%均方误差值endsemilogy(value(1,:),'k-');holdon;semilogy(value(2,:),'b-.');semilogy(value(3,:),'m:');semilogy(value(4,:),'r--');grid;xlabel('thenumberofiteration');ylabel('meansquareerror');legend('W=2.9','W=3.1','W=3.3','W=3.5');4.运行结果：见图7、8不难看出，取得合适的收敛因子μ，对横向LMS算法的好坏有着重要的影响。μ取得较大，则收敛较快，但收敛时的稳态误差也大；μ取得较小，则收敛较慢，但稳态误差也小。若μ过大，则会不收敛；μ过小，会因收敛速度太慢而在指定的时间完不成收敛，它们都是不可取的。图7横向LMS（u=0.025）第页19邱辉1010010218图8横向LMS（u=0.035）㈡基于横向结构的RLS算法1使加权遗忘因子后的误差平方和取得最小值。．设计参数：⑴初始值P(0)1W(0)0e(0)d(0)⑵迭代公式：ˆendnx(n)W(nTPn1xnKnxnn1xPnTPnPn1KnxnPn111TˆˆWnWn1K(n)en00100x(()[(0,,0],Z＝0nnxnT01003．源程序：clear;w=[2.9,3.1,3.3,3.5];M=11;第页20邱辉1010010218sgm2=0.001;T=7;%参考信号延迟N=400;%训练次数iteration=500;L=iteration+T+M-1;%输入信号长度;100是迭代次数lamda=0.899;%迭代步长value=zeros(length(w),L-M+1-T);%-------forww=1:length(w)h=ones(1,3);h(1)=0.5*(1+cos(2*pi/w(ww)));h(3)=h(1);e2=zeros(N,L-M+1-T);%-------forn=1:Nrand('seed',n*N);X=sign(2*rand(1,L)-1);D=zeros(size(X));formm=T+1:LD(mm)=X(mm-T);endU=conv(X,h);randn('seed',n*N);V=randn(size(U)).*sqrt(sgm2);R=U+V;%"R"是自适应滤波器输入%-------W=zeros(M,1);P=eye(M)/0.004;form=1:iterationr=R(T+m+M-1:-1:T+m)';e=D(T+m+M-1)-r'*W;e2(n,m)=e.^2;K=P*r/(lamda+r'*P*r);P=(eye(M)-K*r')*P/lamda;W=W+K.*e;formm=2:M-1r(mm)=r(mm+1);endr(1)=R(T+m+M);endendvalue(ww,:)=mean(e2);end%-------clf;semilogy(value(1,:),'k-');holdon;semilogy(value(2,:),'b-.');semilogy(value(3,:),'m:');semilogy(value(4,:),'r--');grid;v{1}=['\mu=',num2str(mu)];title(v);xlabel('thenumberofiteration');ylabel('meansquareerror');第页21邱辉1010010218legend('W=2.9','W=3.1','W=3.3','W=3.5');4．运行结果：见图9图9横向RLS可见横向RLS收敛速度要比横向LMS块很多。(三)基于格-梯型结构的LMS算法1．设计思路：按照最小均方准则，设计出阶数和时间分别递推的自适应滤波器，可以进行梯型滤波，又具有格型预测前后节独立、收敛快速的优点，但计算量大。2．设计参数：⑴初始值f(n)g(n)x(n)00f(g(g(2)0mmme(n)d(n),e(0)m10mv(0m(0)(0kmm⑵迭代公式第页22邱辉1010010218v(n)wv(n|f(n)||g(n|22mmm1m1f(ngm1(ng(n2)f(n**m1v(nk(n)k(nmmmmmf(n)f(n)k(n)g(nm1mm1mg(n)g(nkfn1*mm1mm1e(n)e(n)g(n)m(n)m1mmd(n)β(nG(n)TMM2e(n)g(n)*β(nβ(n)mmv(n)mmm3.源程序：clear;W=[2.9,3.1,3.3,3.5];w=0.9999;T=7;sgm2=0.001;L=500;%数据长度M=11;N=400;%次数e2=zeros(length(W),L);forww=1:length(W)n=[1,2,3];h(n)=1/2*(1+cos(2*pi*(n-2)/W(ww)));fornn=1:Nkm=ones(L+2,M+1);vm=ones(L+2,M+1);beta=ones(L+2,M+1);fm=ones(L+2,M+1);gm=ones(L+2,M+1);e=ones(L+2,M+2);rand('seed',N*nn);X=sign(2*rand(1,L+T)-1);U=conv(X,h);randn('seed',N*nn);V=randn(size(U))*sqrt(sgm2);R=U+V;form=1:M+1km(2,m)=0;%n=2vm(2,m)=0.8;beta(3,m)=0;fm(2,m)=0;gm(2,m)=0;gm(1,m)=0;end第页23邱辉1010010218forn=3:L+2%m=1e(n,2)=X(n-2);fm(n,1)=R(n-2+T);gm(n,1)=R(n-2+T);endforn=3:L+2form=2:M+1vm(n,m)=w*vm(n-1,m)+(abs(fm(n,m-1)))^2+(abs(gm(n-1,m-1)))^2;km(n,m)=km(n-1,m)+(fm(n-1,m-1)*conj(gm(n-1,m))+conj(gm(n-2,m-1))*fm(n-1,m))/vm(n-1,m);fm(n,m)=fm(n,m-1)-km(n,m)*gm(n-1,m-1);gm(n,m)=gm(n-1,m-1)-conj(km(n,m))*fm(n,m-1);e(n,m+1)=e(n,m)-gm(n,m)*beta(n,m);beta(n+1,m)=beta(n,m)+2*e(n,m+1)*conj(gm(n,m))/vm(n,m);ende2(ww,n-2)=e2(ww,n-2)+e(n,M+2)*e(n,M+2);endendende2=e2/N;clf;semilogy(1:L,e2(1,:),'k-');holdon;semilogy(1:L,e2(2,:),'b-.');semilogy(1:L,e2(3,:),'m:');semilogy(1:L,e2(4,:),'r--');gridon;xlabel('data');ylabel('meansquareerror');legend('W=2.9','W=3.1','W=3.3','W=3.5');4.运行结果：见图10第页24邱辉1010010218图10格-梯形LMS格-梯形LMS的收敛速度很慢，而且其稳定性也不是很好。（四）基于格-梯型结构的RLS算法设计思路：结合格型滤波器阶数可调与横型滤波器时间可调的特性，设计出阶数和时间分别递推的自适应滤波器，使得平方误差最小。1．设计参数：⑴初始值(ne(n)dn(，f(n)g(n)x(n)0000E(n)E(n)(nx(n)2f0b0f0(，(0kmmE(E(0)，(0bfmmm⑵格型预测k(nwk(n2)(n2)f(ng(n2)*m1m1mmk(nH(nfm1E(n2)m1bmk*(nH(nbm1m1E(nfmf(n)f(n)H(ng(nfm1mm1mg(n)g(nH(nf(n)bm1mm1m第页25邱辉1010010218k(n2E(nE(n1ffm(n2)m1mEbmk(n2E(nE(n2)1bm1bm(nmEfm(ngn(22((nnmm(nm1mEbm⑶梯型滤波(nw(2)(nng*nen((mmmm(n(nmE(nmbme(n)e(n)(ng(n)m1mmm2．源程序：clear;W=[2.9,3.1,3.3,3.5];w=0.99;T=7;L=500;%数据长度M=11;N=400;%次数sgm2=0.001;e2=zeros(length(W),L);forww=1:length(W)h=ones(1,3);h(1)=0.5*(1+cos(2*pi/W(ww)));h(3)=h(1);alpha=zeros(L+2,M);k=zeros(L+2,M);kf=zeros(L+2,M);kb=zeros(L+2,M);fm=zeros(L+2,M);gm=zeros(L+2,M);Ef=zeros(L+2,M);Eb=zeros(L+2,M);delta=zeros(L+2,M);kasi=zeros(L+2,M);e=zeros(L+2,M+1);fornn=1:Nrand('seed',N*nn);X=sign(2*rand(1,L+T)-1);D=zeros(size(X));U=conv(X,h);randn('seed',N*nn);V=randn(size(U)).*sqrt(sgm2);R=U+V;%"R"是自适应滤波器输入form=