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;
end;%endofif
end;%endofwhile
plot（accumulate_error,'m'）;
grid;
xlabel（'学习次数'）
ylabel（'误差'）
disp（['计算误差=',num2str（accumulate_error（circle_time））]）;
disp（['迭代次数=',num2str（circle_time）]）;
%测试
a0=double（[00]'）;
n1=w1*a0+b1;
a1=Logistic（n1）;
n2=w2*a1+b2;
a2=Logistic（n2）;
a=a2;
disp（['00=',num2str（a）]）;
a0=double（[01]'）;
n1=w1*a0+b1;
a1=Logistic（n1）;
n2=w2*a1+b2;
a2=Logistic（n2）;
a=a2;
disp（['01=',num2str（a）]）;
a0=double（[10]'）;
n1=w1*a0+b1;
a1=Logistic（n1）;
n2=w2*a1+b2;
a2=Logistic（n2）;
a=a2;
disp（['10=',num2str（a）]）;
a0=double（[11]'）;
n1=w1*a0+b1;
a1=Logistic（n1）;
n2=w2*a1+b2;
a2=Logistic（n2）;
a=a2;
disp（['11=',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）;
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、学习曲线图：
图1MLP
二、试用用奇阶互补法设计两带滤波器组（高、低通互补），进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：
Fp＝1.7KHz,Fr＝2.3KHz,Fs＝8KHz,Armin≥70dB。
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:
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
fori=1:
N1
alpha（i）=2*V（2*i-1）/（1+W（2*i-1）^2）;
end
fori=1:
N2
beta（i）=2*V（2*i）/（1+W（2*i）^2）;
end
fori=1:
N1
a（i）=（1-alpha（i）*Wc+Wc^2）/（1+alpha（i）*Wc+Wc^2）;
end
fori=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））;
forn=1:
length（w）
z=exp（j*w（n）*2*pi）;
H1=1;
fori=1:
N1
H1=H1*（a（i）+z^（-2））/（1+a（i）*z^（-2））;
end
H2=1/z;
fori=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'）;
%holdon;
xlabel（'数字频率'）;
ylabel（'幅度'）;
3、实验结果：
图2两带滤波器
4、四带滤波器组程序：
functionfilterfour
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
b=0;
form=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
fori=1:
N1
alpha（i）=2*V（2*i-1）/（1+W（2*i-1）^2）;
end
fori=1:
N2
beta（i）=2*V（2*i）/（1+W（2*i）^2）;
end
fori=1:
N1
a（i）=（1-alpha（i）*Wc+Wc^2）/（1+alpha（i）*Wc+Wc^2）;
end
fori=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））;
forn=1:
length（w）
z=exp（j*w（n）*2*pi）;
H1=1;
fori=1:
N1
H1=H1*（a（i）+z^（-2））/（1+a（i）*z^（-2））;
end
H21=1;
fori=1:
N1
H21=H21*（a（i）+z^（-4））/（1+a（i）*z^（-4））;
end
H2=1/z;
fori=1:
N2
H2=H2*（b（i）+z^（-2））/（1+b（i）*z^（-2））;
end
H22=1/（z^2）;
fori=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）、由输入数据估计自相关函数，一种渐近无偏估计（称之为取样自相关函数）的公式为：
（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;
end
%设置迭代初值
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）];
end
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；
实验结果：
图4Levinson算法
2、Burg算法
分析：
（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=[];
fork=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=[1a];
sigma2=real（sigma2）;
实验结果：
图5Burg算法
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）
N=length（signal_source）;
M=N-1;
r=xcorr（signal_source,'unbiased'）;
fork=1:
p
R（:
k）=（r（[N+q-k+1:
N+M-k]））.';
end
a1=（-pinv（R）*（r（[N+q+1:
N+M]）.'））.';
a=[1a1];
Y=filter（[1a1],[1],signal_source）;
pp=4*q;
[sigma,a1]=burg（Y,pp）;
sigma2=sigma（2:
end）;
[sigma,a2]=burg（a1（2:
end）,q）;
b=a2;
实验结果：
图6ARMA算法
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
fork=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）;
fork=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'）;
holdon;
实验结果：
图7MUSIC算法
升级会员 