研究生“现代信号处理”课程大型作业（以下四个题目任选三题做）1. 请用多层感知器（MLP ）神经网络误差反向传播（BP ）算法实现异或问题（输入为[00;01;10;11]X T =，要求可以判别输出为0或1），并画出学习曲线。

其中，非线性函数采用S 型Logistic 函数。

2. 试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。

滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。

3. 根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1） Levinson 算法2） Burg 算法3） ARMA 模型法4） MUSIC 算法4. 图1为均衡带限信号所引起失真的横向或格型自适应均衡器(其中横向FIR 系统长M =11), 系统输入是取值为±1的随机序列)(n x ，其均值为零；参考信号)7()(-=n x n d ；信道具有脉冲响应：12(2)[1cos()]1,2,3()20 n n h n W π-⎧+=⎪=⎨⎪⎩其它式中W 用来控制信道的幅度失真(W = 2～4, 如取W = 2.9,3.1,3.3,3.5等)，且信道受到均值为零、方差001.02=v σ(相当于信噪比为30dB)的高斯白噪声)(n v 的干扰。

试比较基于下列几种算法的自适应均衡器在不同信道失真、不同噪声干扰下的收敛情况(对应于每一种情况,在同一坐标下画出其学习曲线):1） 横向/格-梯型结构LMS 算法2） 横向/格-梯型结构RLS 算法并分析其结果。

图1 横向或格-梯型自适应均衡器参考文献[1] 姚天任, 孙洪. 现代数字信号处理[M]. 武汉: 华中理工大学出版社, 2001[2] 杨绿溪. 现代数字信号处理[M]. 北京: 科学出版社, 2007[3] S. K. Mitra. 孙洪等译. 数字信号处理——基于计算机的方法(第三版)[M]. 北京: 电子工业出版社, 2006[4] S.Haykin, 郑宝玉等译. 自适应滤波器原理(第四版)[M].北京: 电子工业出版社, 2003[5] J. G. Proakis, C. M. Rader, F. Y. Ling, etc. Algorithms for Statistical Signal Processing [M].Beijing: Tsinghua University Press, 2003一、请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为[00;01;10;11]，要求可以判别输出为0或1），并画出学习曲线。

其X T中，非线性函数采用S型Logistic函数。

1、原理：反向传播(BP)算法：(1)、多层感知器的中间隐层不直接与外界连接，其误差无法估计。

(2)、反向传播算法：从后向前（反向）逐层“传播”输出层的误差，以间接算出隐层误差。

分两个阶段：正向过程：从输入层经隐层逐层正向计算各单元的输出反向过程：由输出层误差逐层反向计算隐层各单元的误差，并用此误差修正前层的权值。

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;% 收敛条件∑e^2 < thresholdwd1=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 1temp=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 forif accumulate_error(circle_time) <= threshold| circle_time>3001 %then break;end;%end of ifend;%end of whileplot(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二、试用用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。

滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。

1、设计步骤：(1)对Fp 、Fr 进行预畸);();(''Fs Fr tg Fs Fp tg r p ∏=Ω∏=Ω (2)计算'''*r p c ΩΩ=Ω，判断'c Ω是否等于1，即该互补滤波器是否为互补镜像滤波器(3)计算相关系数⎪⎩⎪⎨⎧-==+++=+-=-=ΩΩ=--=偶数)N 为(;21 奇数)N 为 (;;lg /)16/1lg(;150152;1121;1;;])110)(110[(1213090500''02'''211－min 1.0min 1.0i i u q k N q q q q q k k q k k k k r pAr Ap ;)2cos()1(21))12(sin()1(210)1(21'2∑∑∞=∞=+-++-=Ωm m m m m m m i u Nm q u N m q qππ ;42⎥⎦⎤⎢⎣⎡=N N ;221N N N -⎥⎦⎤⎢⎣⎡= ;)/1)(1(2'2'k k v i i i Ω-Ω-= 12'1212,1;12N i v i i i =Ω+=--α 22'22,1;12N i v i i i =Ω+=β (4)互补镜像滤波器的数字实现;22i i i A αα+-= ;22ii i B ββ+-= 1221,1;1)(N i Z A Z A Z H i i i =++=∏-- 22212,1;1)(N i Z B Z B Z Z H i ii =++=∏--- )];()([21)(21Z H Z H Z H L += 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-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;for jj=1:Ma=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%N is odd, u=j endab=0;for m=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));endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfor i=1:N2b(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));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH2=1/z;for i=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,'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-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;for jj=1:Ma=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N); % N is odd, u=j endb=0;for m=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));endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfor i=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));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH21=1;for i=1:N1H21=H21*(a(i)+z^(-4))/(1+a(i)*z^(-4)) ;endH2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endH22=1/(z^2);for i=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,'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)、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：∑--=-≤+=m N n xxN m n m x n x Nm R 10*1),()(1)(ˆ(2)、根据估计所得到的自相关函数，用下面的迭代公式估算AR 模型参数：)()0(*1,2i R a R ki i k k∑=+=σ∑==-+=ki k i k k a i k R a D 00,,0),1(21kkk D σγ=+22121)1(k k k σγσ++-=k i a a a i k k k i k i k ,,2,1,*1,1,,1 =-=-+++γ11,1+++-=k k k a γ(3)、对于AR （p ）模型，按以上述递推公式迭代计算直到p k =+1时为止。