clear all;
clc;
N=5;%表示陣元數(shù)
q=3;%表示信號數(shù)目
ADOA=[50,100,200,320,250]*pi/180;%信號方位角
FDOA=[90,90,90,90,90]*pi/180; %俯仰角；
f=[5e6,4.96e6,5.02e6,5.01e6,4.95e6]; %信號頻率Hz
fc=5e6;
c=3e8;  %光速
r=1.25*c/fc;%圓陣半徑
ep=2*pi*r*fc/c;

sn=512;%采樣數(shù)目
fs=4*max(f);%采樣頻率
snr=10; %信噪比
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%a陣列響應(yīng)向量
A=zeros(N,q); 
for kk=1:q
    yebu=2*pi*f(kk)/c*r*sin(FDOA(kk));
for k=1:N
    A(k,kk)=exp(j*yebu*cos(ADOA(kk)-(k-1)*2*pi/N));
end
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tau=1/fs;
t=tau:tau:sn*tau;
s1=exp(j*2*pi*f(1)*t);%信號源
s2=5*exp(j*2*pi*f(2)*t);%信號源
s3=exp(j*2*pi*f(3)*t);%信號源
s4=exp(j*2*pi*f(4)*t);%信號源
s5=exp(j*2*pi*f(5)*t);%信號源
ss=[s1;s2;s3;s4;s5];%信號源
s=ss(1:q,:);
%cor=abs(corrcoef(s'));
%%%%%%%%%%%%%%%%%%%%%%%
y2=A*s(1:q,:);
y1=awgn(y2,snr,0);
number=q;
R=y1*y1'/sn;%采樣方差矩陣
[v,dd]=eig(R);
if dd(1,1)>dd(2,2)
    En=v(:,number+1:N);
else
    En=v(:,1:N-number);
end
%En代表噪聲子空間
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MUSIC譜
kh=1;
for kf=0:360    
    doa1=kf*pi/180;
    for k3=1:N
        aa(k3,1)=exp(j*ep*sin(pi/2)*cos(doa1-(k3-1)*2*pi/N));
    end 
    pu(kh)=10*log10(abs(aa'*aa/(aa'*En*En'*aa)));
    kh=kh+1;
end
pu(pu<=-3)=-3;

 figure(1);

plot(0:360,pu,'b','linewidth',1.5);
 hold on;
 xlabel('方位角');
 ylabel('C-MUSIC譜');
 for k=1:q
 [k1,k2]=max(pu);
 DOA1(k)=(k2-1);
 pu(k2)=0;
 end
 DOA1