ml.m

此算法为基于最大似然方法的DOA估计 是区别于子空间类的DOA算法 但需要多维搜索 计算量较大

clear
close all
Fs=3e6;%采样频率
N=100;%采样点数
n=0:(N-1);%采样点
n=n/Fs;
f0=1e6;%中心频率
m=16;%阵元数
p=3;
bi=0.5;
j=sqrt(-1);
i=j;
degrad=pi/180;  %将角度制化为弧度制  
angle=[0;4;60];%波达方向
x1=exp(i*(2*pi*n*f0+pi/6));%信号1
x2=exp(i*(2*pi*n*f0/2+pi/3));%信号2
x3=exp(i*(2*pi*n*f0/20+pi/3));%信号3
SNR=[5;5;5];%信噪比
sign=[x1;x2;x1];
nr=randn(m,N);  %建立高斯随机序列
ni=randn(m,N);
E=nr+ni*j;        %建立高斯噪声
Ps=sign*(sign/N)';    %方向向量
Ps1=diag(Ps);
refp=2*10.^(SNR./10);     %2代表的是zao sheng的功率
tmp=sqrt(refp./Ps1);
S1=diag(tmp)*sign;     %求出在给定信噪比下信号的信号
tm=bi*2*i*pi*sin(angle'.*degrad);  %相位差   阵元间距为波长的一半
tmp1=[0:m-1]';        %信号到达阵元的相位差组成的向量
A=tmp1*tm;
A1=exp(A);            %m*p维方向矩阵
X=(A1*S1)+E;           %第k个阵元上的接收信号
Rxx=X*(X/N)';       %接收信号的协方差
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
th2=[-90:0.1:90]'; %线性信号角度范围
%%%%%%%%%%%%%%第一次%%%%%%%%%%%%%%%%

    for k=1:length(th2);
    th=[-2,7,th2(k)];
    tmp2=i*pi*sin(th.*degrad); 
    p=tmp1*tmp2; 
    p=exp(p);
    Pa=p*inv(p'*p)*p';
    ang(k)=sum(diag(Pa*Rxx));
    end;
    [ttt,k]=max(ang);
    k1=th2(k);
    figure(1);
    plot(th2,ang);
    grid on;
    %%%%%%%%%%%%%第二次%%%%%%%%%%%%%
    for k=1:length(th2);
    th=[k1,7,th2(k)];
    tmp2=i*pi*sin(th.*degrad); 
    p=tmp1*tmp2; 
    p=exp(p);
    Pa=p*inv(p'*p)*p';
    ang(k)=sum(diag(Pa*Rxx));
    end;
    [ttt,k]=max(ang);
    k2=th2(k);
    figure(2);
    plot(th2,ang);
    grid on;
    %%%%%%%%%%%第三次%%%%%%%%%%%%%%%%%
    for k=1:length(th2);
    th=[k1,k2,th2(k)];
    tmp2=i*pi*sin(th.*degrad); 
    p=tmp1*tmp2; 
    p=exp(p);
    Pa=p*inv(p'*p)*p';
    ang(k)=sum(diag(Pa*Rxx));
end;
    [ttt,k]=max(ang);
    k3=th2(k);
    figure(3);
    plot(th2,ang);
    grid on;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    th=sort([k1,k2,k3]);