波束形成 Matlab程序教学文稿
- 格式:doc
- 大小:1.21 MB
- 文档页数:33
信息与通信工程学院阵列信号处理实验报告(自适应波束形成Matlab仿真)学号:XXXXXX专业:XXXXXX学生姓名:XXX任课教师:XXX2015年X月题目:自适应波束形成Matlab 仿真1. 算法简述:自适应波束形成,源于自适应天线的一个概念。
接收端的信号处理,可以通过将各阵元输出进行加权求和,将天线阵列波束“导向”到一个方向上,对期望信号得到最大输出功率的导向位置即给出波达方向估计。
波束形成算法是在一定准则下综合个输入信息来计算最优权值的数学方法,线性约束最小方差准则(LCMV )是最重要、最常用的方法之一。
LCMV 是对有用信号形式和来向完全已知,在某种约束条件下使阵列输出的方差最小。
该准则属于广义约束,缺点是需要知道期望分量的波达方向。
准则的代价函数为Rw w w J H )(=,约束条件为H ()θ=w a f ;最佳解为f c R c c R w 11H 1H ][---。
2. 波束形成原理以一维M 元等距离线阵为例,如图1所示,设空间信号为窄带信号,每个通道用一个附加权值系数来调整该通道的幅度和相位。
图1 波束形成算法结构图这时阵列的输出可以表示为:*1()()()Mi i iy t w x t θ==∑ 如果采用矢量来表示各阵元输出及加权系数,即T 12()[()()()]M x t x t x t x t = T 12()[()()()]M w w w w θθθθ=()x n()x n()x n那么,阵列的输出也可以用矢量表示为H ()()()y t t θ=w x为了在某一方向θ上补偿各阵元之间的时延以形成一个主瓣,常规波束形成器在期望方向上的加权矢量可以构成为(1)T ()[1e e ]j j M w ωτωτθ---=观察此加权矢量,发现若空间只有一个来自方向θ的信号,其方向矢量()αθ的表示形式与此权值矢量相同。
则有H H ()()()()()y t t t θαθ==w x x这时常规波束形成器的输出功率可以表示为2H H ()[()]()()()()CBF P E y t θθθαθαθ===w Rw R式中矩阵R 为阵列输出()t x 的协方差矩阵。
1?均匀线阵方向图%8阵元均匀线阵方向图,来波方向为clc; clear all; close all;0度imag=sqrt(_1);element_num=8;% 阵元数为8d_lamda=1/2;%阵元间距d与波长lamda的关系theta=li nspace(-pi/2,pi/2,200);theta0=0;% 来波方向w=exp(imag*2*pi*dl_lamda*sin(theta0)*[0:eleme nt_nu m-1]');for j=1:le ngth(theta)a=exp(imag*2*pi*dd_l amda*si n(theta(j))*[0:eleme nt_nu m-1]'); p(j)=w'*a;end figure;plot(theta,abs(p)),grid on xlabel('theta/radia n')ylabel('amplitude')title('8 阵元均匀线阵方向图')°2 8阵元均匀线阵方向图7654321-15 -1 -0 5 0 06thetaradian1 158当来波方向为45度时,仿真图如下8阵元均匀线阵方向图如下,来波方向为0 度,20log (dB)8阵元均苛銭阵方向图来波方向为0度随着阵元数的增加,波束宽度变窄,分辨力提高:仿真图如下Q d pE =ro 二2. 波束宽度与波达方向及阵元数的关系clcclear allclose allima=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; forj=1:length(theta);fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_n um1*d));psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_n um2*d));beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_ num3*d));endfigure; plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid on xlabel('theta'); ylabel('Width in radians')title(' 波束宽度与波达方向及阵元数的关系') 仿真图如下:。
1.均匀线阵方向图(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=32;%阵元数为8d_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;endpatternmag=abs(p);patternmagnorm=patternmag/max(max(patternmag));patterndB=20*log10(patternmag);patterndBnorm=20*log10(patternmagnorm);figure(1)plot(theta*180/pi,patternmag);grid on;xlabel('theta/radian')ylabel('amplitude/dB')title([num2str(element_num) '阵元均匀线阵方向图','来波方向为' num2str(theta0*180/pi) '度']);hold on;figure(2)plot(theta,patterndBnorm,'r');grid on;xlabel('theta/radian')ylabel('amplitude/dB')title([num2str(element_num) '阵元均匀线阵方向图','来波方向为' num2str(theta0*180/pi) '度']);axis([-1.5 1.5 -50 0]);(2)仿真结果A.来波方向为0°不归一化归一化B.来波方向为45°不归一化归一化C.随着阵元数的增加,波束宽度变窄,分辨力提高,仿真图如下:非归一化归一化不归一化归一化2.波束宽度与波达方向及阵元数的关系(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num1=16;element_num2=128;element_num3=1024;lambda=0.1;d=0.5*lambda;theta=0:0.5:90;for j=1:length(theta)fai(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num1*d)); psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num2*d)); beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num3*d)); endfigureplot(theta,fai,'r',theta,psi,'b',theta,beta,'g');grid on;xlabel('theta');ylabel('width in radians');title('波束宽度与达波方向及阵元数目的关系');legend('N=16','N=128','N=1024');(2)仿真结果结果3. 当阵元间距/2d λ>时,会出现栅瓣,导致空间模糊(1)仿真结果非归一化归一化4. 类似于时域滤波,天线方向图是最优权的傅立叶变换(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=32;source_num=1;d_lambda=0.5;theta=linspace(-pi/2,pi/2,200);theta0=0;w=exp(imag*2*pi*d_lambda*sin(theta0)*[0:element_num-1]');for j=1:length(theta)a=exp(imag*2*pi*d_lambda*sin(theta(j))*[0:element_num-1]'); p(j)=w'*a;endpatternmag=abs(p);patternmagnorm=patternmag/max(max(patternmag));patterndB=20*log10(patternmag); patterndBnorm=20*log10(patternmagnorm); figure(1)subplot(1,2,1);plot(theta,patterndBnorm);grid on;xlabel('theta/radian');ylabel('amplitude/dB');axis([-2.0 2.0 -50 0]);subplot(1,2,2);pfft=fftshift(fft(w,256));pfftmag=abs(pfft);pfftmagnorm=pfftmag/max(max(pfftmag)); pfftdB=20*log10(pfftmagnorm);pfftdBnorm=20*log10(pfftmagnorm);plot(linspace(-pi/2,pi/2,256),pfftdBnorm); grid on;xlabel('theta/radian');ylabel('FFT_amplitude/dB');axis([-2.0 2.0 -50 0]);(2)仿真结果5.最大信噪比准则方向图和功率谱(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_lambda=0.5;%间距为半波长theta=-90:0.5:90;%扫描范围theta0=0;%来波方位theta1=20;%干扰方向L=512;%采样点数for i=1:Lamp0=10*randn(1);amp1=200*randn(1);ampn=1;s(:,i)=amp0*exp(imag*2*pi*0.5*sin(theta0*pi/180)*[0:element_num-1]');j(:,i)=amp1*exp(imag*2*pi*0.5*sin(theta1*pi/180)*[0:element_num-1]');n(:,i)=ampn*exp(randn(element_num,1)+imag*randn(element_num,1)); endRs=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(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*[0:element_num-1]');f(j)=Wopt'*a;p(j)=a'*Rs*a+a'*Rnj*a;endF=20*log10(abs(f)/max(max(abs(f))));P=20*log10(abs(p)/max(max(abs(p))));subplot(1,2,1)plot(theta,F);grid on;hold on;plot(theta0,-80:0,'.'); plot(theta1,-80:0,'.'); xlabel('theta/0'); ylabel('F in dB');title('max-SNR 方向图'); axis([-90 90 -80 0]); hold on;subplot(1,2,2);plot(theta,P,'r');grid on;xlabel('theta/0'); ylabel('功率 in dB'); title('max-SNR 功率谱'); grid on;axis([-90 90 -80 0]); (2)仿真结果6.ASC旁瓣相消----MSE准则(1) matlab 程序clc;close all;clear all;imag=sqrt(-1);M=32;%辅助天线数目d_lambda=0.5;%阵元间距theta0=-30;%来波方向theta1=60;%干扰方向L=512;%采样单元数s=zeros(1,512); %预划分一个区域for ii=1:Lamp0=1*randn(1);%信号的幅度随机产生,保证信号之间是不相关的amp1=200*randn(1);ampn=1;jam(:,ii)=amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:M-1]' )+ampn*(randn(M,1)+imag*randn(M,1)); %干扰+噪声s(ii)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180))+amp1*exp(im ag*2*pi*d_lambda*sin(theta1*pi/180))+ampn*(randn(1,1)+imag*randn( 1,1));%接收信号(信号+干扰+噪声)s0(ii)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180));endRx=1/L*jam*jam';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(imag*2*pi*d_lambda*sin(theta(jj))*[0:M-1]'); f(jj)=Wopt'*a;endF=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/°');ylabel('F/dB');title('MSE准则下的方向图');axis([-90 90 -50 0]);(2)仿真结果7.线性约束最小方差(LCMV)准则(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数d_lambda=0.5;%阵元间距与波长的关系theta=-90:0.5:90; %搜索范围theta0=0; %三个信号源的来波方向theta1=30;theta2=60;L=512;%采样单元数for i=1:Lamp0=10*randn(1);amp1=100*randn(1);amp2=10*randn(1);ampn=10;x(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]')+am p1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]')+amp2*exp(im ag*2*pi*d_lambda*sin(theta2*pi/180)*[0:element_num-1]')+ampn*(randn(element_ num,1)+imag*randn(element_num,1));endRx=1/L*x*x';steer1=exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]');steer2=exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]');steer3=exp(imag*2*pi*d_lambda*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(imag*2*pi*d_lambda*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)))));subplot(1,2,1)plot(theta,F);grid on;hold on;plot(theta0,-20:0,'.');plot(theta1,-20:0,'.');plot(theta2,-20:0,'.');xlabel('theta/°');ylabel('F/dB');title('Capon beamforming 方向图');axis([-90 90 -20 0]);P=20*log10(abs(p)/(max(max(abs(p)))));subplot(1,2,2)plot(theta,P);grid on;hold on;plot(theta0,-20:0,'.');plot(theta1,-20:0,'.');plot(theta2,-20:0,'.');xlabel('theta/°');ylabel('P/dB');title('Capon beamforming 功率谱');axis([-90 90 -20 0]);(2)仿真结果8.Capon beamforming(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数d_lambda=0.5;%阵元间距与波长的关系theta=-90:0.5:90; %搜索范围theta0=0; %三个信号源的来波方向theta1=20;theta2=60;L=1000;%采样单元数for i=1:Lamp0=10*randn(1);amp1=200*randn(1);amp2=200*randn(1);ampn=3;x(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]')+am p1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]')+amp2*exp(im ag*2*pi*d_lambda*sin(theta2*pi/180)*[0:element_num-1]')+ampn*(randn(element_ num,1)+imag*randn(element_num,1));endRx=1/L*x*x';R=inv(Rx);steer=exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]');w=R*steer/(steer'*R*steer);%最优权矢量for j=1:length(theta)a=exp(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*[0:element_num-1]');f(j)=w'*a;p(j)=1/(a'*R*a);endF=20*log10(abs(f)/(max(max(abs(f)))));subplot(1,2,1)plot(theta,F);grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');plot(theta2,-50:0,'.');xlabel('theta/°');ylabel('F/dB');title('Capon beamforming 方向图');axis([-90 90 -50 0]);P=20*log10(abs(p)/(max(max(abs(p)))));subplot(1,2,2)plot(theta,P);grid on;hold on;xlabel('theta/°');ylabel('P/dB');title('Capon beamforming 功率谱');axis([-90 90 -90 0]);(2)仿真结果9.不同方法估计协方差矩阵的Capon波束形成(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_lambda=0.5;%间距为半波长theta=-90:0.5:90;%扫描范围theta0=0;%来波方向theta1=50;%干扰方向L=1024;%采样单元数for i=1:Lamp0=10*randn(1);amp1=50*randn(1);ampn=0.5;s(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]');j(:,i)=amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]'); n(:,i)=ampn*exp(imag*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(imag*2*pi*d_lambda*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(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*[0:element_num-1]');f1(j)=Wopt_Rx'*a;f2(j)=Wopt_Rnj'*a;endF1=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/°');ylabel('F(1,2)/dB');title('不同方法估计协方差矩阵的Capon波束形成');axis([-90 90 -60 0]);(2)仿真结果10.多点约束的Capon波束形成和方向图(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;d_lambda=0.5;theta=-90:0.3:90;theta0=0;theta1=20;theta2=50;L=512;Rx=zeros(element_num,element_num);%产生协方差矩阵for i=1:Lamp0=10*randn(1);amp1=10*randn(1);amp2=50*randn(1);ampn=0.5*randn(1);%噪声的幅度随机产生,保证噪声与信号之间是不相关的j(:,i)=amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]') +amp2*exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*[0:element_num-1]')+ampn*exp (imag*2*pi*randn(1)*[0:element_num-1]');x(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]') +j(:,i);%表示接收信号endRx=1/L*x*x';R=inv(Rx);w=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]')+amp1*ex p(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]')+amp2*exp(imag*2* pi*d_lambda*sin(theta2*pi/180)*[0:element_num-1]');for j=1:length(theta)a=exp(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*[0:element_num-1]');f(j)=w'*a;p(j)=1/(a'*R*a);endF=20*log10(abs(f)/max(max(abs(f))));P=20*log10(abs(p)/max(max(abs(p))));figure;subplot(1,2,1);plot(theta,F);grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');plot(theta2,-50:0,'.');xlabel('theta/°');ylabel('F/dB');title('Capon beamforming方向图'); axis([-90 90 -50 0]);subplot(1,2,2);plot(theta,P);hold on;grid on;plot(theta0,-90:0,'.');plot(theta1,-90:0,'.');plot(theta2,-90:0,'.');xlabel('theta/°');ylabel('P/dB');title('Capon beamforming功率谱');(2)仿真结果11.自适应波束形成方向图(1)matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;c=3e8;f=5e8;lambda=c/f;d_lambda=0.5;theta=-90:0.5:90;theta0=0;theta1=45;theta2=60;L=2048;for i=1:Lamp0=10*randn(1);amp1=100*randn(1);amp2=100*randn(1);ampn=10;x(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]')+am p1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]')+amp2*exp(im ag*2*pi*d_lambda*sin(theta2*pi/180)*[0:element_num-1]')+ampn*(randn(element_ num,1)+imag*randn(element_num,1));endRx=1/L*x*x';steer1=exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*[0:element_num-1]');steer2=exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*[0:element_num-1]');steer3=exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*[0:element_num-1]');C=[steer1 steer2 steer3];F=[1 0 0]';w=inv(Rx)*C*(inv(C'*inv(Rx)*C))*F;for j=1:length(theta)a=exp(imag*2*pi*d_lambda*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)))));subplot(1,2,1)plot(theta,F);grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');plot(theta2,-50:0,'.');xlabel('theta/°');ylabel('F/dB');title('自适应波束形成方向图');axis([-90 90 -50 0]);P=20*log10(abs(p)/(max(max(abs(p))))); subplot(1,2,2)plot(theta,P);grid on;hold on;xlabel('theta/°');ylabel('P/dB');title('功率谱');axis([-90 90 -50 0]);(2)仿真结果(3)GUI界面。
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=0;%来波方向w=exp(imag*2*pi*dl_lamda*si n(theta0)*[0:eleme nt_nu m-1]');for j=1:le ngth(theta)a=exp(imag*2*pi*dd_l amda*si n(theta(j))*[0:eleme nt_nu m-1]'); p(j)=w'*a;endfigure;plot(theta,abs(p)),grid onxlabel('theta/radia n')ylabel('amplitude')title('8阵元均匀线阵方向图')8阵元均匀线阵方向图当来波方向为45度时,仿真图如下:8阵元均匀线阵方向图如下,来波方向为0度,20log (dB )8阵元均苛銭阵方向图来波方向为0度 S =-s ==d E B随着阵元数的增加,波束宽度变窄,分辨力提高:仿真图如下: Qp=二 d E ro2. 波束宽度与波达方向及阵元数的关系clcclear allclose allima=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_n um1*d));psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_n um2*d));beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lamda/(element_ num3*d));endfigure; plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid on xlabel('theta');ylabel('Width in radians')title(' 波束宽度与波达方向及阵元数的关系') 仿真图如下:。
宽带波束形成matlab,关于均匀圆阵MVDR宽带波束形成的程序%⽤均匀圆阵MVDR⽅法进⾏宽带波束形成%% 基本参数f0 = 4*10^6; %信号中⼼频率bandwide = 0.5*10^6; %带宽0.05MHzinterval = 100; %每间隔interval(Hz)形成点频信号,合成宽带NarrowNumber = bandwide/interval; %所需点频个数NarrowF = linspace(f0-bandwide/2, f0+bandwide/2, NarrowNumber); %各点频的频率值R = 50; %阵元半径c = 3*10^8; %传播速度:光速bosu=340; % 波速degrad=pi/180;p=1;seta1=150*degrad; %信号到来⽅位⾓fai1=45*degrad; %信号到来仰⾓kk = 64; %进⾏MVDR所需采样数Nfft = 256; %进⾏FFT⼦带化点数,即进⾏Nfft点快速傅⽴叶变换Nlearn = kk * Nfft; %时域数据快拍数M=40; %均匀圆阵总数MElementFai=(0:(M-1))*9*degrad+4.5*degrad; %各阵元位置向量图MNumber = 10; %选出的阵元个数SelectElement=[ElementFai(1) ElementFai(5) ElementFai(9) ElementFai(13) ElementFai(17) ElementFai(21) ElementFai(25) ElementFai(29) ElementFai(33) ElementFai(37)];A=zeros(MNumber,p);k=[0:MNumber-1]';for t=1:pA(:,t)=exp(-j*2*pi*R*f0/bosu*(cos(seta1*degrad)*cos(2*pi*k/M)+cos(seta1*degrad)*sin(2*pi*k/M)));%圆阵的阵列流型ends=zeros(p,Nlearn);Am=[3,5]; %幅度phi=[20,30]; %初相for ii=1:ps(ii,:)=s(ii,:)+Am(ii)*cos(2*pi*f0*t+phi(ii)*degrad);%宽带信号endF=linspace(f0-bandwide/2, f0+bandwide/2, Nfft);snr=10;n=exp(j*2*pi*randn(MNumber,Nlearn))/snr;%噪声X=A*s + n; % 阵列输出%% 开始处理%对时域数据分段作FFT变换for mn=1:MNumberfor m=1:kkxxF(mn,m,:) = fft( X(mn,((m-1)*Nfft+1):((m-1)*Nfft+Nfft)), Nfft);endendtemp1=cos(seta1-SelectElement);F=linspace(f0-bandwide/2, f0+bandwide/2, Nfft);%求窄带信号的⾃适应最佳权for nfft=1:Nffth=zeros(MNumber);for k=1:kkh=h+xxF(:,k,nfft)*xxF(:,k,nfft)';endRxxf(:,:,nfft)=h./kk; %得到某窄频的⾃相关矩阵windows=ones(MNumber,1);a1=exp(-j*2*pi*R*F(nfft)/c*sin(fai1)*temp1);s=(a1.*windows.').';Wopt(:,nfft)=(inv(Rxxf(:,:,nfft))*s)/(s'*inv(Rxxf(:,:,nfft))*s); %MVDR算法的程序end%% 扫描⽮量NcurveAzi=180; %空间扫描⾓度采样数seta=linspace(0,2*pi,NcurveAzi); %阵元空间扫描指向⾓度for ne=1:Nfftfor na=1:NcurveAzitemp=cos(seta(na)-SelectElement);a=exp(-j*2*pi*R*F(ne)/c*sin(fai1)*temp);%扫描时⽅向⽮量f=(1+temp)/2; %扫描时的阵元⽅向图BF(na,ne)=abs(sum(a.*Wopt(:,nfft)'));end%%仿真结果plot(seta*180/pi,20*log10(BF(:,1)/max(BF(:,1))),'b') gridxlabel('\phi⽅位⾓(度)')ylabel('波束输出/分贝')title('均匀圆阵MVDR算法宽带波束形成图')endhold off。
Ca p o n - 波束形成ma t l a b 仿真(附源代码)Capon波束形成阵列N=16,信号o -30 ,干扰为! 60 , 2 19 , 3 45 ,干扰功率分别为:40dB,35dB,50dB。
Capon波束形成后的方向图和功率谱如下:-10-20-30E -40B0 1D0Capon beamForming \门;[7Capon beftinfonning 勺牛吃fo o o -7-6-9Co o o o o o12 3 4 5-BO-80-100 3 -60 40 '20 0 20 4Q &0tTieta>o-60 40 -20 0 20 40 60 00 100ttiefa/o为了比较接收数据直接估计噪声协方差矩阵和利用干扰+噪声估计协方差矩阵的Capon波束形成的差异,进行如下仿真:可以看出利用干扰+噪声估计协方差矩阵的方向图性能较优于接收数据直接估计噪声协方差矩阵的方向图。
代码:clc;clear all ;close all;ima=sqrt(-1);eleme nt_num=8; %阵元数d_lamda=1/2; %阵元间距与波长的关系theta=-90:0.5:90; %范围theta0=-30; %来波方向theta仁60; %干扰方向1theta2=19; %干扰方向2theta3=45; %干扰方向3■30買E-O-so-90-80 ^0 -40 -20 0 20 40 GO B0theta/O不冋方也佶计协方垫年阵的Capon泼束形成5_-L=1000; %采样单元数for i=1:L;amp0=10*ra ndn( 1);%信号的幅度随机产生,保证信号之间是不相关的amp1=100*ra ndn (1);%输入阵列的噪声amp2=sqrt(10A3.5)*ra ndn( 1);%输入阵列的噪声amp3=sqrt(10A5)*ra ndn⑴;%输入阵列的噪声amp n=3;% 噪声x(:,i)=amp0*exp(ima*2*pi*1/2*si n(theta0*pi/180)*[0:eleme nt_num-1]')+...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]')+...amp3*exp(ima*2*pi*1/2*si n(theta3*pi/180)*[0:eleme nt_num-1]')+...amp n*(ra ndn (eleme nt_nu m,1)+ima*ra ndn (eleme nt_nu m,1));endRx=1/L*x* x';R=i nv(Rx);steer=exp(ima*2*pi*1/2*si n(theta0*pi/180)*[0:eleme nt_num-1]');w=R*steer/(steer'*R*steer);%Capon 最优权矢量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=20*log10(abs(f)/(max(abs (f))));P=20*log10(abs(p)/(max(abs(p))));%此处是功率的对数形式figure;% subplot(121)plot(theta,F),grid on,hold onplot(theta0,-80:0,'.')plot(theta1,-80:0,'.')plot(theta2,-80:0,'.')plot(theta3,-80:0,'.') xlabel('theta/o'); ylabel('F/dB');title('Capon beamforming方向图') % axis([-90 90 -50 0]); % subplot(122)figure;plot(theta,P),grid on xlabel('theta/o'); ylabel('功率/dB'); title('C apon beamformi ng 功率谱')。
在进行均匀圆阵二维波束形成的Matlab编程前,我们首先要了解什么是均匀圆阵二维波束形成。
均匀圆阵是由均匀分布的天线组成,这些天线排列成圆环或者圆盘状,用来传输和接收无线电波。
而二维波束形成是指通过合理控制每个天线的相位和幅度,将传输的信号波束聚焦在特定方向。
在现实生活中,均匀圆阵二维波束形成广泛应用于雷达、通信系统、声纳等领域,能够提高信号的接收性能和抗干扰能力。
接下来,让我们来探讨如何使用Matlab进行均匀圆阵二维波束形成的编程。
我们需要了解Matlab中涉及到的一些基本概念和函数。
在Matlab中,我们可以使用向量和矩阵来表示天线的位置和信号的相位和幅度。
通过Matlab编程,我们可以方便地对天线的布局和相位进行优化,以实现二维波束形成。
另外,Matlab提供了丰富的信号处理工具和绘图函数,可以帮助我们分析和可视化波束形成的效果。
让我们通过一个简单的示例来演示如何在Matlab中实现均匀圆阵二维波束形成。
假设我们有一个均匀圆阵,包含8个天线,我们希望将信号波束聚焦在特定方向。
我们需要使用Matlab创建一个包含8个天线坐标的矩阵,并根据波束形成的需求计算每个天线的相位。
我们可以利用Matlab提供的信号处理函数,对每个天线的输出信号进行加权叠加,从而实现波束形成。
我们可以利用Matlab的绘图函数,将波束形成的效果直观地展示出来。
在实际应用中,除了基本的波束形成技术,我们还可以结合自适应波束形成、干涉波束形成等高级技术,来进一步提高系统的性能。
在Matlab中,这些高级技术同样可以通过编程实现,为我们提供更加灵活和强大的工具来应对复杂的实际场景。
总结回顾起来,Matlab是一个非常强大的工具,能够帮助我们实现均匀圆阵二维波束形成。
通过合理的编程和算法设计,我们可以实现对天线相位和幅度的精确控制,将信号波束聚焦在所需的方向,从而提高系统的性能和可靠性。
在实际工程应用中,我们还可以结合Matlab 提供的其他工具和函数,来进一步完善系统设计和性能优化。
波束形成与智能天线1.均匀线阵方向图%8阵元均匀线阵方向图,来波方向为0度clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_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;endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')当来波方向为45度时,仿真图如下:8阵元均匀线阵方向图如下,来波方向为0度,20log(dB)随着阵元数的增加,波束宽度变窄,分辨力提高:仿真图如下:2.波束宽度与波达方向及阵元数的关系clcclear allclose allima=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));endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下:3.当阵元间距/2dλ>时,会出现栅瓣,导致空间模糊。
波束形成与智能天线1.均匀线阵方向图%8阵元均匀线阵方向图,来波方向为0度clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_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;endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')当来波方向为45度时,仿真图如下:8阵元均匀线阵方向图如下,来波方向为0度,20log(dB)随着阵元数的增加,波束宽度变窄,分辨力提高:仿真图如下:2.波束宽度与波达方向及阵元数的关系clcclear allclose allima=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/-asin(sin(theta(j)*pi/)-lamda/(element_num1*d));psi(j)=theta(j)*pi/-asin(sin(theta(j)*pi/)-lamda/(element_num2*d));beta(j)=theta(j)*pi/-asin(sin(theta(j)*pi/)-lamda/(element_num3*d)); endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下:3. 当阵元间距时,会出现栅瓣,导致空间模糊。
波束形成M a t l a b程序1.均匀线阵方向图%8阵元均匀线阵方向图,来波方向为0度clc;clear all;close all;imag=sqrt(-1);element_num=8;%阵元数为8d_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;endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')当来波方向为45度时,仿真图如下:8阵元均匀线阵方向图如下,来波方向为0度,20log(dB)随着阵元数的增加,波束宽度变窄,分辨力提高:仿真图如下:2.波束宽度与波达方向及阵元数的关系clcclear allclose allima=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)); endfigure;plot(theta,fai,'r',theta,psi,'b',theta,beta,'g'),grid onxlabel('theta');ylabel('Width in radians')title('波束宽度与波达方向及阵元数的关系')仿真图如下:3. 当阵元间距 /2dλ> 时,会出现栅瓣,导致空间模糊。
仿真图如下:4.类似于时域滤波,天线方向图是最优权的傅立叶变换仿真程序和仿真图如下:clcclear allclose allima=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;endfigure;subplot(1,2,1)plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('按定义的方向图')pfft=fftshift(fft(w,128));subplot(1,2,2)plot(linspace(-pi/2,pi/2,128),abs(pfft)),grid onxlabel('theta/radian')ylabel('FFT_amplitude')title('最优权的傅里叶变换')5.%最大信噪比准则方向图和功率谱clc;clear all;close all;ima=sqrt(-1);element_num=8; %阵元数为8d_lamda=1/2; %间距为半波长theta=-90:0.5:90; %范围theta0=0; %来波方向theta1=20; %干扰方向L=512; %采样单元数for i=1:Lamp0=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));endRs=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;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 onplot(theta0,-50:0,'.');plot(theta1,-50:0,'.')xlabel('theta/0');ylabel('F in dB');title('max-SNR 方向图');axis([-90 90 -50 0]);hold onsubplot(122)plot(theta,P,'r');grid onxlabel('theta/0');ylabel('功率 in dB');title('max-SNR功率谱')仿真图如下:6.%ASC旁瓣相消----MSE准则clc;close all;clear allima=sqrt(-1);M=32; %辅助天线的数目d_lamda=.5;theta0=-30; %来波方向theta1=60; %干扰方向L=512; %采样单元数s=zeros(1,512); %预划分一个区域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(thet a1*pi/180))+ampn*(randn(1,1)+ima*randn(1,1));%接收信号(信号+干扰+噪声) 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);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;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轴设置一个极限范围仿真图如下: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));endRx=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);endf=10*log10(abs(f)/(max(max(abs(f)))));figure(1)subplot(121)plot(theta,f),grid on,hold onplot(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 onplot(theta0,-20:0,'.')plot(theta1,-20:0,'.')plot(theta2,-20:0,'.')xlabel('theta/o');ylabel('功率/dB');title('Capon beamforming功率谱')仿真图如下:8. %Capon beamformingClc;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));endRx=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);endF=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 onplot(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 onxlabel('theta/o');ylabel('功率/dB');title('Capon beamforming功率谱') 仿真图如下:9. %不同方法估计协方差矩阵的Capon波束形成clc;clear all;close all;ima=sqrt(-1);element_num=8; %阵元数为8d_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);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_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*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]);仿真图如下: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);%表示接收信号endRx=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);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 onplot(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 onplot(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));%产生阵列信号endRx=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);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('阵列信号方向图')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;时仿真图如下:。