music_xnfx_snr.m

这是用蒙特卡罗方法计算music算法的误差分析

clear all;
clc;
T = 5000;
s=0;
N=256;                            %设定采样点数
M =8;                             %设定阵元数
fs=8;                             %设定采样频率
kd=1/2;                           %阵元间距与波长的比值
f=[4 6];                            %设两个信号元的频率，单位为MHz
theta = [-50 10]*pi/180;              %入射方向角向量
phi =[60 20]*pi*180;                 %两入射信号初始相位
for SNR = -5:20;
    RMS = zeros(1,2);
    for t = 1:T;  
        N=256;                            %设定采样点数
        M =8;                             %设定阵元数
        fs=8;                             %设定采样频率
        kd=1/2;                           %阵元间距与波长的比值
        AMP = [10.^(SNR./20) 10.^(SNR./20)];              
        m=0:M-1;
        n=0:N-1;
        S = repmat(AMP',1,N).*exp(j*(2*pi*f.'/fs*n+phi'*ones(1,N)));  %产生两个远场正弦信号
        A = exp(j*2*pi*kd*m'*sin(theta));                             %方向矩阵

        X=A*S;
        randn('seed',sum(100*clock));
        noise  = randn(M,N) + j*randn(M,N);
        X=X+noise;                                                    %接受数据矩阵

        Rx=X*(X')/N;                                                  %计算协方差矩阵
%*********************************************************************
       [V,D]=eig(Rx);			          %自相关矩阵的特征分解,,得到特征值矩阵和与特征值一一对应的特征向量矩阵
       [Q,index]=sort(diag(D));          %对特征值排序 

       for k=1:M-length(theta)          %将与特征值对应的特征向量也按照特征值的顺序排列     
           En(:,k)=V(:,index(k));        %取排了序的特征向量的前M-length(theta)个来构造噪声子空间En
       end
       L_lim=-60;                        %设置测向左、右界
       R_lim=60;
       step=0.1;                         %设置步长
       theta_axis=L_lim:step:R_lim;
       MD=(R_lim-L_lim)/step+1;
       a_theta=exp(j*2*pi*kd*m'*sin(theta_axis*pi/180));
       for k=1:MD                        %进行谱峰搜索
           d(k)=a_theta(:,k)'*En*En'*a_theta(:,k);
           k=k+1; 
       end
       D2=-10*log10(abs(d)./max(abs(d)));
       plot(theta_axis,D2,'k');
       ylim([0,50])
       grid


%************************************************
       if length(theta) == 1
          [Y, edoaIx] = max(D2);
          ang = theta_axis;
          edoa = ang(edoaIx);
          num = 1;
       else
          edoaIx = 0;
          l=0;
          jumplimit = 100;
          jump = 0;
          startTreshold = .5;
          treshold = startTreshold;
          D2max = max(D2);
    
          while ((length(edoaIx) ~= length(theta)) & (jump < jumplimit))
             edoaIx = 0;
             l=0;
%******************************************%Smoothing
             for k = 1 : length(D2),        
                if D2(k) < D2max*treshold,
                   Q(k)=0;
                else
                   Q(k)=D2(k);
                end
             end % for

             dQ = diff(Q);
             dQ(length(dQ)+1) = 0;
             sdQ = sign(dQ);
         
%**********************************************%Searching
             for k = 2:length(sdQ)-1
                if sdQ(k-1) == 1
                   while (sdQ(k)==0)
                       k=k+1;
                   end % while.
                   if sdQ(k) == -1
                       l=l+1;
                      edoaIx(l) = k;                 %记录峰值对应序号
                   end
                end
             end % for
%***************************************% 数量不同，调整门限
             if length(edoaIx) > length(theta) 
                treshold = treshold*1.1;
                if treshold > 1
                   jump = jumplimit;
                end

             elseif length(edoaIx) < length(theta)
                if treshold < startTreshold
                  jump = jumplimit;
                else
                  treshold = treshold/1.2;
                end
             end
             jump = jump + 1;

           end%while
       end%if
%***********************************************
       if (edoaIx ~= 0)
           ang = theta_axis;
           edoa = ang(edoaIx);
           num = length(edoaIx);
       else
           edoa=0;
           num = 0;
       end%if
       RMS = (edoa-theta*180/pi).^2 + RMS;
    end
    RMSE = sqrt(RMS/T);  
    s = s+1;
    R(s) = RMSE(1);
    K(s) = RMSE(2);
end
SNR=-5:20;
plot(SNR,R,'k*-',SNR,K,'ro-')
xlabel('SNR/(dB)')
ylabel('RMSE(degree)')
legend('signal_1','signal_2');