% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 程序功能:四个输入信号,八单元均匀圆阵的波达方向估计和波束形成 %
% 程序员:王浩 harvey %
% 日期:20071018 %
% 备注:应用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); %总输入信号的协方差矩阵
[s,z]=eig(c); % 求协方差矩阵的特征矢量及特征值
Vn=s(:,[5:8]); % 取出与零特征值对应的特征矢量
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 算法波达方向估计
figure(2)
T1=T*180/pi;
semilogy(T1,P2);grid %Capon 算法波达方向估计
figure(3)
polar(T,p) %绘出应用矩阵运算求解加权系数后的波束形成的方向图
评论0