压缩感知doa matlab,基于压缩感知的DOA估计程序

基于压缩感知的DOA估计程序

程序可运行,有图有真相，MATLAB得事先装好cvx优化包。 clc; clear; close; lambda=1; d=lambda/2; %阵元间距离，取为入射波长的一半 K=500; %采样快拍数 theta=[-5 10]; %入射角度 SignalNum=length(theta); %入射信号数量 Nnum=5; %%阵列阵元数量 SNR1=-10; %%信噪比 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; u2=2; 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); %-----------------------------------------------------------------------