music.m

通过music算法实现对信号源的测向 DOA估计

clear;
   snr=3;    %信噪比
   se1=1;    %选择输入信号的波形
   ta=[5,15,35,80]*pi/180;  %入射信号角度
   ma=12;           %天线阵元个数
   dbc=0.5;   %阵元间距d与信号波长之间的比值
   F=1.8e9;   %信号频率
   f=[F,F,F,F];
   t_max=3.17/1e3;   %对观测时间的限制
   nd=100000;   %对信号的采样点个数
   md=1000;     %?
   istat=0;  %选择确定最小特征值个数的准则
   ta_max=pi/2;   %最大张角
   l1=max(size(ta));
   s1=10.^(snr/20);   %把信噪比由分贝的单位转化为对信号振幅的单位
   s2=1;
   
   %信号的选择
   phape=randn(l1,nd); %生成服从正态分布的随机矩阵（l1*nd）阶
   phap=phape/max(max(phape));              
   t=linspace(0,t_max,nd);      %等差元素向量
   if(se1==1)
     s=sin(2*pi*kron(f',t)+phap);  %选择正弦输入信号（l1*nd）阶
   elseif(se1==2)
     s=sawtooth(2*pi*(kron(f',t)+phap));  %选择锯齿波输入信号
   elseif(se1==3)
     s=phap;   %选择正态分布的随机信号
   end
   s0=s1*s;   %（l1*nd）阶
   
   %构造协方差矩阵
   n=s2*randn(ma,nd);   %噪声模型 生成服从正态分布的随机矩阵（ma*nd）阶
   a1=sin(ta);          %（1*4）阶
   a2=(0:(ma-1))';      %（12*1）阶
   a3=kron(a1,a2);      %信号到达各阵列的相位差（12*4）阶
   a=exp(-1j*2*pi*dbc*a3); %导向矢量阵（12*4）阶
   x=a*s0+n;            %信号模型（ma*nd）阶
   r=x*x'/nd;           %信号协方差矩阵（ma*ma）阶
   [z1,d]=eig(r);       %求协方差矩阵的特征值和特征向量，d为特征值矩阵
   [w,k]=sort(diag(d)); %特征值排序（升序），diag（d）生成对角阵
   w=fliplr(w');        %小特征值矩阵，噪声子空间？颠倒
   z=z1(:,k);                    %？？？？？？？？？？？？？
   %确定最小特征值个数的准则
   fns_w=fliplr(w);           % 为什么不直接用 fns_w=w'？
   fns_s=cumsum(fns_w);       %?每列为累积和
   fns_a=fliplr(fns_s);
   fns_p=cumprod(fns_w);      %每列为乘积
   fns_b=fliplr(fns_p); 
   alr=zeros(ma,1);
   mdl=zeros(ma,1);
   aic=zeros(ma,1);
   for m=0:(ma-1)    
       alr(m+1)=(ma-m)*nd*log(fns_a(m+1)/(fns_b(m+1).^(1/(ma-m)))/(ma-m));   % 似然函数？
       aic(m+1)=alr(m+1)+m*(2*ma-m+1);
       mdl(m+1)=alr(m+1)+m*(2*ma-m+1)*log(nd)/2;
   end
   if istat==1
      fnsmin=min(aic);
      fnscom=(aic==fnsmin);      %选择AIC准则
   else
      fnsmin=min(mdl);
      fnscom=(mdl==fnsmin);      %选择MDL准则          ？？？？？？？
   end
   nx=ma-(min(find(fnscom))-1);  %返回矩阵非零元素的位置     min可以去掉？
   
   %构造噪声矩阵
   pmu_a1=sin(linspace(0,ta_max,md));
   pmu_a2=(0:(ma-1))';
   pmu_a3=kron(pmu_a1,pmu_a2);
   pmu_a=exp(-1j*2*pi*dbc*pmu_a3);
   pmu_en=z(:,1:nx);

   %求空间谱函数
   for i=1:md
      pmu_d(i)=(pmu_a(:,i))'*pmu_en*(pmu_en)'*pmu_a(:,i);
   end
   pmu_v=1./pmu_d;
   pmu=abs(pmu_v);
   
   %空间谱函数作图，并标出谱峰
   j=0;
   for i=2:md-1
      if (pmu(i-1)<=pmu(i))
         if (pmu(i)>=pmu(i+1))
            j=j+1;
            dd(j)=pmu(i);
            df(j)=i;
         end
      end
   end
   pmu_n=j;
   df=df*ta_max/pi*180/md;
   plot(linspace(0,ta_max/pi*180,md),pmu);
   hold on;
   plot(df,dd,'*');
   hold off;
   pmu_n
   dd
   df