%七單元天線陣DOA估計
clear;clc;
d=1;               %天線陣元的間距
lma=2;             %信號中心波長
q1=1*pi/4;q2=1*pi/3;q3=1*pi/6;q4=3*pi/4;%四輸入信號的方向
A1=[exp(-2*pi*j*d*[0:6]*cos(q1)/lma)]';%求陣因子
A2=[exp(-2*pi*j*d*[0:6]*cos(q2)/lma)]';
A3=[exp(-2*pi*j*d*[0:6]*cos(q3)/lma)]';
A4=[exp(-2*pi*j*d*[0:6]*cos(q4)/lma)]';
A=[A1,A2,A3,A4];                       %得出A矩陣
n=1:1900;
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);
    A1a=[exp(-2*pi*j*1*[0:6]*cos(q1a(n))/lma)]';
    Pcapon(n)=(A1a)'*A1a*(inv((A1a)'*Vn*(Vn)'*A1a));%應用Cpon法估計輸出
    T(n)=q1a(n);
    P=abs(Pcapon);
end
figure(1)                                          %繪出各波達方向圖
   polar(T,P);title('Capon-1')
  figure(2);
  T1=T*180/pi;
  semilogy(T1,P);grid;title('Capon-2')