Matlab波束形成程序.doc

word完整版)Matlab波束形成程序
波束形成与智能天线

1.均匀线阵方向图
%8阵元均匀线阵方向图，来波方向为0度
clc；
clear all；
close all;
imag=sqrt（-1）;
element_num=8;%阵元数为8
d_lamda=1/2；%阵元间距d与波长lamda的关系
theta=linspace(—pi/2,pi/2,200）;
theta0=0；%来波方向
w=exp（imag＊2*pi*d_lamda*sin（theta0)＊[0:element_num—1］'）；
for j=1:length(theta)
a=exp(imag*2*pi＊d_lamda*sin（theta(j))*[0:element_num-1］’）；
p(j)=w'＊a；
end
figure；
plot（theta,abs（p））,grid on
xlabel(’theta/radian’)
ylabel（'amplitude’)
title（’8阵元均匀线阵方向图’）

当来波方向为45度时，仿真图如下：

8阵元均匀线阵方向图如下,来波方向为0度，20log（dB）

随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下：

2。 波束宽度与波达方向及阵元数的关系
clc
clear all
close all
ima=sqrt(-1）；
element_num1=16； %阵元数
element_num2=128;
element_num3=1024；
lamda=0。03; ％波长为0.03米
d=1/2*lamda; ％阵元间距与波长的关系
theta=0:0。5：90；
for j=1：length（theta)；
fai(j）=theta（j）*pi/180—asin(sin（theta(j）＊pi/180）-lamda/(element_num1＊d））；
psi（j)=theta(j）＊pi/180-asin（sin(theta(j)*pi/180）-lamda/（element_num2＊d））;
beta（j）=theta（j）*pi/180—asin(sin(theta(j)*pi/180）-lamda/（element_num3＊d))；
end
figure；
plot（theta,fai，'r’，theta,psi'b’，theta,beta,'g’），grid on
xlabel（'theta'）;
ylabel（'Width in radians’)
title('波束宽度与波达方向及阵元数的关系'）
仿真图如下：

3. 当阵元间距 时,会出现栅瓣，导致空间模糊。仿真图如下：

4. 类似于时域滤波，天线方向图是最优权的傅立叶变换
仿真程序和仿真图如下：
clc
clear all
close all
ima=sqrt(-1）;
element_num=32; %阵元数
source_num=1； %信源数
d_lamda=1/2; %阵元间距与波长的关系
theta=linspace(-pi/2,pi/2，200)；
theta0=0; %来波方向
w=exp(ima*2*pi*d_lamda＊sin（theta0）＊［0:element_num—1］')；
for j=1：length（theta);
a=exp（ima*2*pi*d_lamda＊sin（theta（j）)＊［0:element_num-1］’);
p（j)=w'*a；
end
figure；
subplot(1,2，1）
plot(theta,abs（p)）,grid on
xlabel（'theta/radian’)
ylabel（'amplitude')
title（'按定义的方向图’)
pfft=fftshift(fft(w，128））；
subplot(1，2,2)
plot(linspace(-pi/2,pi/2，128），abs（pfft)）,grid on
xlabel（’theta/radian'）
ylabel（’FFT_amplitude’)
title(’最优权的傅里叶变换')

5。 ％最大信噪比准则方向图和功率谱
clc；
clear all;
close all;
ima=sqrt(—1）;
element_num=8; ％阵元数为8
d_lamda=1/2; ％间距为半波长
theta=-90：0。5:90; %范围
theta0=0; %来波方向
theta1=20; %干扰方向
L=512; %采样单元数
for i=1：L
amp0=10*randn（1);
amp1=200＊randn（1）；
ampn=1；
s(：，i）=amp0＊exp（ima*2*pi*1/2*sin(theta0＊pi/180）*［0:element_num—1]'）；
j(:，i)=amp1＊exp（ima＊2*pi*1/2＊sin（theta1＊pi/180）*［0：element_num-1］');
n(：，i)=ampn*(randn(element_num，1）+ima*randn（element_num，1））;
end
Rs=1/L*s＊s'； %信号自相关矩阵
Rnj=1/L＊（j＊j'+n＊n'）; ％干扰+噪声的自相关矩阵
［V，D]=eig(Rs，Rnj）; ％（Rs，Rnj）的广义特征值和特征向量
[D，I]=sort(diag(D））； ％排序
Wopt=V（:,I（8）)； ％最优权矢量
for j=1：length（theta)
a=exp(ima*2*pi*d_lamda＊sin(theta(j）＊pi/180)＊［0：element_num—1]’）;
f(j)=Wopt’＊a；
p（j)=a’*Rs＊a+a'＊Rnj*a；
end
F=20log10（abs(f)/max(max（abs(f）））);
P=20*log10(abs（p）/max(max(abs(p))）)；
subplot（121)
plot（theta,F)；grid on；hold on
plot（theta0，-50：0，'。’);plot(theta1,—50：0,'。'）
xlabel（'theta/0’)；ylabel(’F in dB'）;
title('max—SNR 方向图');
axis（[—90 90 -50 0]);
hold on
subplot(122）
plot(theta，P，'r’）；grid on
xlabel(’theta/0’)；ylabel('功率 in dB’)；
title('max—SNR功率谱’)
仿真图如下：

6. ％ASC旁瓣相消————MSE准则
clc;close all;clear all
ima=sqrt（—1）;
M=32; ％辅助天线的数目
d_lamda=.5;
theta0=-30; ％来波方向
theta1=60； ％干扰方向
L=512； ％采样单元数
s=zeros（1512)； ％预划分一个区域
for ii=1:L；
amp0=1＊randn(1); ％信号的幅度随机产生,保证信号之间是不相关的
amp1=200＊randn(1);
ampn=1；
jam(:,ii)=amp1*exp(ima*2*pi*0.5＊sin(theta1＊pi/180）*［0:M-1］')+ampn*(randn（M，1）+ima＊randn(M，1))； %干扰+噪声
s（ii）=amp0*exp(ima*2*pi＊0。5*sin(theta0*pi/180)）+amp1*exp(ima＊2＊pi*0。5＊sin(theta1*pi/180））+ampn＊（randn(1，1）+ima*randn(1，1）);％接收信号（信号+干扰+噪声）
s0（ii）=amp0*exp(ima*2*pi*0。5＊sin(theta0*pi/180)）；
end
Rx=1/L＊jam＊jam'； ％噪声自相关矩阵，相当于X（t）
r_xd=1/L＊jam＊s';
Wopt=pinv(Rx)*r_xd；
delta=s0-(s-Wopt’*jam）;
delta1=abs（mean(delta.^2）—（mean（delta)）.^2) ％方差
theta=linspace（-pi/2，pi/2，200)；
for jj=1：length（theta）
a=exp(ima*2*pi*。5＊sin(theta(jj）)*［0:M-1］');
f(jj)=Wopt'＊a;
end
F=20＊log10（abs（f)/（max(max（abs(f)）））)；
figure（1)
plot（theta＊180/pi,F),grid on，hold on
plot（theta0,-50：0,’。’）
plot(theta1，—50：0，’。’）
xlabel(’theta/o'）；
ylabel(’F/dB')；
title(’MSE准则下的方向图’)
axis（［-90 90 -50 0])；％可为x轴和y轴设置一个极限范围
仿真图如下:

7。 %线性约束最小方差(LCMV）准则
clc；
clear all ；
close all;
ima=sqrt(-1）；
element_num=8; %阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=—90:0。5:90—0。3; %搜索范围
theta0=0； 三个信号源的来波方向
theta1=30;
theta2=60;
L=512； ％采样单元数
for i=1：L；
amp0=10＊randn(1）;
amp1=100*randn(1）;
amp2=10*randn(1);
ampn=10；
x（：,i)=amp0*exp（ima＊2*pi＊1/2＊sin（theta0＊pi/180）＊［0：element_num-1］'）+.。。
amp1＊exp(ima*2＊pi＊1/2＊sin（theta1*pi/180）＊[0:element_num—1］')+.。.
amp2*exp(ima＊2*pi＊1/2＊sin（theta2＊pi/180）*[0：element_num-1］')+..。
ampn*（randn(element_num，1）+ima*randn(element_num,1）);
end
Rx=1/L*x＊ x';
steer1=exp（ima*2＊pi＊1/2*sin（theta0＊pi/180)＊［0:element_num-1]')
steer2=exp(ima＊2＊pi*1/2*sin（theta1*pi/180）＊[0:element_num—1］’)
steer3=exp(ima＊2*pi＊1/2*sin（theta2*pi/180)*［0:element_num-1］’)
C=[steer1 steer2 steer3］;
F=［1 0 1］’; ％把三个方向都作为来波方向
w=inv(Rx)*C＊(inv(C’*inv（Rx)＊C))＊F;
for j=1：length（theta）;
a=exp(ima*2＊pi*d_lamda*sin(theta(j）＊pi/180）＊[0：element_num—1］’）;
f（j)=w’*a;
p(j)=1/(a'*inv（Rx）*a);
end
f=10＊log10(abs(f)/（max（max(abs(f））)）);
figure(1)
subplot（121)
plot(theta,f）,grid on,hold on
plot(theta0，—20:0，’。'）
plot(theta1,—20：0,’.’）
plot(theta2，—20:0,’.’）
xlabel(’theta/o'）；
ylabel(’F/dB’)；
title（’Capon beamforming方向图'）
axis（［—90 90 -20 0］)；％可为x轴和y轴设置一个极限范围
P=10*log10(abs(p)/(max（max(abs(p)))）);
subplot(122）
plot（theta,P），grid on，hold on
plot（theta0，-20:0,'.’)
plot(theta1，—20:0,'。')
plot(theta2,—20:0,'。'）
xlabel('theta/o')；
ylabel（’功率/dB')；
title(’Capon beamforming功率谱')
仿真图如下：

8. %Capon beamforming
Clc;
clear all ；
close all；
ima=sqrt(—1）;
element_num=8; ％阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=—90：0.5:90； ％范围
theta0=0; ％来波方向
theta1=20; ％干扰方向
theta2=60； %干扰方向
L=1000； ％采样单元数
for i=1:L;
amp0=10*randn（1）;％信号的幅度随机产生,保证信号之间是不相关的
amp1=200＊randn(1);
amp2=200＊randn（1);
ampn=3；
x(:，i)=amp0*exp(ima*2*pi*1/2＊sin(theta0＊pi/180)＊[0:element_num—1]’)+。。.
amp1＊exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]’）+。。。
amp2＊exp(ima*2＊pi*1/2*sin(theta2*pi/180）＊[0:element_num-1］’）+.。。
ampn＊（randn（element_num，1)+ima*randn(element_num，1));
end
Rx=1/L＊x＊ x';
R=inv（Rx）；
steer=exp(ima＊2＊pi＊1/2＊sin（theta0＊pi/180）＊[0:element_num—1]'）;
w=R*steer/(steer'＊R*steer）;％Capon最优权矢量
for j=1:length(theta)；
a=exp(ima*2*pi＊d_lamda*sin（theta(j）＊pi/180)*[0：element_num—1］'）；
f(j)=w'＊a；
p（j)=1/（a'＊R＊a)；
end
F=20＊log10(abs(f)/(max(max(abs(f))）))；
P=20＊log10（abs（p）/（max（max（abs（p）))))；％此处是功率的对数形式
Figure;
subplot(121)
plot（theta,F)，grid on,hold on
plot（theta0，—50:0，'。'）
plot（theta1，—50:0,'。’）
plot(theta2,—50:0,'.’)
xlabel(’theta/o'）;
ylabel(’F/dB');
title（’Capon beamforming方向图'）
axis（［-90 90 -50 0]）;
subplot(122)
plot（theta，P）,grid on
xlabel（'theta/o’);
ylabel(’功率/dB');
title(’Capon beamforming功率谱’)
仿真图如下：

9。 ％不同方法估计协方差矩阵的Capon波束形成
clc;
clear all;
close all;
ima=sqrt（-1);
element_num=8； %阵元数为8
d_lamda=1/2; ％间距为半波长
theta=-90:0.5:90； ％范围
theta0=0; ％来波方向
theta1=50； %干扰方向
L=512;%采样单元数
for i=1:L
amp0=10＊randn（1)；
amp1=50＊randn(1）;
ampn=0。5;
s(:，i）=amp0*exp(ima*2*pi*1/2＊sin(theta0＊pi/180)*［0:element_num-1］’）；
j（:,i）=amp1＊exp(ima＊2＊pi*1/2＊sin(theta1＊pi/180)＊［0:element_num-1］')；
n(:，i）=ampn*exp（ima*2＊pi*randn(1)*［0：element_num—1]’）;
end
Rx=1/L*(s+j+n)＊(s+j+n)';%接收信号自相关矩阵
Rnj=1/L*(j+n）*（j+n)’；％干拢+噪声的自相关矩阵
e=exp（ima＊2*pi＊d_lamda*sin（theta0*pi/180)＊［0：element_num-1］'）;％来波方向信号
Wopt_Rx=inv(Rx）*e/(e'*inv(Rx)＊e);%采用接收信号的权矢量
Wopt_Rnj=inv（Rnj)*e/(e’*inv（Rnj)*e）；％采用干拢+噪声信号的权矢量
for j=1:length(theta)
a=exp(ima*2*pi*d_lamda*sin(theta（j）＊pi/180)＊[0:element_num-1]'）;
f1(j）=Wopt_Rx’*a；
f2(j)=Wopt_Rnj’＊a；
end
F1=20*log10(abs（f1）/max（max(abs(f1)））)；
F2=20*log10（abs(f2）/max（max(abs（f2）)))；
figure;
plot(theta，F1,theta，F2,’r’）,grid on;hold on
plot(theta0,—50:0，'.’);plot(theta1,—50:0，’。')
xlabel('theta/0')；ylabel(’F（1,2）/dB’）；
title(’不同方法估计协方差矩阵的Capon波束形成’）；
axis(［-90 90 —50 0］)；
仿真图如下：

10．%多点约束的Capon波束形成和方向图
clc；
clear all ;
close all;
ima=sqrt（—1）；
element_num=8； %阵元数
d_lamda=1/2; %阵元间距与波长的关系
theta=-90:0.3:90； %搜索范围
theta0=0； %来波方向
theta1=20; ％干扰方向
theta2=50; ％干扰方向
L=512； %采样单元数
Rx=zeros（element_num，element_num）；％产生协方差矩阵
for i=1:L；
amp0=10*randn(1)；
amp1=10*randn(1);
amp2=50＊randn(1）；
namp=0.5＊randn（1);%噪声的幅度随机产生，保证噪声与信号之间是不相关的
J（:，i)=amp1＊exp（ima＊2*pi＊1/2＊sin（theta1＊pi/180）*[0：element_num-1］’)+.。。
amp2*exp(ima＊2＊pi＊1/2*sin（theta2*pi/180)*［0:element_num-1]’)+..。
namp＊exp（ima＊2*pi＊randn（1）＊［0:element_num—1］’);
x（:,i）=amp0＊exp（ima＊2*pi＊1/2*sin(theta0*pi/180)＊[0：element_num-1]'）+.。。
J(：，i）;％表示接收信号
end
Rx=Rx+1/L＊x*x’；
R=inv（Rx);
w=amp0*exp(ima＊2*pi*1/2＊sin(theta0*pi/180）＊[0：element_num—1］'）+。。。
amp2*exp（ima＊2＊pi*1/2*sin（theta2＊pi/180）*[0:element_num—1]’）+...
amp1＊exp（ima*2*pi*1/2＊sin(theta1＊pi/180）*[0：element_num-1］'）;
for j=1：length（theta);
a=exp(ima＊2＊pi*d_lamda*sin(theta（j）*pi/180)*［0:element_num—1]’)；
f（j)=w'＊a；
p（j）=1/（a’＊R*a)；
end
F=10*log10（abs(f）/max(max(f)));％取对数的方向图
P=10＊log10(abs（p)/max(max（p）））；%取对数的功率谱估计
figure;
subplot(121)
plot（theta，F)，grid on，hold on
plot(theta0，-50：0，’。’)；plot(theta1,—50：0，'。’）;plot（theta2,-50：0，'。'）;
xlabel('theta/o’）；
ylabel('F/dB’)；
title(’capon beamforming方向图’)
subplot(122)
plot（theta,P）,grid on,hold on
plot（theta0，—50：0,'。');plot（theta1，-50：0，’.'）；plot(theta2,-50:0，'。’);
xlabel('theta/o'）;
ylabel（'P/dB’）；
title(’capon beamforming功率谱'）
仿真图如下：

大作业
%自适应波束形成方向图
clc；
clear all;
close all;
ima=sqrt（-1）；
element_num=8； %阵元数
c=3e8;
f=500e6；
lamda=c/f；
d=1/2*lamda; ％阵元间距与波长的关系
theta=—90:0.5:90； ％范围
theta0=0; ％来波方向
theta1=45； ％干扰方向
theta2=60; %干扰方向
L=512; ％采样单元数
for i=1:L；
amp0=10＊randn（1）；%信号的幅度随机产生，保证信号之间是不相关的
amp1=100*randn（1);
amp2=100＊randn(1）；
ampn=10;
x(:,i)=amp0＊exp（ima＊2＊pi*1/2＊sin（theta0*pi/180)*［0:element_num-1］’）+.。。
amp1*exp（ima＊2＊pi＊1/2＊sin(theta1＊pi/180)*[0:element_num—1］’）+。.。
amp2*exp（ima*2＊pi*1/2*sin（theta2*pi/180）＊［0:element_num-1］’)+.。.
ampn*(randn(element_num，1）+ima＊randn(element_num,1））；％产生阵列信号
end
Rx=1/L＊x＊ x’;
steer1=exp（ima*2＊pi＊1/2＊sin（theta0＊pi/180）＊［0:element_num—1］'）
steer2=exp（ima*2＊pi＊1/2＊sin（theta1*pi/180）*［0：element_num-1]')
steer3=exp(ima＊2*pi*1/2＊sin（theta2＊pi/180)*［0：element_num-1]'）
C=［steer1 steer2 steer3];
F=［1 0 0］';%把0度方向作为来波方向
w=inv(Rx)＊C*（inv（C'*inv(Rx)＊C）)*F;
for j=1:length(theta)；
a=exp（ima＊2＊pi＊1/2＊sin（theta（j）*pi/180)＊［0:element_num-1]’)；
f（j)=w'＊a；
p（j)=1/（a'＊inv（Rx)*a）;
end
F=20*log10(abs（f)/(max(max（abs（f)））));
P=20*log10（abs（p）/(max(max（abs(p)）)));％此处是功率的对数形式
figure
subplot（121）
plot（theta，P）,grid on,hold on
xlabel（'theta/o'）;
ylabel（'F/dB’);
title('阵列信号方向图’）
axis（[-90 90 -50 0］)；
subplot（122）
plot（theta，F)，grid on,hold on;
plot(theta0,—50:0,'。’)
plot(theta1，-50:0,'。')
plot（theta2,—50:0，’。’）
xlabel('theta/o'）；
ylabel('功率/dB')；axis(［—90 90 —50 0])；
title（'自适应波束形成方向图')
仿真图如下：
：

当采样数L=2048时仿真图如下：

当采样数L=2048时仿真图如下：

当L=512，幅度为
amp0=10＊randn（1）;
amp1=200*randn（1)；
amp2=200＊randn(1）；
ampn=10；时仿真图如下：

当L=512,幅度为
amp0=100*randn(1)；
amp1=200*randn(1）；
amp2=200*randn(1);
ampn=10；时仿真图如下：

当L=512，幅度为
amp0=100＊randn（1)；
amp1=10＊randn（1）;
amp2=10＊randn(1）；
ampn=10;时仿真图如下：

当L=2048，幅度为
amp0=10＊randn(1）；
amp1=50＊randn（1）；
amp2=50＊randn(1)；
ampn=10;时仿真图如下:

当L=2048,幅度
amp0=10＊randn(1）；
amp1=500＊randn（1）;
amp2=500＊randn（1)；
ampn=10；时仿真图如下：