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