互质阵列中稀疏表示理论完成DOA估计算法_matlab

clc; clear all; close all; %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 初始化 %%%%%%%%%%%%%%%%%%%%%%%% lambda=1; d=lambda/2; % 阵元间距离，取为入射波长的一半 N=500; % 采样快拍数 DOA=[-45,-35,-25,0,10,20,40]; % DOA=[-20,10]; SignalNum=length(DOA); % 入射信号数量 ULA_ratio1=3; % 线阵1间距与d的倍率 ULA_ratio2=4; % 线阵2间距与d的倍率 M1=ULA_ratio2; % 线阵1阵元数量 M2=2*ULA_ratio1-1; % 线阵2阵元数量 M=M1+M2; % 总阵元数 SNR=0; Fc=[2*10^3,3*10^3,4*10^3,5*10^3,6*10^3,7*10^3,8*10^3]; % 入射信号频率 fs=20*10^3; % Fc=[2*10^3,2*10^3,2*10^3,2*10^3,2*10^3,2*10^3,2*10^3]; % 相干入射信号频率 thetatest=(-90*pi/180:1*pi/180:90*pi/180); % theta角度搜索范围 thetanum=length(thetatest); %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%信号生成%%%%%%%%%%%%%%%%%%%%%%% T_Vector=(1:N)/fs; mm=[0 3 4 6 8 9 12 16 20]; A1=zeros(M1,SignalNum); A2=zeros(M2,SignalNum); A3=zeros(length(mm),SignalNum); SignalVector1=zeros(SignalNum,N); SignalVector2=zeros(SignalNum,N); SignalVector3=zeros(SignalNum,N); %%%%%%构造线阵1信号 % for k1=1:SignalNum % A1(:,k1)=exp(-1j*(0:M1-1)'*2*pi*d*sin(DOA(k1)*pi/180)/lambda); % SignalVector1(k1,:)=exp(1j*2*pi*Fc(k1).*T_Vector); %信号 % end % Xt1=A1*SignalVector1; % % %%%%%%构造线阵2信号 % for k2=1:SignalNum % A2(:,k2)=exp(-1j*(1:M2)'*2*pi*d*sin(DOA(k2)*pi/180)/lambda); % SignalVector2(k2,:)=exp(1j*2*pi*Fc(k2).*T_Vector); %信号 % end % Xt2=A2*SignalVector2; %%%%%%%构造非均匀线阵信号 for k3=1:SignalNum A3(:,k3)=exp(-1j*mm'*2*pi*d*sin(DOA(k3)*pi/180)/lambda); SignalVector3(k3,:)=exp(1j*2*pi*Fc(k3).*T_Vector); %信号 end Xt3=A3*SignalVector3; %%%% 信号加噪 % for m1=1:M1 % Xt1(m1,:)=awgn(Xt1(m1,:),SNR,'measured'); % end % for m2=1:M1 % Xt2(m2,:)=awgn(Xt2(m2,:),SNR,'measured'); % end for m3=1:length(mm) Xt3(m3,:)=awgn(Xt3(m3,:),SNR,'measured'); end %%%%%%%%%% % %%%%%%%%%%%%%%% 原文方法 S_MUSIC 算法 %%%%%%%%%%%%%%%%%%%%%%% %%% R=Xt3*Xt3'/N; %%%% Vecter_R=R(:); %% -M1*M2:1:M1*M2 字典对应信息 %%%%%%%%%% r=Vecter_R([7,54,35,6,5,53,4,14,3,2,22,12,1,20,30,10,19,38,28,69,37,46,67,78,55]); MM=(length(r)+1)/2; Rx1=zeros(MM,MM); RX=cell(1,MM); for t=1:MM RX{1,t}=r(t:MM-1+t)*r(t:MM-1+t)'; end for tt=1:MM Rx1=Rx1+RX{1,tt}; end Rx=Rx1/MM; %%%%%%%%%%-----特征值分解----%M%%%%%%%%%%% [V, D] = eig(Rx); % X*V = V*D DD = diag(D); % 对角阵变矢量 [DD idx] = sort(DD, 'descend'); % 按从大往小排序特征值 Un = V(:, 1:MM-SignalNum); % 噪声子空间 Us = V(:, MM-SignalNum+1 : end); % 信号子空间 e1=[1,zeros(1,MM-1)].'; sigm_n=min(DD); % 最小特征值作为噪声的估计 sigm_n1=mean(DD((SignalNum+1:M),1)); % 特征值均值 I=eye(M); Pmusic=zeros(length(thetatest),1); for t=1:length(thetatest) aa=exp(1j*2*pi*sin(thetatest(t))*(0:MM-1)'*d/lambda); Pmusic(t)=1/abs((aa'*Un)*(Un'*aa)); end %%% %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %%%%%%%%%%%%%%%%%%%% 改进算法 稀疏求解 %%%%%%%%%%%%%%%%%%%%%%% %%%%%% %%%%%% %%%%计算对应字典的接收矢量 Position=[1 0 0 1 1 0 1 0 1 1 0 0 1 0 0 0 1 0 0 0 1 ]; % 阵元位置信息 for p=0:length(Position)-1 % 平移次数 temp=bitand(Position(p+1:end),Position(1:end-p)); for q=1:length(temp) if(temp(q)==1) qq=length(find(Position(1:q)==0)); % 平移阵列前面0的个数 pp=length(find(Position(1:q+p)==0)); % 固定线阵前面0个数 rr1(p+1)=Xt3(p+q-pp,:)*Xt3(q-qq,:)'/N; break; end end end R=Xt3*Xt3'/N; Vecter_R=R(:); %%% -M1*M2:1:M1*M2 字典对应信息 %%%%%%%%%% rr1=rr1(end-8:-1:1).'; %% 稀疏求解 for k4=1:thetanum % Dict(:,k4)=exp(-1j*(MM-1:-1:-MM+1)'*2*pi*d*sin(thetatest(k4))/lambda); Dict(:,k4)=exp(-1j*(MM-1:-1:0)'*2*pi*d*sin(thetatest(k4))/lambda); end % rr(MM) beita=sqrt(rr1(MM)^2-norm(rr1,2)^2/length(rr1)); cvx_begin variable bb(thetanum,1); minimize(norm(bb,1)); subject to Temp=rr1-Dict*bb; norm(Temp,2)<=beita; cvx_end %%%%%% %%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure; plot(thetatest.*180/pi,Pmusic,'b'); title('S-MUSIC'); figure; plot(thetatest.*180/pi,bb); title('稀疏表示'); beita