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

研究生“现代信号处理”课程大型作业
（以下四个题目任选三题做）
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 横向或格-梯型自适应均衡器
参考文献
[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、流程图：
开始
选择初始值
前向计算求所有神经元的输出
计算输出层
从后向前计算隐层
计算保存权值修正量
修正权值
N
是否收敛？
Y
结束
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)、设输入数据序列为1
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算法。