chapter5_7.m

阵列信号处理中各中方位来波方向估计

%%%%chapter5_7%%%%%
%(1)当信号为相干信号时用music方法和最小模方法求解
%(2) 再用前后向空间平滑法求解
%由仿真结果可知当两信号为相干信号时，music方法和最小模方法失效，用前后向平滑技术估计出来的空间
%谱效果比较好
clc;clear;
P=2;%信号源个数
fai=[30,40]*pi/180;%信号的方向转为弧度
fai=pi*sin(fai);
theta=linspace(-100,100,1000);
L=8;
%%%产生A_theta%%%%
A_theta=zeros(L,2);
for k=1:L
    A_theta(k,1)=exp(j*fai(1)*(k-1));
    A_theta(k,2)=exp(j*fai(2)*(k-1));
end
%%%% 产生信号Sig%%%%
Sig=zeros(2,100);
for k=1:100
    Sig(1,k)=exp(j*pi/4*(k-1));
    Sig(2,k)=exp(j*pi/4*(k-1));%两信号为相干信号时
end
SNR=5;%信噪比为5db
R_guji=zeros(L,L);
y_gce=zeros(L,100);%实际的观测值
sums=zeros(L,L);
randn('state',sum(100*clock))
for k=1:100
    y_gce(:,k)=1/sqrt(2)*sqrt(10.^(SNR/10))*A_theta*Sig(:,k)+randn(L,1);%SNR(i)分别代表三个信噪比的取值
    sums=sums+y_gce*y_gce';
end
R_guji=sums/100;
%%%对估计的相关函数进行svd分解求出特征值，求得的u和v是一样的特征向量
[u,s,v]=svd(R_guji);
%%%求music空间谱估计%%%%%
sums=zeros(L,L);
for k=3:L
    sums=sums+u(:,k)*v(:,k)';
end  
 theta=linspace(-100,100,1000);  
 theta1=theta*pi/180;
 a_theta=zeros(L,1);
 for k=1:1000
     for n=1:L
         a_theta(n,1)=exp(j*pi*sin(theta1(k))*(n-1));
     end
     P_mu(k)=1/((a_theta)'*sums*a_theta);
 end
figure
 plot(theta,abs(P_mu))
  grid on
xlabel('角度');
ylabel('幅度');
title('music方法')
 %%%%%用最小模方法求解%%%%
 En=u(:,P+1:L);%En为噪声子空间，取得是3:8列
 En_=u(2:L,P+1:L);
 c=u(1,P+1:L).';%c为En的第一行元素
  d=zeros(L,1);%d为噪声子空间中最小模的一个向量;
  d(1,1)=1;
  d(2:L)=En_*conj(c)/(c'*c);
 a_theta=zeros(L,1);
 for k=1:1000
     for n=1:L
         a_theta(n,1)=exp(j*pi*sin(theta1(k))*(n-1));
     end
     P_min(k)=1/((a_theta)'*d*d'*a_theta);
 end
  figure

plot(theta,abs(P_min))
 grid on
xlabel('角度');
ylabel('幅度');
title('最小模方法')

%%%%%用前后向空间平滑法求解%%%%
%%%%(2)当M=6时%%%
M=4;%M可取3 4 5 6 7
child_array=L-M+1;%子阵个数；
% A_thetaM=zeros(M,P);
% for k=1:M
%     A_thetaM(k,1)=exp(j*fai(1)*(k-1));
%     A_thetaM(k,2)=exp(j*fai(2)*(k-1));
% end
% %%求D矩阵%%%%
% D=eye(2,2);%维数为pxp维
% for n=1:2
%     D(n,n)=exp(j*fai(n));
% end
% %%求J%%%%%
% J=eye(M,M);
% J=J(:,end:-1:1);%J为倒序矩阵；
%%%前后向的相关矩阵的估计值R_forw和R_back%%%%%
U_forw=zeros(M,child_array);%U每次采样的各阵子信号
R_forw=zeros(M,M);
U_back=zeros(M,child_array);%U每次采样的各阵子信号
R_back=zeros(M,M);
R_fb_guji =zeros(M,M);
for k=1:100
    i=0;
    for n=1:child_array
%         randn_noise=1/sqrt(2)*(randn(M,L-M+1)+j*randn(M,L-M+1));
%      randn_noise=randn(M,L-M+1);
%         U_forw(:,n)=1/sqrt(2)*sqrt(10.^(SNR/10))*A_thetaM*D.^(n-1)*Sig(:,k)+randn_noise(:,n);
%         R_forw=R_forw+U_forw(:,n)*U_forw(:,n)';%前向平滑的相关矩阵
%        U_back(:,n)=J*(1/sqrt(2)*sqrt(10.^(SNR/10))*conj(A_thetaM)*...
%        conj(D.^(child_array-n))*conj(Sig(:,k))+conj(randn_noise(:,child_array+1-n)));
%         R_back= R_back+U_back(:,n)*U_back(:,n)';%后向平滑的相关矩阵
 U_forw(:,n)=y_gce(1+i:M+i,k);
 R_forw=R_forw+U_forw(:,n)*U_forw(:,n)'
 U_back(:,n)=conj(y_gce(L-i:-1:child_array-i,k));
 R_back= R_back+U_back(:,n)*U_back(:,n)';%后向平滑的相关矩阵
  i=i+1;
    end
    R_fb_guji= R_fb_guji +(R_forw+R_back)/(200*child_array);
end

[u,s,v]=svd(R_fb_guji)%u，v两酉阵是相等的；
%%前后向平滑的music方法求解%%%%%
sums_M=zeros(M,M);
for k=3:M
   sums_M=sums_M+u(:,k)*v(:,k)';
end  
 theta=linspace(-100,100,1000);  
 theta1=theta*pi/180;
 a_theta=zeros(M,1);
 for k=1:1000
     for n=1:M
         a_theta(n,1)=exp(j*pi*sin(theta1(k))*(n-1));
     end
     P_mu(k)=1/((a_theta)'*sums_M*a_theta);
 end
 figure
 plot(theta,abs(P_mu))
 grid on
xlabel('角度');
ylabel('幅度');
title('前后向平滑的music方法')
%%%%前后向平滑的最小模方法求解%%%%%
 En_fb=u(:,P+1:M);%En_fb为噪声子空间，取得是3:M列
 Enfb_=u(2:M,P+1:M);
 c=u(1,P+1:M).';%c为En_fb的第一行元素
  d=zeros(M,1);%d为噪声子空间中最小模的一个向量;
  d(1,1)=1;
  d(2:M)=Enfb_*conj(c)/(c'*c);
  a_theta=zeros(M,1);
 for k=1:1000
     for n=1:M
         a_theta(n,1)=exp(j*pi*sin(theta1(k))*(n-1));
     end
     Pfb_min(k)=1/((a_theta)'*d*d'*a_theta);
 end
figure
plot(theta,abs(Pfb_min))
 grid on
xlabel('角度');
ylabel('幅度');
title('前后向平滑的最小模方法')