南邮现代信号处理最后大作业4道题目(含答案).doc

南邮研究生“现代信号处理”期末课程大作业（四个题目任选三题做）1. 请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为，要求可以判别输出为0或1），并画出学习曲线。其中，非线性函数采用S型Logistic函数。2. 试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：Fp1.7KHz, Fr2.3KHz, Fs8KHz, Armin70dB。3. 根据现代数字信号处理（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1） Levinson算法2） Burg算法3） ARMA模型法4） MUSIC算法4. 图1为均衡带限信号所引起失真的横向或格型自适应均衡器(其中横向FIR系统长M=11), 系统输入是取值为1的随机序列，其均值为零；参考信号；信道具有脉冲响应：式中用来控制信道的幅度失真(W = 24, 如取W = 2.9,3.1,3.3,3.5等)，且信道受到均值为零、方差(相当于信噪比为30dB)的高斯白噪声的干扰。试比较基于下列几种算法的自适应均衡器在不同信道失真、不同噪声干扰下的收敛情况(对应于每一种情况,在同一坐标下画出其学习曲线):1） 横向/格-梯型结构LMS算法2） 横向/格-梯型结构RLS算法并分析其结果。图1 横向或格-梯型自适应均衡器一、请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为，要求可以判别输出为0或1），并画出学习曲线。其中，非线性函 数采用S型Logistic函数。1、原理：反向传播(BP)算法：(1)、多层感知器的中间隐层不直接与外界连接，其误差无法估计。(2)、反向传播算法：从后向前（反向）逐层“传播”输出层的误差，以间接算出隐层误差。分两个阶段：正向过程：从输入层经隐层逐层正向计算各单元的输出反向过程：由输出层误差逐层反向计算隐层各单元的误差，并用此误差修正前层的权值。选择初始值开始前向计算求所有神经元的输出 计算输出层从后向前计算隐层计算保存权值修正量修正权值是否收敛？结束YN 2、流程图：j3、程序：%使用了3层结构，第二层隐藏层4个单元。2，3层都使用Logisitic函数。 %训练xor数据。 function mlp()f= fopen(XOR.txt); A = fscanf(f, %g,3 inf); A = A; p = A(1:2, :);%训练输入数据 t = A(3, :);%desire out train_num , input_scale= size(p) ;%规模 fclose(f); accumulate_error=zeros(1,3001); alpha = 0.5;%学习率 threshold = 0.005;% 收敛条件 e2 threshold wd1=0; wd2=0; bd1=0; bd2=0; circle_time =0; hidden_unitnum = 4; %隐藏层的单元数 w1 = rand(hidden_unitnum,2);%4个神经元，每个神经元接受2个输入 w2 = rand(1,hidden_unitnum);%一个神经元，每个神经元接受4个输入 b1 = rand(hidden_unitnum,1); b2 = rand(1,1);while 1 temp=0; circle_time = circle_time +1; for i=1:train_num %前向传播 a0 = double ( p(i,:) );%第i行数据 n1 = w1*a0+b1; a1 = Logistic(n1);%第一个的输出 n2 = w2*a1+b2; a2 = Logistic(n2);%第二个的输出 a = a2; %后向传播敏感性 e = t(i,:)-a; accumulate_error(circle_time) = temp + abs(e)2; temp=accumulate_error(circle_time); s2 = F(a2)*e; %输出层delta值 s1 = F(a1)*w2*s2;%隐层delta值 %修改权值 wd1 = alpha .* s1*a0; wd2 = alpha .* s2*a1; w1 = w1 + wd1; w2 = w2 + wd2; bd1 = alpha .* s1; bd2 = alpha .* s2; b1 = b1 + bd1; b2 = b2 + bd2; end;%end of for if accumulate_error(circle_time) 3001 %then break; end;%end of if end;%end of while plot(accumulate_error,m); grid; xlabel(学习次数) ylabel(误差)disp(计算误差 = ,num2str(accumulate_error(circle_time) );disp(迭代次数 = ,num2str(circle_time);%测试a0 = double (0 0);n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(0 0 = ,num2str(a);a0 = double (0 1);n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(0 1 = ,num2str(a);a0 = double (1 0);n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(1 0 = ,num2str(a);a0 = double (1 1);n1 = w1*a0+b1;a1 = Logistic(n1);n2 = w2*a1+b2;a2 = Logistic(n2);a = a2;disp(1 1 = ,num2str(a);m=0;%-function a= Logistic(n)a = 1./(1+exp(-n);%-function result= F(a)r,c = size(a);result = zeros(r,r);for i =1:rresult(i,i) = (1-a(i)*a(i);end;4、实验结果：计算误差 = 0.0049993迭代次数 = 27060 0 = 0.0231820 1 = 0.9631101 0 = 0.9653901 1 = 0.0433745、学习曲线图：图1.MLP二、试用用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：Fp1.7KHz, Fr2.3KHz, Fs8KHz, Armin70dB。1、设计步骤：(1)对Fp、Fr进行预畸 (2)计算，判断是否等于1，即该互补滤波器是否为互补镜像滤波器(3)计算相关系数 (4)互补镜像滤波器的数字实现 2、程序：function filter2()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-k2);q0=0.5*(1-k1)/(1+k1);q=q0+2*q05+15*q09+150*q013;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;for jj=1:M a=0; for m=0:5 a=a+(-1)m*q(m*(m+1)*sin(2*m+1)*pi*jj/N);%N is odd, u=j enda b=0; for m=1:5 b=b+(-1)m*q(m2)*cos(2*m*pi*jj/N); endb W(jj)=2*q0.25*a/(1+2*b); V(jj)=sqrt(1-k*W(jj)2)*(1-W(jj)2/k);end for i=1:N1 alpha(i)=2*V(2*i-1)/(1+W(2*i-1)2); end for i=1:N2 beta(i)=2*V(2*i)/(1+W(2*i)2); end for i=1:N1 a(i)=(1-alpha(i)*Wc+Wc2)/(1+alpha(i)*Wc+Wc2); end for i=1:N2 b(i)=(1-beta(i)*Wc+Wc2)/(1+beta(i)*Wc+Wc2); endw=0:0.0001:0.5;LP=zeros(size(w);HP=zeros(size(w); for n=1:length(w) z=exp(j*w(n)*2*pi); H1=1; for i=1:N1 H1=H1*(a(i)+z(-2)/(1+a(i)*z(-2) ; end H2=1/z; for i=1:N2 H2=H2*(b(i)+z(-2)/(1+b(i)*z(-2); end LP(n)=abs(H1+H2)/2); HP(n)=abs(H1-H2)/2);end plot(w,LP,k,w,HP,m);%hold on;xlabel(数字频率);ylabel(幅度);3、实验结果：图2.两带滤波器4、四带滤波器组程序：function filterfourFp=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-k2);q0=0.5*(1-k1)/(1+k1);q=q0+2*q05+15*q09+150*q013;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;for jj=1:M a=0; for m=0:5 a=a+(-1)m*q(m*(m+1)*sin(2*m+1)*pi*jj/N); % N is odd, u=j end b=0; for m=1:5 b=b+(-1)m*q(m2)*cos(2*m*pi*jj/N); end W(jj)=2*q0.25*a/(1+2*b); V(jj)=sqrt(1-k*W(jj)2)*(1-W(jj)2/k);end for i=1:N1 alpha(i)=2*V(2*i-1)/(1+W(2*i-1)2); end for i=1:N2 beta(i)=2*V(2*i)/(1+W(2*i)2); end for i=1:N1 a(i)=(1-alpha(i)*Wc+Wc2)/(1+alpha(i)*Wc+Wc2); end for i=1:N2 b(i)=(1-beta(i)*Wc+Wc2)/(1+beta(i)*Wc+Wc2); endw=0:0.0001:0.5; LLP=zeros(size(w);LHP=zeros(size(w); HLP=zeros(size(w);HHP=zeros(size(w);for n=1:length(w) z=exp(j*w(n)*2*pi); H1=1; for i=1:N1 H1=H1*(a(i)+z(-2)/(1+a(i)*z(-2) ; end H21=1; for i=1:N1 H21=H21*(a(i)+z(-4)/(1+a(i)*z(-4) ; end H2=1/z; for i=1:N2H2=H2*(b(i)+z(-2)/(1+b(i)*z(-2); end H22=1/(z2); for i=1:N2H22=H22*(b(i)+z(-4)/(1+b(i)*z(-4); end LP=(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);end plot(w,LLP,k,w,LHP,m,w,HLP,g,w,HHP,b);xlabel(数字频率);ylabel(幅度);5、实验结果：图3.四带滤波器组三、根据现代数字信号处理（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1） Levinson算法2） Burg算法3） ARMA模型法4） MUSIC算法1、 Levinson算法分析：(1)、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：(2)、根据估计所得到的自相关函数，用下面的迭代公式估算AR模型参数：(3)、对于AR（p）模型，按以上述递推公式迭代计算直到时为止。将算出的模型参数代入下式即可得到功率谱估计值：程序：function sigma2,a=levinson(signal_source,p)%阶数由p确定N=length(signal_source);%确定自相关函数for m=0:N-1R(m+1)=sum(conj(signal_source(1:N-m).*signal_source(m+1:N)/N;end%设置迭代初值a1=-R(2)/R(1);sigma2=(1-abs(a1)2)*R(1);gamma=-a1;%开始迭代for k=2:p sigma2(k)=R(1)+sum(a1.*conj(R(2:k); D=R(k+1)+sum(a1.*R(k:-1:2); gamma(k)=D/sigma2(k); sigma2(k)=(1-abs(gamma(k)2)*sigma2(k); a1=a1-gamma(k)*conj(fliplr(a1),-gamma(k);enda=1 a1;%计算功率谱估计值sigma2=real(sigma2);p=15;%p分别为15、30、45、60sigma2,a=Levinson(signal_in_complex,p);%计算功率谱f1=linspace(-0.5,0.5,512); %从-0.5到0.5生成512个等间隔数据for k=1:512S1(k)=10*log10(sigma2(end)/(abs(1+sum(a(2:end).*exp(-j*2*pi*f1(k)*1:p)2); %公式(2.3.7)并以dB表示end；实验结果：图4.Levinson算法2、 Burg算法分析：(1)、设输入数据序列为，对前后向预测误差之和求偏导，得反射系数前后向预测误差递推公式：(2)、重复以上步骤直至k=p，根据迭代得到的AR模型参数计算功率谱，计算功率谱的公式同上面算法。 程序：function sigma2,a=burg(signal_source,p)N=length(signal_source);ef=signal_source; eb=signal_source;sigma2=sum(signal_source*signal_source)/N;a=;for k=1:p efp=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);end;a=1 a;sigma2=real(sigma2);实验结果：图5.Burg算法3、 ARMA算法分析：(1)、用x(n)通过A(z)，得到y(n)。(2)、用一无穷阶的AR模型近似MA模型。用Burg算法可得到此近似AR模型的参数以及激励白噪声的功率。一般此AR模型的阶数应大于MA模型阶数的两倍，以获得较好的近似效果。(3)、可以证明，将上一步求出的近似AR模型参数视为时间序列，则MA模型就可视为一线性预测滤波器，按最小均方误差准则就可以求出MA模型参数，方法同样可用Burg算法。这样，ARMA模型的参数就全部估计出来了，根据