www.pudn.com > twoDomp.zip > twoDomp.m, change:2014-01-03,size:3010b


%-------------参考文献:《基于空间角稀疏表示的二维DOA估计》-李鹏飞------------% 
%-----------------《基于压缩感知的二维DOA估计-陈玉龙》-----------------------% 
%----------------描述的是L型阵列的二维DOA估计问题,利用OMP算法---------------% 
%----------------------------------2014/1/2--------------------------------% 
close all; 
clear all; 
N=9;%阵元个数 
f=1.5e8; 
c=3.0e8; 
lmda=c/f; 
d=lmda/2; 
rad=pi/180; 
snap=128;%快拍数 
snr=25;%信噪比 
theta=[10  11];%俯仰角 
phi=[15  60];%方位角 
K=length(theta);%入射信号个数 
for i=1:N 
    for j=1:K 
        T1(j)=cos(theta(j)*rad)*cos(phi(j)*rad); 
        T2(j)=cos(theta(j)*rad)*sin(phi(j)*rad); 
        Ax(i,j)=exp(-1j*2*pi/lmda*(i-1)*d*T1(j));%X轴上的阵列流型矢量 
        Ay(i,j)=exp(-1j*2*pi/lmda*(i-1)*d*T2(j));%Y轴上的阵列流型矢量 
    end 
end 
S=randn(K,snap); 
X=awgn(Ax*S,snr);%X轴阵列接收到的数据(加噪声) 
Y=awgn(Ay*S,snr);%Y轴阵列接收到的数据(加噪声) 
[U1,L1,V1]=svd(X); 
[U2,L2,V2]=svd(Y); 
[L1,uu1]=sort(diag(L1),'descend'); 
[L2,uu2]=sort(diag(L2),'descend'); 
L1=diag(L1(1:K,:)); 
L2=diag(L2(1:K,:)); 
LL1=inv(L1); 
LL2=inv(L2); 
C=zeros(snap-K,K); 
Dk1=[LL1';C]; 
Dk2=[LL2';C]; 
X1=X*V1*Dk1;   
Y1=Y*V2*Dk2;%降维处理 
%% 对上面的X、Y进行稀疏重构 
%theta1=-90:1:90; 
%phi1=-90:1:90; 
%M=length(theta1); 
%for j=1:M 
%    T11(j)=cos(theta1(j)*rad)*cos(phi1(j)*rad); 
%    T12(j)=cos(theta1(j)*rad)*sin(phi1(j)*rad); 
%end 
% 完备原子库的构造 
T11=-1:1/100:1; 
T12=-1:1/100:1; 
MM=length(T11); 
for i=1:N 
    for j=1:MM 
        Ax1(i,j)=exp(-1j*2*pi/lmda*(i-1)*d*T11(j));%X轴上的阵列流型矢量 
        Ay1(i,j)=exp(-1j*2*pi/lmda*(i-1)*d*T12(j));%Y轴上的阵列流型矢量 
    end 
end 
% OMP算法    
S1=zeros(MM,1); 
S2=zeros(MM,1); 
AX=[]; 
AY=[]; 
p1=[]; 
p2=[]; 
rx=X1; 
ry=Y1; 
r=1; 
for r=1:K 
    yx=Ax1'*rx; 
    yy=Ay1'*ry; 
    for mm=1:size(yx(:,1)); 
        xx1(mm)=norm(yx(mm,:)); 
        yy1(mm)=norm(yy(mm,:)); 
    end 
    [val1 index1]=sort(xx1,'descend'); 
    [val2 index2]=sort(yy1,'descend'); 
    index_1=index1(1); 
    index_2=index2(1); 
    p1=[p1 index_1]; %X轴和Y轴接收信号的稀疏解用相同的结果,怎么表示alpha,bleta 
    p2=[p2 index_2]; 
    AX=Ax1(:,index_1); 
    AY=Ay1(:,index_2); 
     if r==1 
         qx(:,1)=AX/norm(AX); 
         qy(:,1)=AY/norm(AY); 
     else 
         bx=[];by=[]; 
         bx(:,1)=AX; 
         by(:,1)=AY; 
         for qq=2:r 
             bx(:,qq)=bx(:,qq-1)-qx(:,qq-1)'*bx(:,qq-1)*qx(:,qq-1); 
             by(:,qq)=by(:,qq-1)-qy(:,qq-1)'*by(:,qq-1)*qy(:,qq-1); 
         end 
         qx(:,qq)=bx(:,qq)/norm(bx(:,qq)); 
         qy(:,qq)=by(:,qq)/norm(by(:,qq)); 
     end 
    rx=rx-qx(:,r)*qx(:,r)'*rx; 
    ry=ry-qy(:,r)*qy(:,r)'*ry; 
    r=r+1;      
end 
  tt1=-1+1/100*(p1-1); 
  tt2=-1+1/100*(p2-1); 
  for i=1:K 
      %alpha(i)=acos(tt1(i)); 
      %bleta(i)=acos(tt2(i)); 
      phix(i)=atan(tt2(i)/tt1(i)); 
      pp1(i)=sqrt(tt1(i)^2+tt2(i)^2); 
      thetax(i)=acos(pp1(i)); 
  end 
  thetax0=thetax/rad;   %估计theta的值 
  phix0=phix/rad;       %估计phi的值