www.gusucode.com > 协方差矩阵稀疏DOA-源码程序 > 协方差矩阵稀疏DOA/l1_SRACV.m

    clear all;
close all;
clc;
lambda=60;
dd=lambda/2;      %阵元间距离，取为入射波长的一半
Kp=80;            %采样快拍数
Nnum=9; %%阵列1阵元数量
fs=20*10^3;

theta=[6,26];
AA=[2,3,3];
SignalNum=length(theta);
Fc=[2*10^3,5*10^3,2*10^3];     %入射信号频率
SNR=10;  
Aratio=sqrt(10^(SNR/10));    %信号幅度与噪声幅度比值，并假设信号幅度为1
% thetatest=(0*pi/180:1*pi/180:90*pi/180);   %theta角度搜索范围
thetatest=[-60:1:60]*pi/180;
thetanum=length(thetatest);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SignalVector=zeros(SignalNum,Kp);
Xt=zeros(Nnum,Kp);
    M=Nnum;

    T_Vector=(1:Kp)/fs;
d=dd*(0:M-1);
w=2*pi/lambda*sin(thetatest);
A=exp(-1j*d'*w);

% % %  for k2=1:thetanum
% % %     for k1=1:M
% % %         Att(k1)=exp(-j*(k1-1)*2*pi*dd*sin(thetatest(k2))/lambda);
% % %         A(k1,k2) =Att(k1);
% % %     end
% % %   end

%%构造观测矩阵Y
for k2=1:SignalNum
    for k1=1:Nnum                          
        At(k1)=exp(-j*(k1-1)*2*pi*dd*sin(theta(k2)*pi/180)/lambda);
        Atemp(k1,k2)=At(k1);   
    end     
end
%%%构造信号矩阵和噪声矩阵
for k1=1:SignalNum
    SignalVector(k1,:)=exp(j*2*pi*Fc(k1).*T_Vector);  %信号
end
Xtt=Atemp*SignalVector;

%NoiseVector=sqrt(0.5)*(randn(Nnum,K)+j*randn(Nnum,K));
%--------------------加噪声方法---1
for kk=1:Nnum
   Xt(kk,:)=awgn(Xtt(kk,:),SNR,'measured');
end
%-------------------加噪声方法---2
%     for k1=1:Nnum
%         NoiseVector(k1,:)=SignalNum/Aratio.*noisecg(Kp)';
%         Xt(k1,:)=Xtt(k1,:)+NoiseVector(k1,:);
%     end 

 Rx=(Xt*Xt')./Kp;
 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[V, D] = eig(Rx);                %   X*V = V*D 
DD     = diag(D);               % 对角阵变矢量
sigm_n=min(DD);        %最小特征值^作为 的估计

I=eye(M);  
W12=sqrt(Kp)*(kron((Rx^(-0.5)).',Rx^(-0.5)));
Y=Rx-sigm_n*I;
y=W12*Y(:);

Fai=W12*kron(I,A);    

II=ones(1,thetanum);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
belta=sqrt(chi2inv(0.99999,M*M))

% % % cvx_begin
% % %     variable B(thetanum,M) 
% % %      variable b(thetanum,1);    
% % %       minimize(norm(b,1))
% % %     % minimize (norm(sqrt(sum((B.^2)')),1))
% % %    subject to    
% % %    for kk=1:1:thetanum
% % %          norm(B(kk,:))<=b(kk,1);
% % %    end   
% % %      temp0=Rx-A*B-sigm_n*I;
% % %      temp=W12*vec(temp0);
% % %       norm(temp)<=belta ;
% % % %    norm(y-Fai*B(:))<=23 ;
% % %  cvx_end
% % % 
% % %  figure()
% % %  plot(0:1:90,b);
%%%%%%%%%%%%%%%%%%%%%%----B,,W*vec(R-AB-I)---%%%%%%%%%%%%%%%%%%%%  
 belta=sqrt(chi2inv(0.9999,M*M));
  AP=[];
    for i=1:1:M
        AP=blkdiag(AP,A);
    end
    
% I2=eye(thetanum);
%  APP=[];
% cvx_begin
%     variable B(thetanum*M) complex;
%      variable b(thetanum);
%    
%       minimize(norm(b,1))
%     % minimize (norm(sqrt(sum((B.^2)')),1))
%    subject to  
%        for i=1:1:thetanum
%         for j=1:1:M
%             APP=blkdiag(APP,I2(i,:));
%         end
%         b(i)>=norm(APP*B);
%         APP=[];
%     end
%     
%         temp0=vec(Rx-sigm_n*I)-vec(AP*B);
%         temp=W12*temp0;
%           norm(temp)<=belta;
% cvx_end
% 
%  figure(1)
%  plot(0:1:90,b);
 %%%%%%%% %%%%%%%%%%%%%%%%%%----B,r,g---y-Fai*vec(B)---%%%%%%%%%%%%%%%%%%%%%%%%    
 
 I2=eye(thetanum);
 APP=[];
cvx_begin
    variable B(thetanum*M) complex;
    variable r1(thetanum);
    variable g
          minimize(g)
   subject to   
     norm(y-Fai*B)<=belta;
     
   for i=1:1:thetanum
        for j=1:1:M
            APP=blkdiag(APP,I2(i,:));
        end
        r1(i)>=norm(APP*B);
        APP=[];
    end
           II*r1<=g;        
cvx_end

%  figure(2)
%  plot(-90:1:90,r);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

% II=ones(1,thetanum);
%  cvx_begin
%   variable B(thetanum,M) complex;
%   variable r(thetanum);
%   variable g
%       minimize(g)
%  subject to   
%      norm(y-Fai*B)<=belta;
%      
%    for i=1:1:thetanum
%         for j=1:1:M
%             APP=blkdiag(APP,I2(i,:));
%         end
%         r(i)>=norm(APP*B);
%         APP=[];
%     end
%            II*r<=g;        
% cvx_end
%  
%   figure()
%   plot(-90:1:90,r);