现代信号处理大型作业题目+答案Word文件下载.docx
《现代信号处理大型作业题目+答案Word文件下载.docx》由会员分享，可在线阅读，更多相关《现代信号处理大型作业题目+答案Word文件下载.docx（20页珍藏版）》请在冰豆网上搜索。
一、请用多层感知器（MLP）神经网络误差反向传播（BP）算法实现异或问题（输入为
其中，非线性函数采用S型Logistic函数。
1、原理：
反向传播（BP）算法：
（1）、多层感知器的中间隐层不直接与外界连接，其误差无法估计。
（2）、反向传播算法：
从后向前（反向）逐层“传播”输出层的误差，以间接算出隐层误差。
分两个阶段：
正向过程：
从输入层经隐层逐层正向计算各单元的输出
反向过程：
由输出层误差逐层反向计算隐层各单元的误差，并用此误差修正前层的权值。
2、流程图：
3、程序：
%使用了3层结构，第二层隐藏层4个单元。
2，3层都使用Logisitic函数。
%训练xor数据。
functionmlp（）
f=fopen（'
XOR.txt'
）;
A=fscanf（f,'
%g'
[3inf]）;
A=A;
p=A（1:
2,:
）'
;
%训练输入数据
t=A（3,:
%desireout
[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）;
while1
temp=0;
circle_time=circle_time+1;
fori=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;
%endoffor
ifaccumulate_error（circle_time）<
=threshold|circle_time>
3001%thenbreak;
%endofif
%endofwhile
plot（accumulate_error,'
m'
grid;
xlabel（'
学习次数'
）
ylabel（'
误差'
disp（['
计算误差='
num2str（accumulate_error（circle_time））]）;
迭代次数='
num2str（circle_time）]）;
%测试
a0=double（[00]'
n1=w1*a0+b1;
a1=Logistic（n1）;
n2=w2*a1+b2;
a2=Logistic（n2）;
a=a2;
00='
num2str（a）]）;
a0=double（[01]'
01='
a0=double（[10]'
10='
a0=double（[11]'
11='
m=0;
%----------------------------------------------------------
function[a]=Logistic（n）
a=1./（1+exp（-n））;
function[result]=F（a）
[r,c]=size（a）;
result=zeros（r,r）;
fori=1:
r
result（i,i）=（1-a（i））*a（i）;
end;
4、实验结果：
计算误差=0.0049993
迭代次数=2706
00=0.023182
01=0.963110
10=0.965390
11=0.043374
5、学习曲线图：
图1.MLP
二、试用用奇阶互补法设计两带滤波器组（高、低通互补），进而实现四带滤波器组；
1、设计步骤：
（1）对Fp、Fr进行预畸
（2）计算
，判断
是否等于1，即该互补滤波器是否为互补镜像滤波器
（3）计算相关系数
（4）互补镜像滤波器的数字实现
2、程序：
functionfilter2（）
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;
forjj=1:
M
a=0;
form=0:
5
a=a+（-1）^m*q^（m*（m+1））*sin（（2*m+1）*pi*jj/N）;
%Nisodd,u=j
end
a
b=0;
form=1:
b=b+（-1）^m*q^（m^2）*cos（2*m*pi*jj/N）;
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
N1
alpha（i）=2*V（2*i-1）/（1+W（2*i-1）^2）;
N2
beta（i）=2*V（2*i）/（1+W（2*i）^2）;
a（i）=（1-alpha（i）*Wc+Wc^2）/（1+alpha（i）*Wc+Wc^2）;
b（i）=（1-beta（i）*Wc+Wc^2）/（1+beta（i）*Wc+Wc^2）;
w=0:
0.0001:
0.5;
LP=zeros（size（w））;
HP=zeros（size（w））;
forn=1:
length（w）
z=exp（j*w（n）*2*pi）;
H1=1;
H1=H1*（a（i）+z^（-2））/（1+a（i）*z^（-2））;
H2=1/z;
H2=H2*（b（i）+z^（-2））/（1+b（i）*z^（-2））;
LP（n）=abs（（H1+H2）/2）;
HP（n）=abs（（H1-H2）/2）;
end
plot（w,LP,'
k'
w,HP,'
%holdon;
xlabel（'
数字频率'
ylabel（'
幅度'
3、实验结果：
图2.两带滤波器
4、四带滤波器组程序：
functionfilterfour
%Nisodd,u=j
LLP=zeros（size（w））;
LHP=zeros（size（w））;
HLP=zeros（size（w））;
HHP=zeros（size（w））;
H1=H1*（a（i）+z^（-2））/（1+a（i）*z^（-2））;
H21=1;
H21=H21*（a（i）+z^（-4））/（1+a（i）*z^（-4））;
H2=H2*（b（i）+z^（-2））/（1+b（i）*z^（-2））;
H22=1/（z^2）;
H22=H22*（b（i）+z^（-4））/（1+b（i）*z^（-4））;
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）;
plot（w,LLP,'
w,LHP,'
w,HLP,'
g'
w,HHP,'
b'
5、实验结果：
图3.四带滤波器组
三、根据《现代数字信号处理》（姚天任等，华中理工大学出版社，2001）第四章附录提供的数据（pp.352-353），试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：
1、Levinson算法
分析：
（1）、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：
（2）、根据估计所得到的自相关函数，用下面的迭代公式估算AR模型参数：
（3）、对于AR（p）模型，按以上述递推公式迭代计算直到
时为止。
将算出的模型参数代入下式即可得到功率谱估计值：
程序：
function[sigma2,a]=levinson（signal_source,p）
%阶数由p确定
N=length（signal_source）;
%确定自相关函数
form=0:
N-1
R（m+1）=sum（conj（signal_source（[1:
N-m]））.*signal_source（[m+1:
N]））/N;
%设置迭代初值
a1=-R
（2）/R
（1）;
sigma2=（1-abs（a1）^2）*R
（1）;
gamma=-a1;
%开始迭代
fork=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）];
a=[1a1];
%计算功率谱估计值
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个等间隔数据
fork=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）、设输入数据序列为
，对前后向预测误差之和求偏导，得反射系数
前后向预测误差递推公式：
（2）、重复以上步骤直至k=p，根据迭代得到的AR模型参数计算功率谱，计算功率谱的公式同上面算法。
function[sigma2,a]=burg（signal_source,p）
ef=signal_source;
eb=signal_source;
sigma2=sum（signal_source*signal_source'
）/N;
a=[];
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）];
a=[1a];
图5.Burg算法
3、ARMA算法
（1）、用x（n）通过A（z），得到y（n）。
（2）、用一无穷阶的AR模型近似MA模型。
用Burg算法可得到此近似AR模型的参数以及激励白噪声的功率。
一般此AR模型的阶数应大于MA模型阶数的两倍，以获得较好的近似效果。
（3）、可以证明，将上一步求出的近似AR模型参数视为时间序列，则MA模型就可视为一线性预测滤波器，按最小均方误差准则就可以求出MA模型参数，方法同样可用Burg算法。
这样，ARMA模型的参数就全部估计出来了，根据以下公式即可算出功率谱：
function[a,b,sigma2]=arma（signal_source,p,q）
M=N-1;
r=xcorr（signal_source,'
unbiased'
R（:
k）=（r（[N+q-k+1:
N+M-k]））.'
a1=（-pinv（R）*（r（[N+q+1:
N+M]）.'
））.'
Y=filter（[1a1],[1],signal_source）;
pp=4*q;
[sigma,a1]=burg（Y,pp）;
sigma2=sigma（2:
[sigma,a2]=burg（a1（2:
end）,q）;
b=a2;
图6.ARMA算法
4、MUSIC算法
（1）构造自相关阵
自相关函数可用有偏估计代替。
（2）计算R的特征向量
，i=1,2,…,N。
（3）估计R的维数M，方法有AIC和MDL法。
（4）根据剩余特征向量计算MUSIC谱。
functionS=music（X,p,ii）
N=length（X）;
r=xcorr（X,'
biased'
clearR
N
k）=（r（[N-（k-1）:
2*N-k]））.'
[V,D]=eig（R）;
f=linspace（-0.5,0.5,512）;
S0=fft（V（:
p+1:
end）,512）;
S（k）=10*log10（1/（S0（k,:
）*S0（k,:
））;
S=[S（257:
512）S（1:
256）];
subplot（2,2,ii）;
plot（f,S,'
title（['
MUSIC:
'
int2str（p）,'
维'
]）;
归一化频率'
）,ylabel（'
功率谱/dB'
holdon;
图7.MUSIC算法
升级会员 