波束形成_Matlab程序

1、1均匀线阵方向图%8阵元均匀线阵方向图，来波方向为0度clc;clear all;close all;imag=sqrt(_1);element_num=8;% 阵元数为 8d_lamda=1/2;%阵元间距d与波长lamda的关系theta=li nspace(-pi/2,pi/2,200);theta0=45/180*pi;%来波方向(我觉得应该是天线阵的指向)%theta0=0;%来波方向w=exp(imag*2*pi*dl_lamda*si n(theta0)*0:eleme nt_nu m-1');for j=1:le ngth(theta)%(我认为是入射角度，即来波方向，

2、计算阵列流形矩阵A)a=exp(imag*2*pi*damda*si n(theta(j)*0:eleme nt_nu m-1');p(j)=w'*a; %(matlab中的默认为共轭转置，如果要计算转置为w.'*a)endfigure;plot(theta,abs(p),grid onxlabel('theta/radia n')ylabel('amplitude')title('8阵元均匀线阵方向图')见张小飞的书阵列信号处理的理论和应用 2.3.4节阵列的方向图8阵元均勻线阵方向图T TTT-1-0.500.511.

3、52th etaradian当来波方向为45度时，仿真图如下:8阵元均匀线阵方向图来液方向为4E度8阵元均匀线阵万向图如下，来波万向为0度，20log (dB)吕阵元均勾銭阵方向图.来液方向为0産mp u- ap=d 常随着阵元数的增加，波束宽度变窄，分辨力提高：仿真图如下:11=32IAIiIIti4I1 5-1-0 500 511 5theta'radiano o o co- CN| II2. 波束宽度与波达万向及阵元数的关系clcclear allclose allima=sqrt(_1);eleme nt_nu m1=16;% 阵元数eleme nt_n um2=128;ele

4、me nt_nu m3=1024;lamda=0.03;%波长为 0.03 米d=1/2*lamda;%阵元间距与波长的关系theta=0:0.5:90;for j=1:le ngth(theta);fai(j)=theta(j)*pi/180-asi n(sin (theta(j)*pi/180)-lamda/(eleme nt_n um1*d); psi(j)=theta(j)*pi/180-asi n(s in (theta(j)*pi/180)-lamda/(eleme nt_n um2*d);beta(j)=theta(j)*pi/180-asi n(sin (theta(j)*pi

5、/180)-lamda/(eleme nt_num3*d); endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下：Q.70.6波東宽痩与波达方向及阵元数的黄乘5 4 3 2 10.0.0.0.0.i-I7060503et出O3020当阵元间距 d /2时，会出现栅瓣，导致空间模

6、糊。仿真图如下:2205060S0505 -1订-2-2胡-344 HP 匸一lamda为0一9倍波长0 1 theta/radianiheta/radian23. 类似于时域滤波，天线方向图是最优权的傅立叶变换 仿真程序和仿真图如下： 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; %来波方向 (ayy 应该是阵列指向方向 ) w=exp(ima*2*pi*d_l

7、amda*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;endfigure;subplot(1,2,1) plot(theta,abs(p),grid on xlabel('theta/radian') ylabel('amplitude') title(' 按定义的方向图 ') pfft=fftshift(fft(w,128);sub

8、plot(1,2,2) plot(linspace(-pi/2,pi/2,128),abs(pfft),grid on xlabel('theta/radian') ylabel('FFT_amplitude') title(' 最优权的傅里叶变换 ')o0 1llietaadianSOS3 3 22Q-8n=_d era按定义的方向图theta/radian4. %最大信噪比准则方向图和功率谱clc;clear all; close all; ima=sqrt(-1); eleme nt_num=8; d_l amda=1/2; theta=-

9、90:0.5:90; thetaO=O;theta1=20; L=512;for i=1:L%阵元数为8%间距为半波长%范围%来波方向%干扰方向%采样单元数Rs=1/L*s*s'Rnj=1/L*(j*j'+n*n'); V,D=eig(Rs,Rnj); D,l=sort(diag(D); Wopt=V(:,l(8);for j=1:le ngth(theta)amp0=10*ra ndn (1);amp1=200*ra ndn (1);amp n=1;s(:,i)=amp0*exp(ima*2*pi*1/2*si n(theta0*pi/180)*0:eleme nt_

10、num-1'); j(:,i)=amp1*exp(ima*2*pi*1/2*si n(theta1*pi/180)*0:eleme nt_num-1');n (:,i)=amp n*(ra ndn (eleme nt_nu m,1)+ima*ra ndn (eleme nt_nu m,1);end%信号自相关矩阵%干扰+噪声的自相关矩阵%（ Rs,Rnj）的广义特征值和特征向量%排序%最优权矢量a=exp(ima*2*pi*d_lamda*si n(theta(j)*pi/180)*0:eleme nt_num-1'); f(j)=Wopt'*a; p(j)=a

11、'*Rs*a+a'*Rnj*a;endF=20*log10(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)

12、;hold onsubplot(122)plot(theta,P,'r');grid on xlabel('theta/0');ylabel('功率 in dB'); title('max-SNR 功率谱')仿真图如下：5. %ASC旁瓣相消-MSE准则%辅助天线的数目clc;close all;clear all ima=sqrt(-1);%来波方向%干扰方向 %采样单元数 %预划分一个区域M=32; d_lamda=.5; theta0=-30; theta1=60;L=512; s=zeros(1,512); for ii=

13、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(thet a1*pi/180)+ampn*(randn(1,1)+ima*randn(1,1);% 接收信号(

14、信号 +干扰 +噪声) s0(ii)=amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180);endRx=1/L*jam*jam'%噪声自相关矩阵，相当于 X(t)r_xd=1/L*jam*s'Wopt=pinv(Rx)*r_xd;delta=s0-(s-Wopt'*jam);delta仁abs(mea n(delta42)-(mea n( delta)42)% 方差theta=linspace(-pi/2,pi/2,200);for jj=1:length(theta) a=exp(ima*2*pi*.5*sin(theta(jj)*0:M-

15、1'); f(jj)=Wopt'*a;endF=20*log10(abs(f)/(max(max(abs(f);figure(1)plot(theta*180/pi,F),grid on,hold onplot(theta0,-50:0,'.')plot(theta1,-50:0,'.')xlabel('theta/o');ylabel('F/dB');title('MSE 准则下的方向图 ')axis(-90 90 -50 0);%可为x轴和y轴设置一个极限范围仿真图如下：MSF准则下的方向圈-3

16、0-6040-2002040 GO SOthata/o6. %线性约束最小方差（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*

17、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); endRx=1/L*x* x' steer1=exp(ima*2*pi*1/2*sin(theta0*pi/180)*0:e

18、lement_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)*

19、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');titl

20、e('Capon beamforming 方向图 ') axis(-90 90 -20 0);%可为 x 轴和 y 轴设置一个极限范围 P=10*log10(abs(p)/(max(max(abs(p);subp_oo22)pofahefpp)込id opho-d on p-ofaheQr?20o< p-ofaheQrl20o< p-ofaheQr220o_) x-abeQhefa、o-)八 y_abe_(&4M、dB) m-e(oapon beamforming&4MSS-)ss.F/dBIXJ ocq mhG onihCDno0功率/dE45cm

21、pcln beamfclmrnngft30 swleu/o008. %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;%阵元数%阵元间距与波长的关系%范围%来波方向%干扰方向%干扰方向% 采样单元数ampO=1O*ra ndn;信号的幅度随机产生，保证信号之间是不相关的amp1=200*randn(1);amp2=2OO*randn(1);a

22、mpn=3;x(:,i)=ampO*exp(ima*2*pi*1/2*sin(thetaO*pi/18O)*O:element_num-1')+.amp1*exp(ima*2*pi*1/2*sin(theta1*pi/18O)*O:element_num-1')+.amp2*exp(ima*2*pi*1/2*sin(theta2*pi/18O)*O:element_num-1')+. ampn*(randn(element_num,1)+ima*randn(element_num,1);endRx=1/L*x* x'R=inv(Rx);steer=exp(ima

23、*2*pi*1/2*sin(thetaO*pi/18O)*O: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/18O)*O:element_num-1'); f(j)=w'*a;p(j)=1/(a'*R*a);endF=2O*log1O(abs(f)/(max(max(abs(f);P=20*log10(abs(p)/(max(max(abs(p);%此处是功率的对数形式

24、 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('t

25、heta/o');ylabel('功率 /dB');title('C apon beamformi ng 功率谱') 仿真图如下：Capon波束形成9. %不同方法估计协方差矩阵的 clc;clear all;close all;%阵元数为 8%间距为半波长 %范围 %来波方向 %干扰方向ima=sqrt(-1); element_num=8; d_lamda=1/2; theta=-90:0.5:90; theta0=0; theta1=50;L=512;%采样单元数for i=1:Lamp0=10*randn(1);amp1=50*randn(1);

26、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');endRx=1/L*(s+j+n)*(s+j+n)'% 接收信号自相关矩阵Rnj=1/L*(j+n)*(j+n)'% 干拢+噪声的自相关矩阵 e=exp(ima*2*pi*d

27、_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;endF1=20

28、*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 onplot(theta0,-50:0,'.');plot(theta1,-50:0,'.')xlabel('theta/0');ylabel('F(1,2)/dB');title(' 不同方法估计协方差矩阵的 Capon 波束形成 ');axis(-90 90 -50 0

29、);仿真图如下：二-10-1520-25-30-35-40-45不同方法估计协方差矩阵的4“门波束形成-80*60*40-20020406080theta/0-5010. %多点约束的Capon波束形成和方向图 clc;clear all ; close all;ima=sqrt(-1); eleme nt_num=8; d_l amda=1/2; theta=-90:0.3:90; theta0=0; theta1=20; theta2=50;L=512;%阵元数%阵元间距与波长的关系%搜索范围%来波方向%干扰方向%干扰方向%米样单兀数Rx=zeros(eleme nt_nu m,eleme

30、 nt_nu m);%P 生协方差矩阵for i=1:L;amp0=10*ra ndn (1);amp1=10*ra ndn (1);amp2=50*ra ndn (1);namp=0.5*randn(1);%噪声的幅度随机产生，保证噪声与信号之间是不相关的J(:,i)=amp1*exp(ima*2*pi*1/2*si n(theta1*pi/180)*0:eleme nt_num-1')+.amp2*exp(ima*2*pi*1/2*si n(theta2*pi/180)*0:eleme nt_num-1')+.n amp*exp(ima*2*pi*ra ndn (1)*0:

31、eleme nt_nu m-1'); x(:,i)=amp0*exp(ima*2*pi*1/2*si n(theta0*pi/180)*0:eleme nt_num-1')+.J(:,i);%表示接收信号endRx=Rx+1/L*x*x'R=i nv(Rx);w=amp0*exp(ima*2*pi*1/2*si n(theta0*pi/180)*0:eleme nt_num-1')+. amp2*exp(ima*2*pi*1/2*si n(theta2*pi/180)*0:eleme nt_num-1')+.amp1*exp(ima*2*pi*1/2*s

32、i n(theta1*pi/180)*0:eleme nt_num-1');for j=1:le ngth(theta);a=exp(ima*2*pi*d_lamda*si n( theta(j)*pi/180)*0:eleme nt_num-1'); f(j)=w'*a;p(j)=1/(a'*R*a);endF=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

33、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

34、:0,'.'); xlabel('theta/o');ylabel('P/dB');title('cap on beamformi ng 功率谱') 仿真图如下：capon beamforiniiig 向圏capon bemmfcirming功率谱theta-otheta/o大作业%自适应波束形成方向图clc;clear all;close all;ima=sqrt(-1);element_num=8;%阵元数c=3e8;f=5OOe6;lamda=c/f;d=1/2*lamda;%阵元间距与波长的关系theta=-9O:O.5:

35、9O;%范围thetaO=O;%来波方向theta1=45;%干扰方向theta2=6O;%干扰方向L=512;% 采样单元数for i=1:L;ampO=1O*ra ndn;信号的幅度随机产生，保证信号之间是不相关的amp1=100*randn(1);amp2=1OO*randn(1);ampn=1O;x(:,i)=ampO*exp(ima*2*pi*1/2*sin(thetaO*pi/18O)*O:element_num-1')+. amp1*exp(ima*2*pi*1/2*sin(theta1*pi/18O)*O:element_num-1')+. amp2*exp(i

36、ma*2*pi*1/2*sin(theta2*pi/18O)*O:element_num-1')+.amp n*(ra ndn (eleme nt_nu m,1)+ima*ra ndn (eleme nt_nu m,1);%产 生阵列信号 endRx=1/L*x* x'steer1=exp(ima*2*pi*1/2*sin(thetaO*pi/18O)*O: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(t

37、heta2*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);endF=20*log10(abs(f)/(max(max(abs(f);P=20*log10(abs(p)/(max(max(abs(p);%此处是功率的对数形式figuresubplot(121)plot(theta,P),grid on,hold onxlabel('theta/o');ylabel('F/dB');title(' 阵列信号方向图 &#