%九單元線陣DOA估計 單信號加噪時Capon法%
 %信號處理部分                          
d=1; %天線陣元的間距
a=2;  %信號中心波長
b=20;%信號幅度
q1=30*(pi/180);%輸入信號的方向
A1=[exp(-2*pi*j*d*[0:8]*cos(q1)/a)]';%求陣因子
A=[A1];%得出A 矩陣
n=1:2000;
v1=.015;%信號的頻率

d=[b*cos(v1*n)];%輸入信號矢量
randn('state',1);
en=randn(9,2000);%高斯白噪聲
U=A*d;%總的輸入信號


%算法部分
U1=(U)';
c=cov(U*U1);%總輸入信號的協方差矩陣
ci=inv(c);%求協方差矩陣的逆矩陣
q1b=[pi/180:pi/180:pi];
for n2=1:length(q1b)
    q1a(n2)=q1b(n2);
    A1a=[exp(-2*pi*j*1*[0:8]*cos(q1a(n2))/a)]';
    Pcap(n2)=inv((A1a)'*ci*(A1a));%應用Capon 法估計輸出
    T(n2)=q1a(n2);
    P2=abs(Pcap);
end

%繪出應用 法估Capon 計的波達方向極坐標圖
subplot(1,2,1)
polar(T,P2)
xlabel('Capon 法估計的波達方向極坐標圖')
%繪出應用Capon 法估計的波達方向圖
subplot(1,2,2)
T1=T*180/pi;
semilogy(T1,P2);grid
xlabel('Capon 法估計的波達方向圖')