www.pudn.com > SSR_DOA.rar > SSR_doa.m, change:2014-03-13,size:3436b


clc; 
clear; 
close; 
 
lambda=1; 
d=lambda/2;      %阵元间距离,取为入射波长的一半 
K=500;            %采样快拍数 
theta=[-50,40];  %入射角度   
SignalNum=length(theta);     %入射信号数量 
Nnum=5;                 %%阵列阵元数量 
SNR1=-8;          %%信噪比  
Aratio=sqrt(10^(SNR1/10));    %信号幅度与噪声幅度比值,并假设信号幅度为1 
 
Fs=5*10^3;               %信号频率 
Fc=[2*10^3,5*10^3,8*10^3];       %入射信号频率 
fs=20*10^3; 
 
thetatest=(-90*pi/180:1*pi/180:90*pi/180);   %theta角度搜索范围   
thetanum=length(thetatest); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算信号协方差矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
T_Vector=(1:K)/fs; 
A=zeros(Nnum,SignalNum); 
SignalVector=zeros(SignalNum,K); 
%NoiseVector=zeros(Nnum,K); 
Xt=zeros(Nnum,K); 
%%构造A矩阵 
for k2=1:SignalNum 
    for k1=1:Nnum                          %1:12 
        At(k1)=exp(j*(k1-1)*2*pi*d*sin(theta(k2)*pi/180)/lambda); 
        A(k1,k2)=At(k1);    
    end      
end 
%%%构造信号矩阵和噪声矩阵 
for k1=1:SignalNum 
    SignalVector(k1,:)=exp(j*2*pi*Fc(k1).*T_Vector);  %信号 
end 
Xtt=A*SignalVector; 
 
%NoiseVector=sqrt(0.5)*(randn(Nnum,K)+j*randn(Nnum,K)); 
for kk=1:Nnum 
   Xt(kk,:)=awgn(Xtt(kk,:),SNR1,'measured'); 
end 
Rx=(Xt*Xt')./K; 
 
Rs=(SignalVector*SignalVector')./K; 
sigm_s=Rs(:,1); 
% %%%%%%%%%%%%%%%%%%%%%%%%%%-----特征值分解----%M%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
[V, D] = eig(Rx);                %   X*V = V*D  
DD     = diag(D);               % 对角阵变矢量 
% [DD idx] = sort(DD, 'descend'); % 按从大往小排序特征值 
Un     = V(:, 1:Nnum-SignalNum);         % 噪声子空间 
Us     = V(:, Nnum-SignalNum+1 : end);  % 信号子空间 
e1=[1,zeros(1,Nnum-1)].';  
sigm_n=min(DD);        %最小特征值^作为 的估计 
I=eye(Nnum); 
 
 for k1=1:thetanum 
        Atemp0=exp(j*2*pi*d/lambda*sin(thetatest(k1))*[0:Nnum-1]).'; 
    S(k1)=1/(Atemp0'*Un*Un'*Atemp0);   
  end 
% figure(1) 
% plot(thetatest.*180./pi,10*log10(abs(S)/(max(abs(S)))));%输出功率(dB) 
% grid on; 
% % grid on; 
% title('Music') 
% xlabel('方位角(度)')                   
% ylabel('输出功率(dB)') 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%----构造--G--selection矩阵%%%%%%%%%%%%%% 
 
M=Nnum; 
G=zeros(M*M,2*M-1); 
J0=eye(M); 
G(:,M)=J0(:); 
 %for k=1:M-1 
    for i=1: M-1 
      J=[zeros(M-i,i),eye(M-i);zeros(i,i),zeros(i,M-i)]; 
      G(:,M-i)=J(:); 
      J1=J'; 
      G(:,M+i)=J1(:); 
    end 
 
 
%%%%%%----Bthita------%%%% 
 Bthita=zeros(2*M-1,thetanum); 
 Bt=zeros(1,2*M-1); 
for k2=1:thetanum                    %相当于文章thita1---thitaQ 
     for k1=1:M                          
         Bt(1,k1+M-1)=exp(-j*(k1-1)*2*pi*d*sin(thetatest(k2))/lambda);   
         Bt(1,k1)=exp(j*(M-k1)*2*pi*d*sin(thetatest(k2))/lambda);   
         Bthita(:,k2)=Bt';    
     end   
end 
%%%----u---K稀疏矢量---- 
u=zeros(1,thetanum); 
for z=1:SignalNum 
    u(1,theta(z)+ 90+1)=sigm_s(z);    %应该是等于sigm^2,每个信号的噪声方差???? 
     
end 
 u=u'; 
%%%%%%%%%%%%%%%%%%%%%%%%%%----cvx运算-------%%%%%%%%%%%%%%%%%%%%% 
y=Rx(:); 
 
W12=sqrt(K)*(kron((Rx^(-0.5)).',Rx^(-0.5))); 
Q=W12*G*Bthita;               
S1=W12*(y-sigm_n*I(:))-Q*u; 
 
beita=sqrt(chi2inv(0.999999,M*M));               %卡方分布M*M 
 cvx_begin 
 variable u2(181,1) 
 minimize( norm(u2,1)) 
 subject to 
 S=W12*(y-sigm_n*I(:))-Q*u2;       
%  S=y-G*Bthita*u2-sigm_n*I(:); 
 norm(S) <=beita ; 
 cvx_end 
 figure() 
 plot(-90:1:90,u2); 
  
%-----------------------------------------------------------------------