% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% 程序功能:四個輸入信號,八單元均勻圓陣的波達(dá)方向估計和波束形成    %
% 程序員:王浩 harvey                                              %
% 日期:20071018                                                   %
% 備注:應(yīng)用MUSIC和Capon方法                                       %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
clear;
q4=pi*5.8/6;q3=pi*3/2;q2=pi*1/3;q1=pi*2/3;%四個輸入信號的方向
E=1;
lam=1; %信號波長
r=1; %天線陣的半徑
m1=1;m2=1;m3=1;m4=1;m5=1;m6=1;m7=1;m8=1;%信號幅度
t=q1;
ps1=(pi*r*cos(t)/lam); % 求各個陣元相對圓心的相位差
ps2=(pi*cos(t+(pi/4))/lam);
ps3=(pi*cos(t+(pi/2))/lam);
ps4=(pi*r*cos(t+(.75*pi))/lam);
ps5=(pi*r*cos(t+(pi))/lam);
ps6=(pi*r*cos(1.25*pi+t)/lam);
ps7=(pi*r*cos((1.5)*pi+t)/lam);
ps8=(pi*cos((pi/4)-t)/lam);
aa1=m1*(exp(-j*(ps1)));
aa2=m2*(exp(-j*(ps2)));
aa3=m3*(exp(-j*(ps3)));
aa4=m4*(exp(-j*(ps4)));
aa5=m5*(exp(-j*(ps5)));
aa6=m6*(exp(-j*(ps6)));
aa7=m7*(exp(-j*(ps7)));
aa8=m8*(exp(-j*(ps8)));
A1=[aa1;aa2;aa3;aa4;aa5;aa6;aa7;aa8]; %求陣因子
t=q2;
ps1=(pi*r*cos(t)/lam);
ps2=(pi*cos(t+(pi/4))/lam);
ps3=(pi*cos(t+(pi/2))/lam);
ps4=(pi*r*cos(t+(.75*pi))/lam);
ps5=(pi*r*cos(t+(pi))/lam);
ps6=(pi*r*cos(1.25*pi+t)/lam);
ps7=(pi*r*cos((1.5)*pi+t)/lam);
ps8=(pi*cos((pi/4)-t)/lam);
ab1=m1*(exp(-j*(ps1)));
ab2=m2*(exp(-j*(ps2)));
ab3=m3*(exp(-j*(ps3)));
ab4=m4*(exp(-j*(ps4)));
ab5=m5*(exp(-j*(ps5)));
ab6=m6*(exp(-j*(ps6)));
ab7=m7*(exp(-j*(ps7)));
ab8=m8*(exp(-j*(ps8)));
A2=[ab1;ab2;ab3;ab4;ab5;ab6;ab7;ab8];
t=q3;
ps1=(pi*r*cos(t)/lam);
ps2=(pi*cos(t+(pi/4))/lam);
ps3=(pi*cos(t+(pi/2))/lam);
ps4=(pi*r*cos(t+(.75*pi))/lam);
ps5=(pi*r*cos(t+(pi))/lam);
ps6=(pi*r*cos(1.25*pi+t)/lam);
ps7=(pi*r*cos((1.5)*pi+t)/lam);
ps8=(pi*cos((pi/4)-t)/lam);
ac1=m1*(exp(-j*(ps1)));
ac2=m2*(exp(-j*(ps2)));
ac3=m3*(exp(-j*(ps3)));
ac4=m4*(exp(-j*(ps4)));
ac5=m5*(exp(-j*(ps5)));
ac6=m6*(exp(-j*(ps6)));
ac7=m7*(exp(-j*(ps7)));
ac8=m8*(exp(-j*(ps8)));
A3=[ac1;ac2;ac3;ac4;ac5;ac6;ac7;ac8];
t=q4;
ps1=(pi*r*cos(t)/lam);
ps2=(pi*cos(t+(pi/4))/lam);
ps3=(pi*cos(t+(pi/2))/lam);
ps4=(pi*r*cos(t+(.75*pi))/lam);
ps5=(pi*r*cos(t+(pi))/lam);
ps6=(pi*r*cos(1.25*pi+t)/lam);
ps7=(pi*r*cos((1.5)*pi+t)/lam);
ps8=(pi*cos((pi/4)-t)/lam);
ad1=m1*(exp(-j*(ps1)));
ad2=m2*(exp(-j*(ps2)));
ad3=m3*(exp(-j*(ps3)));
ad4=m4*(exp(-j*(ps4)));
ad5=m5*(exp(-j*(ps5)));
ad6=m6*(exp(-j*(ps6)));
ad7=m7*(exp(-j*(ps7)));
ad8=m8*(exp(-j*(ps8)));
A4=[ad1;ad2;ad3;ad4;ad5;ad6;ad7;ad8];
A=[A1 A2 A3 A4]; %得出A 矩陣
n=1:1900;
v1=.06;
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; %總的輸入信號
U1=(U)';
c=cov(U*U1); %總輸入信號的協(xié)方差矩陣
[s,z]=eig(c); % 求協(xié)方差矩陣的特征矢量及特征值
Vn=s(:,[5:8]); % 取出與零特征值對應(yīng)的特征矢量
ci=inv(c);
bb=[1 0 0 0]'; % 波束形成的條件:信號1增益為1,抑制信號2,3,4
Wopte=A'\bb; % 求解線性方程組（求Wopte）
q1b=[2*pi:-2*pi/180:2*pi/180];
for n=1:length(q1b)
h(n)=q1b(n);
ps1=(pi*r*cos(h(n))/lam);
ps2=(pi*cos(h(n)+(pi/4))/lam);
ps3=(pi*cos(h(n)+(pi/2))/lam);
ps4=(pi*r*cos(h(n)+(.75*pi))/lam);
ps5=(pi*r*cos(h(n)+(pi))/lam);
ps6=(pi*r*cos(1.25*pi+h(n))/lam);
ps7=(pi*r*cos((1.5)*pi+h(n))/lam);
ps8=(pi*cos((pi/4)-h(n))/lam);
ar1=m1*(exp(-j*(ps1)));
ar2=m2*(exp(-j*(ps2)));
ar3=m3*(exp(-j*(ps3)));
ar4=m4*(exp(-j*(ps4)));
ar5=m5*(exp(-j*(ps5)));
ar6=m6*(exp(-j*(ps6)));
ar7=m7*(exp(-j*(ps7)));
ar8=m8*(exp(-j*(ps8)));
A1a=[ar1;ar2;ar3;ar4;ar5;ar6;ar7;ar8];
Pmusic(n)=(A1a)'*A1a*(inv((A1a)'*Vn*(Vn)'*A1a));
Pcap(n)=inv((A1a)'*ci*(A1a));
Pcbf(n)=(A1a)'*c*(A1a);
T(n)=q1b(n);
Pcbf(n)=(A1a)'*c*(A1a);
P3=abs(Pcbf);
P1=abs(Pmusic); %Music 算法估計
P2=abs(Pcap); % Capon 算法估計
Ye(n)=Wopte'*A1a;
p=abs(Ye);
P3=abs(Pcbf);
end
figure(1)
T1=T*180/pi;
semilogy(T1,P1);grid % Music 算法波達(dá)方向估計
figure(2)
T1=T*180/pi;
semilogy(T1,P2);grid %Capon 算法波達(dá)方向估計
figure(3)
polar(T,p) %繪出應(yīng)用矩陣運(yùn)算求解加權(quán)系數(shù)后的波束形成的方向圖