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

程序可运行,有图有真相，MATLAB得事先装好cvx优化包。clc;clear;close;lambda=1;d=lambda/2;%阵元间距离，取为入射波长的一半K=500;%采样快拍数theta=[-510];%入射角度SignalNum=length(theta);%入射信号数量Nnum=5;%%阵列阵元数量SNR1=-10;%%信噪比Aratio=sqrt(10^(SNR1/10));%信号幅度与噪声幅度比值，并假设信号幅度为1Fs=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矩阵fork2=1:SignalNumfork1=1:Nnum%1:12At(k1)=exp(j*(k1-1)*2*pi*d*sin(theta(k2)*pi/180)/lambda);A(k1,k2)=At(k1);endend%%%构造信号矩阵和噪声矩阵fork1=1:SignalNumSignalVector(k1,:)=exp(j*2*pi*Fc(k1).*T_Vector);%信号endXtt=A*SignalVector;%NoiseVector=sqrt(0.5)*(randn(Nnum,K)+j*randn(Nnum,K));forkk=1:NnumXt(kk,:)=awgn(Xtt(kk,:),SNR1,'measured');endRx=(Xt*Xt')./K;Rs=(SignalVector*SignalVector')./K;sigm_s=Rs(:,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%-----特征值分解----%M%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[V,D]=eig(Rx);%X*V=V*DDD=diag(D);%对角阵变矢量%[DDidx]=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);fork1=1:thetanumAtemp0=exp(j*2*pi*d/lambda*sin(thetatest(k1))*[0:Nnum-1]).';S(k1)=1/(Atemp0'*Un*Un'*Atemp0);endfigure(1)plot(thetatest.*180./pi,10*log10(abs(S)/(max(abs(S)))));%输出功率（dB）gridon;gridon;title('Music')xlabel('方位角（度）')ylabel('输出功率（dB）')%%%%%%%%%%%%%%%%%%%%%%%%%%----构造--G--selection矩阵%%%%%%%%%%%%%%M=Nnum;G=zeros(M*M,2*M-1);J0=eye(M);G(:,M)=J0(:);%fork=1:M-1fori=1:M-1J=[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);fork2=1:thetanum%相当于文章thita1---thitaQfork1=1:MBt(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';endend%%%----u---K稀疏矢量----u=zeros(1,thetanum);forz=1:SignalNumu(1,theta(z)+90+1)=sigm_s(z);%应该是等于sigm^2,每个信号的噪声方差？？？？endu=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*Mcvx_beginvariableu2(181,1)minimize(norm(u2,1))subjecttoS=W12*(y-sigm_n*I(:))-Q*u2;%S=y-G*Bthita*u2-sigm_n*I(:);norm(S)=beita;cvx_endfigure()plot(-90:1:90,u2);%-----------------------------------------------------------------------