%---------西安电子科技大学ADSP仿真作业第一题-----%
%------本题目要填写四个表格，每个表格附加一个图像--%
%------填写表格的时候，要改动四个代码块，在程序中标出来了--------%
%--表格1和表格3的填写需要取消注释奇数代码块1和3，注释掉代码块2和4-%
%--表格2和表格4的填写需要取消注释奇数代码块2和4，注释掉代码块1和3-%
%-------自定义数据集可以根据自己的喜好随意设置----%
%---------作者：小何-------------------------%
%---------请关注【科研果园】微信公众号----------%
%---------打开微信搜索"outputing"------------%
clear all;
clc;
%---------自定义的数据，a1,a2,Vw,M=信号长度----%
a1=0.1997;
a2=0.313;
Vw=0.19970313;
M=10000;
Vv_set=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0];
N_set=[1 2 3 4 5 6 7 8 9 10];
%--------------输出的预定义数据集合------------%
SNR_data=zeros(1,10);
cost_data=zeros(1,10);
cost_1_data=zeros(1,10);
for i=1:10
%-代码块1：第1个表和第3个表，用这个，注释掉下面的代码块2---%
% Vv=Vv_set(i);N=8;
%-代码块2：第2个表和第4个表，用这个，注释掉上面的代码块1---%
N=N_set(i);Vv=0.5;
%------------------这下面的代码不要改动--------------%
% 产生发射信号S(n)
randn('state',0);
W=sqrt(Vw)*randn(1,M);
S=zeros(1,M);
S(1)=W(1);
S(2)=W(2);
for n=3:M
S(n)=a1*S(n-1)+a2*S(n-2)+W(n);
end
%figure;n=1:M;stem(n,S);
%title('The source signal S(n)');xlabel('n');ylabel('S(n)');
%-----------产生噪声信号X(n)=S(n)+V(n)--------%
randn('state',1);
V= Vv*randn(1,M);
X=S+V;
%figure;n=1:M;stem(n,X(n));
%title('The received signal X(n)');xlabel('n');ylabel('X(n)');
%-----------计算信噪比SNR--------------------%
S_Power=mean(S.^2);
N_Power=mean(V.^2);
SNR=10*log10(abs(S_Power/N_Power));% SNR->信噪比
%-----------构造长度为N的Wiener滤波器----------%
r=zeros(1,N+1);
for n=0:N
temp=0;
for c=(n+1):M
temp=temp+X(c)*X(c-n);
end
r(n+1)=(1/M)*temp;
temp=0;
end
%r(1)=(1/M)(X(1)*X(1)+(X(2)*X(2)+(X(3)*X(3)+...+(X(M)*X(M))
%r(2)=(1/M)(X(2)*X(1)+(X(3)*X(2)+(X(4)*X(3)+...+(X(M)*X(M-1))
%r(3)=(1/M)(X(3)*X(1)+(X(4)*X(2)+(X(5)*X(3)+...+(X(M)*X(M-2))
%r(4)
%r(5)
%...
%r(n+1)=(1/M)(X(n+1)*X(1)+(X(n+2)*X(2)+(X(n+3)*X(3)+...+(X(M)*X(M-n))
%R为信号X(n)的自相关矩阵，维数为(N-1)*(N-1)
%R是由r构造而成的N*N的矩阵
R=zeros(N,N);
R(1,:)=r(1:N);
for row=2:N
R(row,1)=r(row);
R(row,2:N)=R(row-1,1:N-1);
end
%-----------P是X(n)与S(n)的互相关---------%
P=zeros(1,N);
X_temp=X;
for k=0:N-1
X_temp(1:M-k)=X(1+k:M);
X_temp(M-k+1:M)=0;
P(k+1)=1/M*X_temp*S';
end
Wiener=zeros(1,N);
Wiener=inv(R)*P';
Wiener=Wiener';
%----由Wiener滤波器进行滤波处理，得到Y(n)---%
Y_temp=conv(Wiener,X);
Y=Y_temp(1:M);
%-----------计算代价函数cost-------------%
E=zeros(1,M);
E=S-Y;
cost=mean(E.^2);
%-----------一步线性预测----------------%
N_P=N;
A_P=zeros(N_P);
r_p=zeros(1,N_P);
r_p=r(2:N+1);
A_P=A_P';
A_P=inv(R)*(r_p');
Y_P_temp=conv(A_P,X);
Y_P=zeros(1,M);
Y_P(3:M)=Y_P_temp(1:M-2);
%-----------计算一步线性预测后的cost-----%
E_P=zeros(1,M);
E_P=S-Y_P;
cost_1=mean(E_P.^2);
SNR_data(i)=SNR;
cost_data(i)=cost;
cost_1_data(i)=cost_1;
end
%-代码块3：第1个表和第3个表，用这个，注释掉下面的代码块4---%
% disp('a1=');disp(a1);
% disp('a2=');disp(a2);
% disp('N=');disp(N);
% disp('Var_v=');disp(Vv_set);
% disp('SNR=');disp(SNR_data);
% disp('cost_set=');disp(cost_data);
% disp('cost_1_set=');disp(cost_1_data);
% figure;
% plot(SNR_data,cost_data);title('The relationship between cost funtion & SNR');
% xlabel('SNR(dB)');ylabel('The cost');
% figure;
% plot(SNR_data,cost_1_data);title('The relationship of SNR & cost after step1 prediction');
% xlabel('SNR(dB)');ylabel('One cost after one step prediction');
%------------------------------------------%
%-代码块1：第2个表和第4个表，用这个，注释掉上面的代码块3---%
disp('a1=');disp(a1);
disp('a2=');disp(a2);
disp('Var_v=');disp(Vv);
disp('Var_w=');disp(Vw);
disp('N=');disp(N_set);
disp('SNR=');disp(SNR);
disp('cost_set=');disp(cost_data);
disp('cost_1_set=');disp(cost_1_data);
figure;
stem(N_set,cost_data);
title('The relationship of Filter length N & the cost');
xlabel('Filter length N');ylabel('The cost');
figure;
stem(N_set,cost_1_data);
title('The relationship of Filter length N & the cost_1');
xlabel('Filter length N');ylabel('The cost_1');
%------------------------------------------%
网友评论