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

南邮研究生“现代信号处理”期末课程大作业

（四个题目任选三题做）

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 ；信道具有脉冲响应：

1

2(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 横向或格-梯型自适应均衡器

一、请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为[00;01;10;11]

X 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 %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: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

二、试用用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带

滤波器组；并画出其频响。滤波器设计参数为：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 为(;2

1

奇数)N 为 (;;lg /)16/1lg(;150152;

1121;

1;

;

])110

)(110[(1213

090500'

'

02'''2

11－min

1.0min

1.0i i u q k N q q q q q k k q k k k k r

p Ar Ap

;)

2cos()1(21))12(sin(

)

1(21

)1(21

'2

∑∑∞=∞

=+-++-=

Ωm m

m m m m m

i u N

m q u N

m q q ππ

;42??????=N N ;221N N N -???

???=

;)/1)(1(2'2'k k v i i i Ω-Ω-= 12

'1

21

2,1;12N i v i i i =Ω+=

--α 22'22,1;12N i v i

i

i =Ω+=

β

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

;22i i i A αα+-=

;22i

i

i B ββ+-=

1221,1;1)(N i Z A Z A Z H i i i =++=∏-- 22

21

2,1;1)(N i Z

B Z B Z Z H i i i =++=∏--- )];

()([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: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

N n xx

N m n m x n x N

m R 10*1),

()(1)(?

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

)

()0(*1

,2i R a R k

i i k k

∑=+=σ

∑==-+=k

i k i k k a i k R a D 0

0,,0),

1(

2

1k k

k D σγ=

+

22

121)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时为止。将算出的模型参数代入下式即可得到功率谱估计值：

2

1,21)(?∑=-+=p

i jwi

i p p

j xx

e a e S σω

程序：

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

∑∑-=-

-+--=--+--+-=

1

2

1211

*

11)

)1()(()1()(2N k

n k k N k n k k k n e n e

n e n e γ

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

???

? ??-???? ??--=???? ??--+

--+

)1()(11)()(11*n e n e n e n e k k k

k k k γγ 1,,...,3,2,1,0,1,1,1,==-=----k i k k k i k i k a k i a a a γ

(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 模型的参数就全部估计出来了，根据以下公式即可算出功率谱：

2

1

2

1211)(?∑∑=-=-++=p

i jwi

i q

i jwi

i p

j xx e a e b e S σω

程序：

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)构造自相关阵

?????

?

? ??-------=)0()2()1())2(()0()

1())1(()1()

0(R N R N R N R R R N R R R

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

(2)计算R 的特征向量i v

?，i =1,2,…,N 。 (3)估计R 的维数M ，方法有AIC 和MDL 法。

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

∑+==

N

M i f P 12

MUSIC ?1

)(?i H

v

e

程序：

function S=music(X,p,ii) N=length(X);

r=xcorr(X,'biased'); clear R for k=1:N

R(:,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);

for k=1:512

S(k)=10*log10(1/(S0(k,:)*S0(k,:)'));

end

S=[S(257:512) S(1:256)];

subplot(2,2,ii);

plot(f,S,'b');

title(['MUSIC: ',int2str(p),' 维']);

xlabel('归一化频率'), ylabel('功率谱/dB');

hold on;

实验结果：

图7.MUSIC算法