clear;
q1=pi*2/3;q2=pi*1/3;q3=pi*3/2;q4=pi*5.8/6;                  %四個輸入信號的方向
E=1;                                                        
lam=1;                                                      %信號波長
r=1;                                                        %天線陣的半徑
m1=1;m2=1;m3=1;m4=1;m5=1;m6=1;m7=1;m8=1;                    %信號幅度
a=0:7;%a的取值
for n=1:length(a)
    A1(n)=[exp(-j*pi*cos(q1-(pi*a(n)/4))/lam)];%第n個陣元對第1個信號的影響
    A2(n)=[exp(-j*pi*cos(q2-(pi*a(n)/4))/lam)];%第n個陣元對第2個信號的影響
    A3(n)=[exp(-j*pi*cos(q3-(pi*a(n)/4))/lam)];%第n個陣元對第3個信號的影響
    A4(n)=[exp(-j*pi*cos(q4-(pi*a(n)/4))/lam)];%第n個陣元對第4個信號的影響
end
A=[A1',A2',A3',A4'];%得出A矩陣——8x4                                         
n=1:1900;%n的取值（抽樣）
v1=.06;%v1-v4代表四個信號的頻率
v2=.02;
v3=.03;
v4=.073;
D=[1*cos(v1*n);1*sin(v2*n);1*sin(v3*n);1*square(v4*n)];      %四個輸入信號
U=A*D;   %8x1900                                             %總的輸入信號
U1=(U)'; %1900x8
c=cov(U*U1);   %8階方陣                                      %總輸入信號的協(xié)方差矩陣
[s,z]=eig(c)                                                %求協(xié)方差矩陣的特征矢量及特征值
Vn=s(:,[1:4]);                                               %取出與零特征值對應(yīng)的特征矢量
ci=inv(c); %求逆矩陣，得到8階方陣
%bb=[1 0 0 0]';
%Wopte=A'\bb;                                                 %求解線性方程組（求wopte)
q1b=[2*pi:-2*pi/180:2*pi/180];  %0-2pi內(nèi)均勻分割——1x180
b=0:7; %b的取值——1x8
for t=1:length(q1b)%開始第一重循環(huán)
    h(t)=q1b(t);%0-2*pi內(nèi)均勻取值
    for m=1:length(b)
        Aa(m)=[exp(-j*pi*cos(h(t)-(pi*a(m)/4))/lam)];%1x8矩陣——陣列流型
    end
    A1a=Aa';%求轉(zhuǎn)置矩陣
    Pmusic(t)=(A1a)'*A1a*(inv((A1a)'*Vn*(Vn)'*A1a));
    Pcap(t)=inv((A1a)'*ci*(A1a));
    T(t)=q1b(t); %將q1b取值反應(yīng)到T坐標(biāo)軸上                     
    P1=abs(Pmusic);                                           %MUSIC算法估計(jì) 
    P2=abs(Pcap);                                             %CAPON算法估計(jì)
    %Ye(t)=Wopte'*A1a;
    %p=abs(Ye);
end
figure name 'MUSIC'
T1=T*180/pi;
semilogy(T1,P1);grid                                          %MUSIC算法波達(dá)方向估計(jì)
figure name 'Polar plans-MUSIC'%MUSIC極坐標(biāo)圖
polar(T1,P1)
figure name 'CAPON'%CAPON極坐標(biāo)圖
T1=T*180/pi;
semilogy(T1,P2);grid                                          %CAPON算法波達(dá)方向估計(jì)
figure name 'Polar plans-CAPON'
polar(T1,P2)
%figure(3)
%polar(T,p)                                                    %繪出應(yīng)用矩陣運(yùn)算求解加權(quán)系數(shù)后的波束