clear all
d=1;                                 %天線陣元間距
lma=2;                               %信號中心波長
q1=1*pi/4;
q2=pi/3;q3=pi/6;q4=3*pi/4;           %四個輸入信號的方向
A1=[exp(-2*pi*i*d*[0:6]*cos(q1)/lma)]';      %求陣因子
A2=[exp(-2*pi*i*d*[0:6]*cos(q2)/lma)]';
A3=[exp(-2*pi*i*d*[0:6]*cos(q3)/lma)]';
A4=[exp(-2*pi*i*d*[0:6]*cos(q4)/lma)]';
A=[A1,A2,A3,A4];                     %得出A矩陣
n=1:2000;
v1=.015;                             %四個信號的頻率
v2=.05;
v3=.02;
v4=.035;
d=[1.3*cos(v1*n);1*sin(v2*n);1*sin(v3*n);1*sin(v4*n)];  %輸入信號矢量
U=A*d;                               %總的輸入信號
U1=(U)';
c=cov(U*U1);                         %總輸入信號的協方差矩陣
[s,h]=eig(c);                        %求協方差矩陣的特征矢量及特征值
Vn=s(:,[5:7]);                       %取出與零特征值對應的特征矢量
ci=inv(c);                           %求協方差矩陣的逆矩陣
q1b=[pi/180:pi/180:pi];
for n=1:length(q1b)
    q1a(n)=q1b(n);
    Ala=[exp(-2*pi*i*1*[0:6]*cos(q1a(n))/lma)]';
Pmusic(n)=(Ala)'*Ala*((inv(Ala)'*Vn*(Vn)'*Ala));  %應用MUSIC法估計輸出
Pcap(n)=inv((Ala)'*ci*(Ala));
T(n)=qla(n);
P1=abs(Pmusic);
P2=Abs(Pcap);
end
figure(1);                            %繪出應用MUSIC法估計的波達方向圖
polar(T,P1);
figure(2)                           %繪出應用MUSIC法估計的波達方向圖
T1=T*180/pi
semilogy(T1,P1);grid
figure(3)                            %繪出應用Capon法估計的波達方向圖
polar(T,P2)
figure(4)                            %繪出應用Capon法估計的波達方向圖
T1=T*180/pi;
semilogy(T1,P2);grid