波束形成-Matlab程序
自适应波束形成及Matlab程序代码注解
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(theta 1*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(imag*2*pi*d_lambda*sin(th eta1*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]')+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*(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]')+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*(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_lam bda*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*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]'); 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]')+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*(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界面。
均匀圆阵波束形成matlab
均匀圆阵波束形成matlab
在MATLAB中，可以使用meshgrid函数生成均匀圆阵波束的坐标矩阵，然后使用波数表达式计算波束的强度。
以下是生成均匀圆阵波束并绘制其强度图的示例代码：
设置波长lambda, 波数k, 中心点x0, y0, 和圆阵半径R
lambda = 1;
k = 2;
x0 = 0; y0 = 0;
R = 5;
生成均匀圆阵的网格坐标
[x, y] = meshgrid(-R:R, -R:R);
计算波数强度
E = exp(1i * (k * (x - x0) + lambda * (y - y0)));
绘制波数强度图
surf(x, y, abs(E)); % 绘制波数强度的模值图
xlabel('x');
ylabel('y');
zlabel('E(x, y)');
title('Uniform Circle Array Wavefront Strength');
设置颜色映射为热图
colormap(jet);
添加色条
colorbar;
这段代码首先定义了波长lambda、波数k以及圆阵的中心点x0、y0和半径
R。
然后使用meshgrid函数生成网格坐标。
接着根据波数表达式计算每个点上的波数强度E。
最后，使用surf函数绘制波数强度图，并使用热图作为颜色映射。
LMS波束形成代码
LMS波束形成代码（matlab）LMS算法的仿真程序：%lms 算法clear allclose allhold off%系统信道权数sysorder = 5 ;%抽头数N=1000;%总采样次数inp = randn(N,1);%产生高斯随机系列n = randn(N,1);[b,a] = butter(2,0.25);Gz = tf(b,a,-1);%逆变换函数h= [0.0976;0.2873;0.3360;0.2210;0.0964;];%信道特性向量y = lsim(Gz,inp);%加入噪声n = n * std(y)/(10*std(n));%噪声信号d = y + n;%期望输出信号totallength=size(d,1);%步长N=60 ; %60节点作为训练序列%算法的开始w = zeros ( sysorder , 1 ) ;%初始化for n = sysorder : Nu = inp(n:-1:n-sysorder+1) ;% u的矩阵y(n)= w' * u;%系统输出e(n) = d(n) - y(n) ;%误差if n < 20mu=0.32;elsemu=0.15;endw = w + mu * u * e(n) ;%迭代方程end%检验结果for n = N+1 : totallengthu = inp(n:-1:n-sysorder+1) ;y(n) = w' * u ;e(n) = d(n) - y(n) ;%误差endhold onplot(d)plot(y,'r');title('系统输出') ;xlabel('样本')ylabel('实际输出')figuresemilogy((abs(e))) ;% e的绝对值坐标title('误差曲线') ;xlabel('样本')ylabel('误差矢量')figure%作图plot(h, 'k+')hold onplot(w, 'r*')legend('实际权矢量','估计权矢量') title('比较实际和估计权矢量') ;axis([0 6 0.05 0.35])2.NLMS算法的仿真程序：%lms 算法clear allclose allhold off%系统信道权数sysorder = 5 ;%抽头数N=1000;%总采样次数inp = randn(N,1);%产生高斯随机系列n = randn(N,1);[b,a] = butter(2,0.25);Gz = tf(b,a,-1);%逆变换函数h= [0.0976;0.2873;0.3360;0.2210;0.0964;];%信道特性向量y = lsim(Gz,inp);%加入噪声n = n * std(y)/(10*std(n));%噪声信号d = y + n;%期望输出信号totallength=size(d,1);%步长N=60 ; %60节点作为训练序列%算法的开始w = zeros ( sysorder , 1 ) ;%初始化for n = sysorder : Nu = inp(n:-1:n-sysorder+1) ;% u的矩阵y(n)= w' * u;%系统输出r(n)=u'*u;%自相关矩阵e(n) = d(n) - y(n) ;%误差fai=.0001;%修正参数，为防止u'*u过小导致步长值太大而设置的if n < 20mu=0.32;elsemu=0.15;endw = w + mu * u * e(n)/(r(n)+fai) ;;%迭代方程end%检验结果for n = N+1 : totallengthu = inp(n:-1:n-sysorder+1) ;y(n) = w' * u ;e(n) = d(n) - y(n) ;%误差endhold onplot(d)plot(y,'r');title('系统输出') ;xlabel('样本')ylabel('实际输出')figuresemilogy((abs(e))) ;% e的绝对值坐标title('误差曲线') ;xlabel('样本')ylabel('误差矢量')figure%作图plot(h, 'k+')hold onplot(w, 'r*')legend('实际权矢量','估计权矢量') title('比较实际和估计权矢量') ;axis([0 6 0.05 0.35])3.RLS算法的仿真程序：% RLS算法randn('seed', 0) ;rand('seed', 0) ;NoOfData = 8000 ; % Set no of data points used for training Order = 32 ; % 自适应滤波权数Lambda = 0.98 ; % 遗忘因子Delta = 0.001 ; % 相关矩阵R的初始化x = randn(NoOfData, 1) ;%高斯随机系列h = rand(Order, 1) ; % 系统随机抽样d = filter(h, 1, x) ; % 期望输出% RLS算法的初始化P = Delta * eye ( Order, Order ) ;%相关矩阵w = zeros ( Order, 1 ) ;%滤波系数矢量的初始化% RLS Adaptationfor n = Order : NoOfData ;u = x(n:-1:n-Order+1) ;%延时函数pi_ = u' * P ;%互相关函数k = Lambda + pi_ * u ;K = pi_'/k;%增益矢量e(n) = d(n) - w' * u ;%误差函数w = w + K * e(n) ;%递归公式PPrime = K * pi_ ;P = ( P - PPrime ) / Lambda ;%误差相关矩阵w_err(n) = norm(h - w) ;%真实估计误差end ;% 作图表示结果figure ;plot(20*log10(abs(e))) ;%| e |的误差曲线title('学习曲线') ;xlabel('迭代次数') ;ylabel('输出误差估计') ;figure ;semilogy(w_err) ;%作实际估计误差图title('矢量估计误差') ;xlabel('迭代次数') ;ylabel('误差权矢量') ;4.自适应均衡器的仿真程序：% Illustration of the conventional RLS algorithm close all;W=2.9;Nexp=10;N=2000;Nmc=1; % Number of ensemble realizations M=11;%抽头系数lambda=0.99;%遗忘因子varv=0.001;%噪声方差h=zeros(3,1);%h的初始化er=zeros(N,Nmc);%er的初始化h(1)=0.5*(1+cos(2*pi*(1-2)/W));h(2)=0.5*(1+cos(2*pi*(2-2)/W));h(3)=0.5*(1+cos(2*pi*(3-2)/W));% 学习曲线hc=[0 h(1) h(2) h(3)]';n0=7;t=(1:N)';for i=1:Nmcy=sign(rand(N,1)-0.5);%输入信号v=sqrt(varv)*randn(N,1);%噪声信号x=filter(hc,1,y)+v;%信号混合x=[zeros(M-1,1);x];%x矩阵yd=zeros(N+M-1,1); %延迟信号初始化e=yd;yd(n0+M-1:N+M-1)=y(1:N-n0+1);% CRLS 算法% Initializationlambda=0.98;P=(10^-3)*eye(M,M);c=zeros(M,1);g=c;glambda=g;% 迭代范围for n=M:M+N-1xn=flipud(x(n-M+1:n));glambda=P*xn;alphal=lambda+conj(glambda')*xn;g=glambda/lambda;a(n)=1-conj(g')*xn; P=(P-g*conj(glambda'))/lambda;P=(P+P')/2;e(n)=yd(n)-conj(c')*xn;c=c+g*conj(e(n));endendeplot=e(M:M+N-1).^2;subplot(2,1,1), plot(t,abs(eplot)) ylabel('|e(n)|^2');xlabel('n');subplot(2,1,2), plot(t,a(M:M+N-1)); ylabel('\alpha(n)');xlabel('n');5.自适应陷波器的仿真程序：N=400; %总采样长度t=0:N-1; %时间的变化范围s=sin(2*pi*t/20); %输入信号A=0.5; %干扰信号的幅值fai=pi/3;%干扰信号的相移n=A*cos(2*pi*t/10+fai);%干扰信号x=s+n;%信号混合subplot(2,2,1);%作第一子图plot(t,s);subplot(2,2,2); %作第二子图plot(t,x);x1=cos(2*pi*t/10);x2=sin(2*pi*t/10);%初始化w1=0.1;w2=0.1;e=zeros(1,N);y=0;u=0.05;%迭代步长for i=1:Ny=w1*x1(i)+w2*x2(i);e(i)=x(i)-y;%误差信号w1=w1+u*e(i)*x1(i);%迭代方程w2=w2+u*e(i)*x2(i);%迭代方程endsubplot(2,2,3); %作第三子图plot(t,e);subplot(2,2,4); %作第四子图plot(t,s-e);。
(完整word版)自适应波束形成与Matlab程序代码注解
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界面。
波束成形matlab
1111
波束成形是指通过控制输入信号的相位和振幅来实现发射或接收方向性辐射的技术。
在MATLAB中，可以通过以下步骤实现波束成形：
1. 确定要形成的波束方向和形状。
一种常用的方法是使用线性阵列，通过调整各个阵元的相位来控制波束方向。
在MATLAB中，可以通过创建一个表示阵列的矢量，然后使用复数表示各个阵元的相位来实现。
2. 使用MATLAB的信号处理工具箱中的FFT函数对输入信号进行频谱分析，并利用其频率响应进行波束成形处理。
3. 在实现波束成形算法时，需要注意信号传播路径和损耗。
MATLAB提供了有关传播路径和损耗的模型和函数，例如通过计算信号在不同介质中的传播速度和衰减来模拟实际情况。
4. 除了传播路径和损耗外，还可以考虑到干扰信号的处理。
MATLAB提供了各种信号处理方法，如滤波、降噪和干扰抑制等，可以帮助进行波束成形时的干扰信号处理。
总体而言，MATLAB可以通过使用其信号处理工具箱和丰富的数值计算功能来实现波束成形。
在实现波束成形算法时，还可以考虑干扰信号处理等问题，以提高波束成形的效果和精度。
（完整word版）自适应波束形成与Matlab程序代码注解
（完整word版）⾃适应波束形成与Matlab程序代码注解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;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)); 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(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');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;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]')+amp1*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)))));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]')+amp1*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）仿真结果。
(完整版)波束形成-Matlab程序
(完整版)波束形成-M a t l a b程序-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN1.均匀线阵方向图%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=45/180*pi;%来波方向 (我觉得应该是天线阵的指向)%theta0=0;%来波方向w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');for j=1:length(theta) %(我认为是入射角度，即来波方向,计算阵列流形矩阵A) a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]');p(j)=w'*a; %(matlab中的'默认为共轭转置，如果要计算转置为w.'*a)endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')见张小飞的书《阵列信号处理的理论和应用2.3.4节阵列的方向图》当来波方向为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 λ> 时，会出现栅瓣，导致空间模糊。
波束形成 Matlab程序教学文稿
波束形成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λ> 时，会出现栅瓣，导致空间模糊。
球面波束形成matlab
球面波束形成matlab
在MATLAB中，你可以使用以下步骤来生成球面波束：
步骤1，首先，你需要定义球面波束的一些基本参数，比如波
束的中心位置、半径、波长、以及波束的振幅等。
这些参数将决定
波束的形成和特性。
步骤2，接下来，你可以创建一个三维网格来表示波束的传播
空间。
你可以使用meshgrid函数来生成X、Y和Z方向上的坐标点。
步骤3，然后，根据球面波束的数学表达式，你可以计算每个
网格点上的波束振幅。
这通常涉及到使用波束的振幅函数以及传播
距离等参数来计算。
步骤4，最后，你可以使用surf函数将计算得到的波束振幅数
据在三维空间中进行可视化。
这将帮助你直观地观察到球面波束的
形成情况。
需要注意的是，生成球面波束涉及到一定的数学计算和编程技巧，需要对MATLAB的基本语法和数学运算有一定的了解。
另外，根
据具体的波束形成需求，你可能需要对上述步骤进行适当的调整和扩展。
希望这些步骤能够帮助你在MATLAB中成功生成球面波束。
自适应波束形成与Matlab程序代码注解讲诉
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界面。
波束成形music算法matlab程序
波束成形music算法matlab程序一、程序概述本程序实现了一种基于Music算法的波束成形算法，可用于信号处理领域中的方向性信号处理。
Music算法是一种基于信号统计特性的自适应算法，能够根据信号的方向信息自动调整波束方向，从而实现信号的高灵敏度和低噪声接收。
二、程序代码```matlab%参数初始化fs=1000;%采样频率time=-1:1/fs:1;%时间向量theta=0;%信号方向N=length(time);%样本数d=10;%接收器间距theta_max=atan2(sin(theta),cos(theta));%最大信号方向mu=0.05;%最大信号功率因子lambda=0.3;%自适应步长因子theta_ini=theta-pi/2;%初始搜索范围R=sqrt(d^2+theta_ini^2);%距离矩阵U=zeros(N,N);%方向矩阵p=zeros(N,1);%接收器位置向量%Music算法主程序fork=1:Np(k)=R(k,round(rand()*N));%随机选择一个接收器位置endU(:,k)=exp(-1i*pi*R(:,k)'*p(:)/p(:).^2);%方向矩阵初始化forn=1:100%自适应迭代次数fork=theta_ini:pi:pi/lambda%搜索范围ifabs(U(:,k).^2)>mu%如果当前方向上信号功率大于最大信号功率因子，则更新方向矩阵U(:,k)=exp(-1i*pi*R(:,k)'*p(:)/p(:).^2)*U(:,k);%更新方向矩阵元素endendendU=U/sum(U);%归一化方向矩阵%波束成形输出S=U*signal;%将信号乘以方向矩阵得到输出信号```三、程序说明程序首先初始化了一些参数，包括采样频率、时间向量、信号方向、接收器间距、最大信号方向、最大信号功率因子、自适应步长因子等。
mVDR波束形成matlab程序
close allclear allclcc=1500;fs=10e3;T = 1;t = 0:1/fs:T;L=length(t);f=2000;w=2*pi*f;k=w/c;M=11; %阵元个数Nmid=1; %参考点d=3;%阵元间距m=[0:1:M-1];yi=zeros(M,1);% 返回一个M*1维的零矩阵zi=zeros(M,1);xi=m*d;xi=xi.';%各阵元坐标y1=12;x1=10;z1=10;% 声源位置，y轴指向声源平面Ric1=sqrt((x1-xi).^2+(y1-yi).^2+(z1-zi).^2);%声源至各阵元的距离M*1维Rn1=Ric1-Ric1(Nmid);%声源至各阵元与参考阵元的声程差矢量M*1维s1=cos(w*t);%参考阵元接收到的信号1*L维snr =20;Am= 10^(-snr/20);n1=Am*(randn(M,L)+j*randn(M,L));%各阵元噪声矢量p1=zeros(M,L);%M*L维for k1=1:Mp1(k1,:)=Ric1(Nmid)/Ric1(k1)*s1.*exp(-j*w*Rn1(k1)/c); %各阵元经过幅度衰减和相位延迟后接收到的信号，M*L维endp=p1+n1;%各阵元接收的声压信号矩阵M*LR=p*p'/L;%接收数据的自协方差矩阵M*MRP=inv(R);%求R的逆矩阵% ----------------------------------------------------------% 扫描范围step_x=0.1;step_z=0.1;y=y1;x=[0:step_x:20];z=[0:step_z:20];for k1=1:length(z) % 纵坐标for k2=1:length(x)Ri=sqrt((x(k2)-xi).^2+(y-yi).^2+(z(k1)-zi).^2); %该扫描点至各阵元的聚焦距离矢量Rn=Ri-Ri(Nmid);%扫描点到各阵元与参考阵元的程差矢量M*1b=exp(-j*w*Rn/c);%声压聚焦方向矢量M*1Pmvdr(k1,k2)=1/abs(b'*RP*b);endend% 归一化for k1=1:length(z)pp1(k1)=max(Pmvdr(k1,:));% Pmvdr的第k1行的最大元素的值endPmvdr=Pmvdr/max(pp1);% Pmvdr所有元素除以其最大值，归一化幅度figure;mesh(x,z,Pmvdr);xlabel('x(m)'),ylabel('z(m)')title('mvdr单亮点声图')figure;pcolor(x,z,Pmvdr);shading interp;xlabel('x(m)'),ylabel('z(m)')title('mvdr单亮点声图')colorbar[a,b] = find(Pmvdr==max(max(Pmvdr))); xm1=b/length(x)*20, zm1=a/length(z)*20, [a b]=find(Pmvdr > 0.5*max(max(Pmvdr))); deta_x=(max(b)-min(b))/length(x)*20 deta_z=(max(a)-min(a))/length(z)*20。
常规波束形成matlab
常规的波束形成方法通常使用MATLAB实现。
以下是一个简单的MATLAB代码示例，演示了如何使用常规的波束形成方法：```matlab设置参数c = 343; 光速，单位为m/sf = 1000; 信号频率，单位为HzN = 128; 阵元数目d = c/f; 阵元间距，单位为mtheta = 1:1:180; 方向角，单位为度t = 0:1/f:1-1/f; 时间，单位为sA = zeros(1,length(theta)); 信号幅度A(5) = 5; 在第5个方向上有目标W = zeros(N,length(theta)); 阵列流形计算阵列流形for i = 1:1:Nfor j = 1:length(theta)W(i,j) = exp(-2*pi*f*(i-1)*d*sind(theta(j))/c);endend计算波束形成输出S = zeros(length(theta),length(t));for j = 1:length(theta)S(j,:) = sin(2*pi*f.*t).*A(j);endS = W*S;画图figure;plot(theta,abs(S));xlabel('方向角（度）');ylabel('幅度');title('波束形成输出');```这段代码使用了常规的波束形成方法，首先设定了光速、信号频率、阵元数目、阵元间距、方向角和时间等参数。
然后，根据给定的信号幅度和阵列流形计算波束形成输出。
最后，使用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=45/180*pi;%来波方向(我觉得应该是天线阵的指向)%theta0=0;%来波方向w=exp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]');for j=1:length(theta) %(我认为是入射角度，即来波方向,计算阵列流形矩阵A) a=exp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]');p(j)=w'*a; %(matlab中的'默认为共轭转置，如果要计算转置为w.'*a)endfigure;plot(theta,abs(p)),grid onxlabel('theta/radian')ylabel('amplitude')title('8阵元均匀线阵方向图')见张小飞的书《阵列信号处理的理论和应用2.3.4节阵列的方向图》当来波方向为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 λ> 时，会出现栅瓣，导致空间模糊。
CAPON波束形成-Matlab程序教学文稿
C A P O N波束形成-M a t l a b程序CAPON波束形成器仿真1．实验原理波束形成就是从传感器阵列重构源信号。
（1）、通过增加期望信源的贡献来实现；（2）、通过抑制掉干扰源来实现。
经典的波束形成需要观测方向(期望信源的方向)的知识。
盲波束形成试图在没有期望信源方向信息的情况下进行信源的恢复。
波束形成技术的基本思想是：通过将各阵元输出进行加权求和，在一时间内将天线阵列波束“导向”到一个方向上，对期望信号得到最大输出功率的导向位置即给出DOA估计。
虽然阵列天线的方向图是全方向的，但阵列的输出经过加权求和后，却可以被调整到阵列接收的方向增益聚集在一个方向上，相当于形成了一个”波束”。
这就是波束形成的物理意义所在。
在智能天线中，波束形成是关键技术之一，是提高信噪比、增加用户容量的保证，能够成倍地提高通信系统的容量，有效地抑制各种干扰，并改善通信质量。
波束形成器的最佳权向量w取决于阵列方向向量)(ak，而在移动通信里用户的方向向量一般未知，需要估计(称之为DOA估计)。
因此，在计算波束形成的最佳权向量之前，必须在已知阵列几何结构的前提下先估计期望信号的波达方向。
Capon 波束形成器求解的优化问题可表述为w arg min P(w)θ=其约束条件为 1)(a w H =θCapon 波束形成器在使噪声和干扰所贡献的功率为最小的同时，保持了期 望信号的功率不变。
因此，它可以看作是一个尖锐的空间带通滤波器。
最优加 权向量w 可以利用Lagrange 乘子法求解，其结果为 )(a R ˆ)(a )(a R ˆw 1H 1CAP θθθ--= 当μ不取常数，而取作)(a R ˆ)(a 11H θθμ-=时，最佳权向量就转变成Capon 波束形成器的权向量。
空间谱为)(a R ˆ)(a 1)(P 1-H CAP θθθ=2．变量定义 M ：均匀线阵列数目P ：信号源个数nn ：快拍数angle1、angle2、angle3：信号来波角度u ：复高斯噪声Ps ：信号能量refp ：信噪比（实值）X：接收信号Rxx：接收信号的相关矩阵doa：波达方向估计3．仿真结果采用上述算法进行仿真，结果如图所示。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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. 当阵元间距 时，会出现栅瓣，导致空间模糊。
仿真图如下：4. 类似于时域滤波，天线方向图是最优权的傅立叶变换 仿真程序和仿真图如下： clcclear all close 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]');/2dλ>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 on xlabel('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;时仿真图如下：。