求根MUSIC算法

rand（‘seed‘0）;
randn（‘seed‘0）;
clear all;  
N=6;   %%The number of elements
dd=0.5;%%Equal to d/lamda
Snapshot=512;
SNR=30;%%Input SNR （dB）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta=linspace（-50503）;
% theta=[-30 -5 30];
D=length（theta）;  %%The number of the DOA

V_theta=exp（1j*2*pi*dd*（[0:N-1].‘*sin（theta*pi/180）））;

error_gain=0.8+0.4*rand（N1）;
error_gain=error_gain./error_gain（1）;
error_phase=exp（1j*2*pi*rand（N1））;
% error_phase=exp（1j*2*pi*rand（N1））;
error_phase=error_phase./error_phase（1）;
error_sensors=error_gain.*error_phase;
error_sensors=diag（error_sensors.‘）;%% Diagnol matrix

for jj=1:D
    for ii=1:Snapshot
        power_signal_source（jjii）=10*（randn（1）+1j*randn（1））/sqrt（2）;
    end
end
clear jj;clear ii;

% X=V_theta*S;
% X=V_theta*power_signal_source;

X=error_sensors*V_theta*power_signal_source;

X1=awgn（XSNR‘measured‘）;
% Rxxm=X*X‘/Snapshot;
Rxxm=X1*X1‘/Snapshot;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[U DD]=eig（Rxxm）;

Noise_space=U（:1:N-D）;
iteration=1;
F=eye（N）;
% 
for ii=1:iteration
    
   X1=inv（F）*X1;
   Rxxm=X1*X1‘/Snapshot;
   [U DD]=eig（Rxxm）; 
   Noise_space=U（:1:N-D）;
% 
   syms z;
   pz=z.^（[0:N-1]‘）;  
   pz1=（z^（-1））.^（[0:N-1]）; 
   fz=z^（N-1）*pz1*Noise_space*Noise_space‘*pz; 
   a=sym2poly（fz）;
   r=roots（a）;
   r1=abs（r）;
   for jj=1:2*D     
       [YI（jj）]=min（abs（r1-1））;     
       r1（I（jj））=inf;  
   end

   for jj=1:2*D      
       theta1（jj）=asin（angle（r（I（jj）））/pi）*180/pi;%
   end
   
   theta2=floor（theta1*10^5）/10^5;
   theta_rootMUSIC=unique（theta2）;
   
   w=zeros（N1）;
   w（11）=1;
   Qk=zeros（NN）;
   
   for jj=1:D
       V_theta_rootMUSIC=exp（1j*2*pi*dd*（[0:N-1].‘*sin（theta_rootMUSIC（jj）*pi/180）））;
       diag_V_theta_rootMUSIC=diag（V_theta_rootMUSIC.‘）;
       Qk=Qk+（diag_V_theta_rootMUSIC‘*Noise_space*Noise_space‘*diag_V_theta_rootMUSIC）;
   end
   
   delta_kPlus1=pinv（Qk）*w/（w‘*pinv（Qk）*w）;
   
   F=diag（delta_kPlus1.‘）;
end

% F=eye（N）;

% Search_theta=-90:0.2:90;
% Co_Music=zeros（1length（Search_theta））;
% 
% for n=1:length（Search_theta）
%     A_theta=exp（1j*2*pi*dd*（[0:N-1].‘*sin（Search_theta（n）*pi/180）））;
%     Co_Music（n）=1/（A_theta‘*F‘*Noise_space*Noise_space‘*F*A_theta）;
% %     Co_Music（n）=1/（A_theta‘*Noise_space*Noise_space‘*A_theta）;
% end
% 
% Co_Music=（abs（Co_Music）./max（abs（Co_Music）））;
% 
% figure（1）;
% plot（Search_thetaCo_Music）;
% hold on;
% axis（[-90 90 0 1]）;
% fontsize（12）; 
% vline（theta）;
% xlabel（‘\theta（deg）‘）;
% ylabel（‘Spectrum（Normalized）‘）;
% hold off;
% 
% Peak_theta=[-29.8 -4.2 34];
% % Peak_theta=[-50 -18 3.4 26 62.8];
% Peak_theta=[2.6 21.4 71.6];
% 
% Qk=zeros（NN）;
% for jj=1:D
%     V_Peak_theta=exp（1j*2*pi*dd*（[0:N-1].‘*sin（Peak_theta（jj）*pi/180）））;