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]);