现代信号处理大作业题目+答案

研究生“现代信号处理”课程大型作业

（以下四个题目任选三题做）

1. 请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为，并画出学习曲线。其中，非线性函数X [0 0;0 1;1 0;1 1]T，要求可以判别输出为0或1）采用S型Logistic函数。

2. 试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：Fp＝1.7KHz, Fr＝2.3KHz, Fs＝8KHz, Armin≥70dB。

3. 根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据(pp.352-353)，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线： 1） Levinson算法 2） Burg算法 3） ARMA模型法 4） MUSIC算法

4. 图1为均衡带限信号所引起失真的横向或格型自适应均衡器(其中横向FIR系统长M=11), 系统输入是取值为±1的随机序列x(n)，其均值为零；参考信号d(n) x(n 7)；信道具有脉冲响应：

2 (n 2) 1

)] n 1,2,3 [1 cos(

h(n) 2W

0 其它

式中W用来控制信道的幅度失真(W = 2～4, 如取W = 2.9,3.1,3.3,3.5等)，且信道受到均

2

值为零、方差 v 0.001(相当于信噪比为30dB)的高斯白噪声v(n)的干扰。试比较基

于下列几种算法的自适应均衡器在不同信道失真、不同噪声干扰下的收敛情况(对应于每一种情况,在同一坐标下画出其学习曲线): 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）算法实现异或问题（输入为X [0 0;0 1;1 0;1 1]T，要求可以判别输出为0或1），并画出学习曲线。其中，非线性函 数采用S型Logistic函数。

1、原理：

反向传播(BP)算法：

(1)、多层感知器的中间隐层不直接与外界连接，其误差无法估计。

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

正向过程：从输入层经隐层逐层正向计算各单元的输出

反向过程：由输出层误差逐层反向计算隐层各单元的误差，并用此误

差修正前层的权值。 2、流程图：

j

3、程序：

%使用了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 < 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) <= threshold| circle_time>3001 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;

%then

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:r

result(i,i) = (1-a(i))*a(i); end;

4、实验结果：

计算误差 = 0.0049993 迭代次数 = 2706 0 0 = 0.023182 0 1 = 0.963110 1 0 = 0.965390 1 1 = 0.043374 5、学习曲线图：

图1.MLP

二、试用用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。滤波器设计参数为：Fp＝1.7KHz, Fr＝2.3KHz, Fs＝8KHz, Armin≥70dB。

1、设计步骤：

(1)对Fp、Fr进行预畸

Fp

);Fs

Fr

'r tg();

Fs

'p tg(

(2)计算 'c 'p* 'r，判断 'c是否等于1，即该互补滤波器是否为互补镜像滤波器

(3)计算相关系数

1－12

k1 [(10k

'p 'r

;

0.1Apmin

1)(10

0.1Armin

1)];

k' k2;

11 k'

q0 ;

21 k'

5913

q q0 2q0 15q0 150q0;N lg(k12/16)/lgq; i;( N为奇数) u 1

i ;(N为偶数) 2

2q i'

1

2m 0

( 1)

m 1

m

qm(m 1)sin(

2

1 2 ( 1)mqm

(2m 1)

u)N

;

2m cos(u)

N

N N

N2 ; N1 N2;

4 2 vi (1 k i'2)(1 i'2/k); i

2v2i 1

;i 1, N1 '2

1 2i 1

i

2v2i

;i 1, N2

1 '22i

2 i2 i

; Bi ; 2 i2 i

(4)互补镜像滤波器的数字实现

Ai

Ai Z 2Bi Z 2 1

H1(Z) i 1, N1 H2(Z) Z ;i 1, N2 2 2

1 AZ1 BZiiii

1

HL(Z) [H1(Z) H2(Z)];

2

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: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 a

b=0;

for m=1:5

b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N); end b

W(jj)=2*q^0.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+Wc^2)/(1+alpha(i)*Wc+Wc^2); end

for i=1:N2

b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2); end

w=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 filterfour

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: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^(m^2)*cos(2*m*pi*jj/N); end

W(jj)=2*q^0.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+Wc^2)/(1+alpha(i)*Wc+Wc^2); end

for i=1:N2

b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2); end

w=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:N2

H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));

end

H22=1/(z^2); for i=1:N2

H22=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)、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：

(m) 1Rxx

N

N 1 m

x(n)x(m n),

*n 0k

m N 1

(2)、根据估计所得到的自相关函数，用下面的迭代公式估算AR模型参数：

R(0) ak,iR*(i)

2k

i 1

ak,0 0

Dk ak,iR(k 1 i),

i 0

k

k 1

Dk

k2

2

k2 1 (1 k 1) k2

*

ak 1,i ak,i k 1ak,k 1 i,

i 1,2, ,k

ak 1,k 1 k 1

(3)、对于AR（p）模型，按以上述递推公式迭代计算直到k 1 p时为止。将算出的模型参数代入下式即可得到功率谱估计值：

(e) Sxx

j

2p

ap,ie jwi

i 1p

2

程序：

function [sigma2,a]=levinson(signal_source,p)

%阶数由p确定

N=length(signal_source); %确定自相关函数 for m=0:N-1

R(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)]; end

a=[1 a1];

%计算功率谱估计值 sigma2=real(sigma2);

p=15;%p分别为15、30、45、60

[sigma2,a]=Levinson(signal_in_complex,p); %计算功率谱

f1=linspace(-0.5,0.5,512); %从-0.5到0.5生成512个等间隔数据

for k=1:512

S1(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)、设输入数据序列为x(n)，0 n N 1，对前后向预测误差之和求偏导，得反射系数

*

2 ek(n)e 1k 1(n 1)n kN 1

k

(e

n k

N 1

k 1

(n) ek 1(n 1))

22

前后向预测误差递推公式：

ek(n) 1 e(n) * k k

k ek 1(n) e (n 1) 1 k 1

ak,i ak 1,i kak 1,k i,i 1,2,3,...,k,ak 1,0 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模型的参数就全部估计出来了，根据以下公式即可算出功率谱：

(ej ) 2Sxxp

bie jwi

i 1p

2

q

2

aie jwi

i 1

程序：

function [a,b,sigma2]=arma(signal_source,p,q) N=length(signal_source); M=N-1;

r=xcorr(signal_source,'unbiased'); for k=1:p

R(:,k)=(r([N+q-k+1:N+M-k])).'; end

a1=(-pinv(R)*(r([N+q+1:N+M]).')).'; a=[1 a1];

Y=filter([1 a1],[1],signal_source); pp=4*q;

[sigma,a1]=burg(Y,pp); sigma2=sigma(2:end);

[sigma,a2]=burg(a1(2:end),q); b=a2;

实验结果：

图6.ARMA算法

4、 MUSIC算法

分析：

(1)构造自相关阵

R( 1) R( (N 1)) R(0)

R(0) R( (N 2)) R(1)

R

R(N 1)R(N 2) R(0)

自相关函数可用有偏估计代替。

i，i=1,2,…,N。 (2)计算R的特征向量v

(3)估计R的维数M，方法有AIC和MDL法。

(4)根据剩余特征向量计算MUSIC谱。