music_line.m

实现基本MUSIC算法,加入高斯白噪声,完成信号来波方向的精确定位

        clear all;
        clc;
        d = 0.11;
        v = 0.35;
        b = 0.75;
        N = 500;                            %设定采样点数
        SNR = [0 0];                         
        fs=8;                             %设定采样频率
                                 %阵元间距与波长的比值
        f=[0.9 0.6];                            %设两个信号元的频率，单位为MHz
        theta = [-30 60]*pi/180;              %入射方向角向量
        phi =[60 20]*pi*180;                 %两入射信号初始相位
        M = 4;                             %设定阵元数
        AMP = 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*d*m'*(sin(theta)*b./v));                             %方向矩阵
        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*d*b/v*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');
       xlabel('DOA(degree)')
       ylabel('P(dB)')
       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 = .1;
          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