ap2.m

相干源对最大似然算法测向的影响 设置了两个相干的信号

    clear
    S=5;     %信噪比
    M=7;     %阵元数
    N=2;     %信号数
    nd=100;    %对信号的采样点个数
    f=[1000000 1000000];        %信号频率矩阵
    t_max=5/1e4;       %对观测时间的限制
    ta=[15 30];%信号入射角度
    ta=ta*pi/180;
    dbc=0.5;       %阵元间距d与信号波长之间的比值
    ta_max=pi/2;   %最大张角
    s1=10^(S/20);    %把信噪比由分贝的单位转化为对信号振幅的单位
    s2=1;
    %phape=randn(N,nd); %生成服从正态分布的随机矩阵（N*nd）阶  
    %phap=phape/max(max(phape));
    t=linspace(0,t_max,nd);           %等差元素向量
    s=exp(j*2*pi*kron(f',t));%+j*phap);  %输入信号（N*nd）阶
    s0=s1*s;                          %（N*nd）阶
    n=s2*randn(M,nd);   %噪声模型 生成服从正态分布的随机矩阵（M*nd）阶
    a1=sin(ta);          
    a2=(0:(M-1))';      
    a3=kron(a1,a2);          %信号到达各阵列的相位差
    a=exp(-1*j*2*pi*dbc*a3); %导向矢量阵
    X=a*s0+n;                %信号模型（M*nd）阶
    R=X*X'/nd;          
     %确定初始值 
         th1=[-90:0.1:90]';
         th=th1/180*pi;
         q=zeros(1,N);
         Q=zeros(1,length(th));
for          i=1:N
             p1=sin(q);   
             p2=(0:(M-1))';
             p3=kron(p1,p2);
             p4=exp(-1*j*2*pi*dbc*p3);%第K(K>1)个信号的导向矢量
        for  k=1:length(th)
             pmu_a1=sin(th(k));
             pmu_a2=(0:(M-1))';
             pmu_a3=pmu_a1*pmu_a2; 
             pmu_a=exp(-1*j*2*pi*dbc*pmu_a3);
         if         i==1
                    P=pmu_a*pinv(pmu_a'*pmu_a)*pmu_a';
                    qk=trace(P*R);
                    Q(k)=qk;
         else   
                    pmu_b=zeros(M,i);
            for     m=1:(i-1)
                    pmu_b(:,m)=p4(:,m);
            end
                    pmu_b(:,i)=pmu_a(:,1);
                    P=pmu_b*pinv(pmu_b'*pmu_b)*pmu_b';
                    qk=trace(P*R);
                    Q(k)=qk;
         end   
        end
                    b=max(Q);
         for        k=1:length(th)      
                if (b-Q(k))==0
                     c=k;
                    q(i)=th(c);
                end
         end
end

     %迭代
             d=zeros(1,N);
for          p=1:N
     for     i=1:N
             p1=sin(q);   
             p2=(0:(M-1))';
             p3=kron(p1,p2);
             p4=exp(-1*j*2*pi*dbc*p3);%第K(K>1)个信号的导向矢量
        for  k=1:length(th)
             pmu_a1=sin(th(k));
             pmu_a2=(0:(M-1))';
             pmu_a3=pmu_a1*pmu_a2; 
             pmu_a=exp(-1*j*2*pi*dbc*pmu_a3);
             p4(:,i)=pmu_a(:,1);
             P=p4*pinv(p4'*p4)*p4';
             qk=trace(P*R);
             Q(k)=qk;
       end   
                      b=max(Q);
                      g(i)=abs(b);
           for        k=1:length(th)    
                if (b-Q(k))==0
                      c=k;
                      d(i)=th(c);
                end
           end      
     end
                q=d;
end
              q=sort((q)*180/pi)
              figure(1);
              plot(th1,abs(Q));
              grid on;
              hold on;
              plot(q,g,'*');
              hold off;