采用最大似然交替投影迭代的方法对信号进行DoA估计，很实用

采用最大似然交替投影迭代的方法对信号进行DoA估计，很实用的

    clear all
    clc
    
    S=5;     %信噪比
    M=7;     %阵元数
    N=3;     %信号数
    nd=100;    %对信号的采样点个数
    f=[1000000 1000000 1000000];        %信号频率矩阵
    t_max=5/1e4;       %对观测时间的限制
    ta=[0 40 60];%信号入射角度
    ta=ta*pi/180;
    dbc=0.5;       %阵元间距d与信号波长之间的比值
    ta_max=pi/2;   %最大张角
    s1=10^（S/20）;    %把信噪比由分贝的单位转化为对信号振幅的单位
    s2=1;
    j=sqrt（-1）;
    %phape=randn（Nnd）; %生成服从正态分布的随机矩阵（N*nd）阶  
    %phap=phape/max（max（phape））;
    t=linspace（0t_maxnd）;           %等差元素向量
    s=exp（j*2*pi*kron（f‘t））;%+j*phap）;  %输入信号（N*nd）阶
    s0=s1*s;                          %（N*nd）阶
    n=s2*randn（Mnd）;   %噪声模型 生成服从正态分布的随机矩阵（M*nd）阶
    a1=sin（ta）;          
    a2=（0:（M-1））‘;      
    a3=kron（a1a2）;          %信号到达各阵列的相位差
    a=exp（-1*j*2*pi*dbc*a3）; %导向矢量阵
    X=a*s0+n;                %信号模型（M*nd）阶
    R=X*X‘/nd;
%     [USV]=svd（R）;     %奇异值分解，计算主特征值并保留秩，并保留相应的特征向量矩阵
% disp（S）;
% Un=U（:N+1:M）;
% Gn=Un*Un‘;
% seaching_doa=-90:0.1:90 ;%现在的搜寻范围为-90到90
% for i=1:length（seaching_doa）
%     a_theta=exp（-j*（0:M-1）‘*pi*sin（pi*seaching_doa（i）/180））;   %角度
%     Pmusic（i）=1./abs（（a_theta）‘*Gn*（a_theta））; %功率谱函数
% end
% plot（seaching_doa10*log（Pmusic）‘r‘）
% xlabel（‘入射角/degree‘）;
% ylabel（‘空间谱/dB‘）;
% legend（‘Music Spectrum‘）;
% title（‘经典MUSIC估计‘）;
% grid on;
     %确定初始值 
         th1=（-90:0.1:90）‘;
         th=th1/180*pi;
         q=zeros（1N）;
         Q=zeros（1length（th））;
for          i=1:N            % 其中N表示信号个数
             p1=sin（q）;   
             p2=（0:（M-1））‘;
             p3=kron（p1p2）;
             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（Mi）;
            for     m=1:（i-1）             %表示其前i-1个数
                    pmu_b（:m）=p4（:m）;
            end
                    pmu_b（:i）=pmu_a（:1）;   %表示a（i）
                    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（1N）;
for p=1:N
     for     i=1:N
             p1=sin（q）;   
             p2=（0:（M-1））‘;
             p3=kron（p1p2）;
             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*p

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3668  2011-09-01 22:49  AP_ML.m