資源描述:
《基于壓縮感知的doa估計(jì)程序》由會(huì)員上傳分享,免費(fèi)在線閱讀,更多相關(guān)內(nèi)容在工程資料-天天文庫(kù)。
1、程序可運(yùn)行,有圖有真相,MATLAB得事先裝好cvx優(yōu)化包。clc;clear;close;lambda=1;d=lambda/2;%陣元間距離,取為入射波長(zhǎng)的一半K=500;%采樣快拍數(shù)theta=[-510];%入射角度SignalNum=length(theta);%入射信號(hào)數(shù)量Nnum=5;%%陣列陣元數(shù)量SNR1=-10;%%信噪比Aratio=sqrt(10^(SNR1/10));%信號(hào)幅度與噪聲幅度比值,并假設(shè)信號(hào)幅度為1Fs=5*10^3;%信號(hào)頻率Fc=[2*10^3,5*10^3,8*10^3];%入射信號(hào)頻率fs=20*10^3
2、;thetatest=(-90*pi/180:1*pi/180:90*pi/180);%theta角度搜索范圍thetanum=length(thetatest);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%計(jì)算信號(hào)協(xié)方差矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%T_Vector=(1:K)/fs;A=zeros(Nnum,SignalNum);SignalVector=zeros(SignalNum,K);%NoiseVector=zeros(Nnum,K);Xt=zeros(Nnum,K);%%構(gòu)
3、造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%%%構(gòu)造信號(hào)矩陣和噪聲矩陣fork1=1:SignalNumSignalVector(k1,:)=exp(j*2*pi*Fc(k1).*T_Vector);%信號(hào)endXtt=A*SignalVector;%NoiseVector=sqrt(0.5)*(randn(Nnum,K)+j*randn(Nnum,K));fo
4、rkk=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);%對(duì)角陣變矢量%[DDidx]=sort(DD,'descend');%按從大往小排
5、序特征值Un=V(:,1:Nnum-SignalNum);%噪聲子空間Us=V(:,Nnum-SignalNum+1:end);%信號(hào)子空間e1=[1,zeros(1,Nnum-1)].';sigm_n=min(DD);%最小特征值^作為的估計(jì)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,
6、10*log10(abs(S)/(max(abs(S)))));%輸出功率(dB)gridon;gridon;title('Music')xlabel('方位角(度)')ylabel('輸出功率(dB)')%%%%%%%%%%%%%%%%%%%%%%%%%%----構(gòu)造--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
7、,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%相當(dāng)于文章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))/lambd
8、a);Bthita(:,k2)=Bt';endend%%%----u---K稀疏矢量----u=zeros(