现代信号处理大型作业汇总
合集下载
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
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';
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)]);
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;
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
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);
%测试
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;
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));
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('幅度');
2013年现代信号处理大型作业
(选择1,2,3三道题)
一、请用多层感知器(MLP)神经网络误差反向传播(BP)算法实现异或问题(输入为 ,要求可以判别输出为0或1),并画出学习曲线。其中,非线性函数采用S型Logistic函数。
1、原理:
反向传播(BP)算法:
(1)、多层感知器的中间隐层不直接与外界连接,其误差无法估计。
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;
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)) ;
(2)计算 ,判断 是否等于1,即该互补滤波器是否为互补镜像滤波器
(3)计算相关系数
(4)互补镜像滤波器的数字实现
2、程序:
function filter2()
Fp=1700;Fr=2300;Fs=8000;
Wp=tan(pi*Fp/Fs);
Wr=tan(pi*Fr/Fs);
Wc=sqrt(Wp*Wr);
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;
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
xlabel('数字频率');
ylabel('幅度');
5、实验结果:
图3四带滤波器组
三、根据《现代数字信号处理》(姚天任等,华中理工大学出版社,2001)第四章附录提供的数据(pp.352-353),试用如下方法估计其功率谱,并画出不同参数情况下的功率谱曲线:
1)Levinson算法
2)Burg算法
3)ARMA模型法
4)MUSIC算法
1、Levinson算法
分析:
(1)、由输入数据估计自相关函数,一种渐近无偏估计(称之为取样自相关函数)的公式为:
(2)、根据估计所得到的自相关函数,用下面的迭代公式估算AR模型参数:
(3)、对于AR(p)模型,按以上述递推公式迭代计算直到 时为止。将算出的模型参数代入下式即可得到功率谱估计值:
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);
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
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);
程序:
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);
m=0;
%----------------------------------------------------------
function [a]= Logistic(n)
a = 1./(1+exp(-n));
%----------------------------------------------------------
%训练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) ;%规模
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');
(2)、反向传播算法:从后向前(反向)逐层“传播”输出层的误差,以间接算出隐层误差。分两个阶段:
正向过程:从输入层经隐层逐层正向计算各单元的输出
反向过程:由输出层误差逐层反向计算隐层各单元的误差,并用此误差修正前层的权值。
2、流程图:
3、程序:
%使用了3层结构,第二层隐藏层4个单元。2,3层都使用Logisitic函数。
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);
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;%隐藏层的单元数
1 0 = 0.965390
1 1 = 0.043374
5、学习曲线图:
图1MLP
二、试用用奇阶互补法设计两带滤波器组(高、低通互补),进而实现四带滤波器组;并画出其频响。滤波器设计参数为:Fp=1.7KHz,Fr=2.3KHz,Fs=8KHz,Armin≥70dB。
1、设计步骤:
(1)对Fp、Fr进行预畸
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;
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)]);
for i=1:train_num
%前向传播
a0 = double ( p(i,:)' );%第i行数据
n1 = w1*a0+b1;
a1 = Logistic(n1);%第一个的输出
n2 = w2*a1+b2;
a2 = Logistic(n2)Fra bibliotek%第二个的输出
a = a2;
%后向传播敏感性
e = t(i,:)-a;
temp=accumulate_error(circle_time);
s2 = F(a2)*e;%输出层delta值
s1 = F(a1)*w2'*s2;%隐层delta值
%修改权值
wd1 = alpha .* s1*a0';
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)]);
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;
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
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);
%测试
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;
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));
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('幅度');
2013年现代信号处理大型作业
(选择1,2,3三道题)
一、请用多层感知器(MLP)神经网络误差反向传播(BP)算法实现异或问题(输入为 ,要求可以判别输出为0或1),并画出学习曲线。其中,非线性函数采用S型Logistic函数。
1、原理:
反向传播(BP)算法:
(1)、多层感知器的中间隐层不直接与外界连接,其误差无法估计。
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;
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)) ;
(2)计算 ,判断 是否等于1,即该互补滤波器是否为互补镜像滤波器
(3)计算相关系数
(4)互补镜像滤波器的数字实现
2、程序:
function filter2()
Fp=1700;Fr=2300;Fs=8000;
Wp=tan(pi*Fp/Fs);
Wr=tan(pi*Fr/Fs);
Wc=sqrt(Wp*Wr);
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;
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
xlabel('数字频率');
ylabel('幅度');
5、实验结果:
图3四带滤波器组
三、根据《现代数字信号处理》(姚天任等,华中理工大学出版社,2001)第四章附录提供的数据(pp.352-353),试用如下方法估计其功率谱,并画出不同参数情况下的功率谱曲线:
1)Levinson算法
2)Burg算法
3)ARMA模型法
4)MUSIC算法
1、Levinson算法
分析:
(1)、由输入数据估计自相关函数,一种渐近无偏估计(称之为取样自相关函数)的公式为:
(2)、根据估计所得到的自相关函数,用下面的迭代公式估算AR模型参数:
(3)、对于AR(p)模型,按以上述递推公式迭代计算直到 时为止。将算出的模型参数代入下式即可得到功率谱估计值:
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);
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
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);
程序:
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);
m=0;
%----------------------------------------------------------
function [a]= Logistic(n)
a = 1./(1+exp(-n));
%----------------------------------------------------------
%训练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) ;%规模
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');
(2)、反向传播算法:从后向前(反向)逐层“传播”输出层的误差,以间接算出隐层误差。分两个阶段:
正向过程:从输入层经隐层逐层正向计算各单元的输出
反向过程:由输出层误差逐层反向计算隐层各单元的误差,并用此误差修正前层的权值。
2、流程图:
3、程序:
%使用了3层结构,第二层隐藏层4个单元。2,3层都使用Logisitic函数。
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);
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;%隐藏层的单元数
1 0 = 0.965390
1 1 = 0.043374
5、学习曲线图:
图1MLP
二、试用用奇阶互补法设计两带滤波器组(高、低通互补),进而实现四带滤波器组;并画出其频响。滤波器设计参数为:Fp=1.7KHz,Fr=2.3KHz,Fs=8KHz,Armin≥70dB。
1、设计步骤:
(1)对Fp、Fr进行预畸
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;
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)]);
for i=1:train_num
%前向传播
a0 = double ( p(i,:)' );%第i行数据
n1 = w1*a0+b1;
a1 = Logistic(n1);%第一个的输出
n2 = w2*a1+b2;
a2 = Logistic(n2)Fra bibliotek%第二个的输出
a = a2;
%后向传播敏感性
e = t(i,:)-a;