简介这份MATLAB波束赋形Beamforming仿真代码文档面向无线通信、雷达、声纳等信号处理方向的学习者和工程技术人员可作为课堂教学、课程设计或工程预研的参考资料。文档以doc格式完整收录了均匀线阵方向图绘制、16/128/1024元阵元数对波束宽度和分辨率的影响、阵元间距大于半波长引发栅瓣等仿真代码并给出了按定义方向图与最优权FFT结果的对比以及最大信噪比准则方向图、功率谱和ASC旁瓣相消MSE准则等示例可直观比较不同天线数量和多类波束成形算法的效果。资源包共1个doc文件大小1.23MB已有335人浏览学习。读者既能直接复制运行各小节MATLAB脚本并观察仿真曲线也能借助分节注释梳理阵列流型、加权系数和干扰抑制等关键知识点是快速上手波束赋形实验的实用参考。1. 波束形成的起点把天线阵列搬进 MATLAB雷达在 0 度方向盯住目标20 度方向却有一部干扰机在持续压制常规接收机只能把增益做高而波束形成要做的是给 8 个阵元各分配一个复数权系数让期望方向同相叠加、干扰方向相位相消。这份 MATLAB 资源把这条链路完整串了起来从均匀线阵方向图、阵元数与波束宽度的定量关系到最大信噪比max-SNR、MSE 旁瓣相消、LCMV 和 CaponMVDR这些自适应准则最后落到协方差矩阵不同估计方法对方向图的实际影响。适合正在做阵列信号处理课程设计、雷达或 5G 物理层仿真以及想把手写方向图程序升级成自适应波束形成器的人。整套代码是纯 MATLAB 脚本不依赖额外工具箱拿过来改参数就能跑。2. 均匀线阵方向图阵列流形、半波长间距与波束宽度2.1 导向矢量与阵列流形均匀线阵ULA里最核心的概念是导向矢量。假设 N 个阵元沿直线等间距排列间距为 d一个从 θ 方向来的窄带平面波到达第 n 个阵元时相对参考阵元的相位延迟是 2πnd·sinθ/λ。把这 N 个相位延迟写成列向量就得到该方向上的导向矢量a exp(1i*2*pi*d_lamda*sin(theta)*[0:element_num-1]);其中d_lamda d/λ是归一化阵元间距theta是来波方向[0:element_num-1]生成从 0 到 N-1 的阵元索引。这个向量描述了阵列对某个方向入射信号的幅度和相位响应而方向图本质上就是扫描所有 θ 时加权矢量 w 与导向矢量 a(θ) 内积的模值。d_lamda1/2是均匀线阵里最常出现的配置。间距取半波长有两个原因一是保证整个可见区-90 度到 90 度内不出现栅瓣二是让阵列孔径在物理尺寸受限时尽量大从而获得更窄的主瓣。间距小于半波长会牺牲分辨率大于半波长则会引入空间模糊后面第 2.4 节专门说这个问题。原代码里用了imagsqrt(-1)这种老式复数写法在新版 MATLAB 里imag是内置函数名建议统一换成1i否则跑完这段再调用imag()会出错。2.2 方向图仿真8 阵元均匀线阵先看资源里第一个完整的可运行程序它生成 8 阵元、来波方向 0 度时的均匀线阵方向图% 8 阵元均匀线阵方向图来波方向 0 度 clc; clear all; close all; imagsqrt(-1); % 复数单位新版建议用 1i element_num8; % 阵元数 d_lamda1/2; % 阵元间距 d / 波长 lamda thetalinspace(-pi/2,pi/2,200); theta00; % 来波方向 wexp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]); for j1:length(theta) aexp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]); p(j)w*a; end figure; plot(theta,abs(p)), grid on xlabel(theta/radian) ylabel(amplitude) title(8 阵元均匀线阵方向图)w是对准 theta0 方向的导向矢量a是扫描方向的导向矢量p(j)w*a是两者的内积。因为w和a都是列向量w是共轭转置内积结果是一个复数取abs后就是该扫描角上的幅度响应。这里严格说是“常规波束形成器CBF”的方向图也就是均匀加权情况下的阵列响应主瓣在 0 度第一旁瓣大约 -13.3 dB这是均匀线阵均匀加权的经典理论值跑完代码看 y 轴幅度比就能对上。2.3 波束宽度与阵元数的关系资源里第二段代码用 16、128、1024 三个阵元数做对比核心公式藏在asin(sin(theta) - λ/(N·d))这个反解里element_num116; % 阵元数对比组 element_num2128; element_num31024; lamda0.03; % 波长 0.03 米对应 10 GHz d1/2*lamda; % 半波长间距 theta0:0.5:90; for j1: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)); endfai、psi、beta的实际含义是当波束指向 theta 时由于阵元数有限导致的主瓣偏移量单位是弧度。lamda/(element_num*d)这一项就是波束宽度的尺度因子N 越大这个值越小主瓣越窄。在 theta0 处16 阵元的偏移约 7.2 度128 阵元降到不足 1 度1024 阵元则接近 0.1 度量级。把三条曲线画在一张图里能直观看到波束宽度随阵元数增加的收紧过程。实际系统里阵元数决定阵列孔径孔径越大角分辨力越高但成本和通道数也线性上涨这是个基本权衡。2.4 间距超过半波长时栅瓣与空间模糊资源里对栅瓣只给了现象描述“当阵元间距 d λ/2 时会出现栅瓣导致空间模糊”。要复现这个现象改一行参数就行把第 2.2 节代码里的d_lamda从1/2改成0.75或1方向图会在可见区出现第二个与主瓣等高的峰。栅瓣的物理原因是相位差 2πd·sinθ/λ 在 θ 从 -90 度扫到 90 度时变化范围超过 2π阵列无法区分两个不同方向的来波产生空间模糊。这也是为什么均匀线阵的半波长间距是“最大无栅瓣间距”而不是经验值。做稀布阵或稀疏阵设计时栅瓣抑制是排布优化的重要约束这里的原理同样适用。3. 方向图是权矢量的傅里叶变换从常规波束到最大信噪比准则3.1 为什么方向图能跟 FFT 对应对均匀线阵导向矢量的第 n 个元素是 exp(j·2π·d_lamda·sinθ·n)。把 sinθ 看成空间频率轴的变量导向矢量就是一组空间复指数权矢量 w 与 a(θ) 的内积本质上就是 w 的空间傅里叶变换在某个 sinθ 处的值。换句话说方向图不是“算出来的曲线”而是权矢量的频谱。这个视角的意义在于设计方向图可以等价为设计一个空间滤波器时域 FIR 滤波器的窗函数设计、频域采样设计方法在这里都能平移过来用。这段代码里的fft(w,128)就是直接对权矢量做傅里叶变换得到方向图。3.2 定义法与 FFT 法的结果对比element_num32; % 阵元数 source_num1; % 信源数 d_lamda1/2; thetalinspace(-pi/2,pi/2,200); theta00; % 来波方向 wexp(ima*2*pi*d_lamda*sin(theta0)*[0:element_num-1]); for j1:length(theta) aexp(ima*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]); p(j)w*a; end subplot(1,2,1) plot(theta,abs(p)), grid on % 按定义扫描计算的方向图 pfftfftshift(fft(w,128)); subplot(1,2,2) plot(linspace(-pi/2,pi/2,128),abs(pfft)), grid onfft(w,128)先把 32 点权矢量补零到 128 点再做 FFT补零相当于频域插值让方向图曲线更平滑。fftshift把零频分量移到数组中间对应把 sinθ0也就是 0 度方向放到横轴中心。两幅图主瓣位置和形状一致差别只在横轴映射方式左边是 sinθ 的等间隔采样右边是 FFT 频点的等间隔采样二者在空间频率域是同一回事。这里有个经常被忽视的细节FFT 法的横轴是linspace(-pi/2,pi/2,128)它把 128 个频点线性映射到 -90 度到 90 度。这个映射只对 sinθ 是线性的对 θ 本身不是。所以在大角度区域FFT 画出来的方向图和按定义扫描的结果会有轻微形变角度越大差异越明显。3.3 最大信噪比准则广义特征分解最大信噪比准则的目标是让输出信噪比最大化。设阵列输出 y wx其中 x 包含期望信号 s、干扰 j 和噪声 n则输出信噪比可以写成 Rayleigh 商形式SINR (wRs w) / (wRnj w)其中 Rs 是信号自相关矩阵Rnj 是干扰加噪声自相关矩阵。最大化这个比值等价于求解广义特征分解Rs v λ Rnj v取最大广义特征值对应的特征向量。element_num8; % 阵元数 d_lamda1/2; theta-90:0.5:90; theta00; % 来波方向 theta120; % 干扰方向 L512; % 采样单元数 for i1:L amp010*randn(1); % 信号幅度随机保证各快拍不相关 amp1200*randn(1); % 干扰幅度功率远大于信号 ampn1; % 噪声幅度 s(:,i)amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]); j(:,i)amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]); n(:,i)ampn*(randn(element_num,1)ima*randn(element_num,1)); end Rs1/L*s*s; % 信号自相关矩阵 Rnj1/L*(j*jn*n); % 干扰噪声自相关矩阵 [V,D]eig(Rs,Rnj); % 广义特征分解 [D,I]sort(diag(D)); % 特征值升序排列 WoptV(:,I(8)); % 取最大特征值对应的特征向量amp1200*randn(1)意味着干扰幅度标准差是信号的 20 倍干扰功率约为信号的 400 倍这是故意设置的高干信比场景。eig(Rs,Rnj)做的是广义特征分解不是普通特征分解求解的是 Rs·v λ·Rnj·v。sort(diag(D))默认升序取I(8)就是最后一个、也就是最大的特征值对应的列它在数值上等价于I(end)。这个权矢量会自动在 20 度干扰方向形成深零陷同时保持 0 度方向增益效果比固定波束形成好得多。3.4 幅度配比如何影响方向图原代码里amp0、amp1、amp2的配比直接决定方向图质量。干扰幅度远超信号时max-SNR 权会把主要自由度用于抑制干扰零陷很深反过来如果信号幅度远大于干扰算法会倾向于保持主瓣增益零陷相对变浅。这是自适应波束形成的固有行为——自由度有限时权矢量会优先抑制功率大的分量。做仿真时调整这几个幅度参数比改任何其他参数都更能直观理解自适应算法的取舍逻辑。4. 自适应波束形成的三条主线MSE 旁瓣相消、LCMV 与 Capon4.1 三种准则的选型对比max-SNR 在第 3 章已经讲过它需要分别估计 Rs 和 Rnj实际系统里往往分不开。工程上更常用的是下面这三种做法准则优化目标需要的先验信息典型实现MSE最小均方误差输出与参考信号误差最小参考信号或期望信号维纳解R⁻¹r_xdLCMV线性约束最小方差约束方向增益输出功率最小约束矩阵 C 和响应向量 F闭式解R⁻¹C(CR⁻¹C)⁻¹FCapon / MVDR期望方向增益为 1输出功率最小期望方向导向矢量闭式解R⁻¹a/(aR⁻¹a)MSE 是“有参考信号”的准则LCMV 是“多点约束”的准则Capon 是“单点约束”的 MVDR 特例。下面分别拆开看。4.2 MSE 准则旁瓣相消器旁瓣相消SLC的结构是一个主通道用全向天线接收“信号干扰噪声”M 个辅助通道主要接收干扰通过加权辅助通道把主通道里的干扰分量对消掉。M32; % 辅助天线数目 d_lamda.5; theta0-30; % 来波方向 theta160; % 干扰方向 L512; % 采样单元数 for ii1:L amp01*randn(1); amp1200*randn(1); ampn1; jam(:,ii)amp1*exp(ima*2*pi*0.5*sin(theta1*pi/180)*[0:M-1]) ... ampn*(randn(M,1)ima*randn(M,1)); % 辅助通道干扰噪声 s(ii)amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180)) ... amp1*exp(ima*2*pi*0.5*sin(theta1*pi/180)) ... ampn*(randn(1,1)ima*randn(1,1)); % 主通道信号干扰噪声 s0(ii)amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180)); % 纯信号分量 end Rx1/L*jam*jam; % 辅助通道自相关矩阵 r_xd1/L*jam*s; % 辅助通道与主通道的互相关向量 Woptpinv(Rx)*r_xd; % 维纳解用伪逆避免奇异Rx是辅助通道的协方差矩阵r_xd是辅助通道与主通道的互相关。维纳解 Wopt pinv(Rx)·r_xd 的意义是找到一组辅助通道权让加权后的辅助输出尽量逼近主通道里的干扰分量然后从主通道减去。这里用pinv而不是inv很关键M32、L512 时协方差矩阵可能接近奇异伪逆能给出最小范数解避免数值爆炸。代码最后还计算了delta s0 - (s - Wopt*jam)的方差这是衡量对消残差的指标如果对消理想delta应该逼近纯信号 s0方差接近于噪声功率。这个验证思路可以沿用到任何自适应算法里不只是 SLC。4.3 LCMV 准则多点约束的闭式解LCMV 的思路是对某些方向施加硬约束增益固定为指定值在满足约束的前提下最小化输出功率。约束写成矩阵形式 Cw F其中 C 的每一列是一个约束方向的导向矢量F 是对应方向的期望响应。element_num8; theta00; theta130; theta260; % ... 信号生成部分与前面相同三个方向的幅度分别为 amp0/amp1/amp2 ... C[steer1 steer2 steer3]; % 约束矩阵0/30/60 度导向矢量 F[1 0 1]; % 期望响应0 度和 60 度增益为 130 度为 0 winv(Rx)*C*(inv(C*inv(Rx)*C))*F;F[1 0 1]表示把 0 度和 60 度都当作“要保留的方向”30 度方向约束为零增益。这个配置下方向图会在 0 度和 60 度保持响应在 30 度形成零陷。如果改成F[1 0 0]那就是大作业里的经典做法只保留 0 度来波方向45 度和 60 度两个干扰方向全部约束为零方向图会同时在这两个角度出现深零陷。注意theta-90:0.5:90-0.3这种写法终点取 89.7 而不是 90是为了让扫描角度避开边界防止在 ±90 度附近出现数值异常。扫描角度与约束角度重合时LCMV 方向图会在该点出现尖峰或凹陷这是约束点与扫描点重合后的正常现象。4.4 Capon / MVDR单点约束的最小方差Capon 波束形成是 LCMV 在单点约束下的特例只约束期望方向增益为 1其余全部交给“最小化输出功率”去处理。Rx1/L*x*x; % 接收信号自相关矩阵 Rinv(Rx); steerexp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]); wR*steer/(steer*R*steer); % Capon 最优权矢量 for j1:length(theta) aexp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]); f(j)w*a; % 方向图 p(j)1/(a*R*a); % Capon 功率谱 end闭式解里的steer*R*steer是标量整个式子保证steer*w 1。方向图由f(j)w*a给出而p(j)1/(a*R*a)是 Capon 谱——它估计的是“从 θ 方向入射的信号功率”而不是阵列输出功率。Capon 谱比常规周期图有更高的角度分辨力代价是每次谱估计都要做一次矩阵求逆。原代码里干扰方向取 20 度和 60 度信号幅度 10、干扰幅度 200方向图上两个干扰位置的谱峰会被自适应权压下去而 Capon 谱则在这两个角度显示干扰功率的估计值。5. 协方差矩阵的工程细节多点约束、快拍数与幅度配比5.1 用 Rx 还是用 Rnj结果差很多资源第 9 段代码做了件很有工程意义的事分别用接收信号自相关矩阵 Rx 和干扰加噪声自相关矩阵 Rnj 计算 Capon 权对比方向图差异。Wopt_Rxinv(Rx)*e/(e*inv(Rx)*e); % 用接收信号自相关矩阵 Wopt_Rnjinv(Rnj)*e/(e*inv(Rnj)*e); % 用干扰噪声自相关矩阵用 Rx 时期望信号分量也被算进了协方差矩阵Capon 会把信号本身当作“要抑制的成分”导致期望方向的方向图出现凹陷信噪比越高凹陷越明显。用 Rnj 时协方差矩阵里只有干扰和噪声期望方向的增益被完整保留。实际系统里拿不到纯 Rnj常见的近似做法是先在没有期望信号的时候估计 Rnj或者在协方差矩阵里对角加载牺牲一点信噪比换方向图保真。5.2 快拍数 L 与对角加载大作业部分对比了 L512 和 L2048 两组参数规律很明确快拍数翻四倍协方差矩阵估计方差减半零陷更深、位置更准、旁瓣更干净。快拍数太少时Rx 的特征值分散小特征值对应的方向会出现随机高旁瓣LCMV 的求逆过程也会放大噪声。工程上我一般会加一层保护Rx_loaded Rx 1e-6 * eye(element_num); % 对角加载 w inv(Rx_loaded) * C * (inv(C*inv(Rx_loaded)*C)) * F;对角加载相当于给协方差矩阵的所有特征值加了一个底噪把病态条件数拉下来代价是零陷深度略微变浅。加载量取1e-3到1e-6倍的主对角线均值具体看干扰强度和数值精度。5.3 直接叠加导向矢量做权是固定波束而非自适应资源第 10 段代码里有个值得警惕的写法w amp0*steer0 amp2*steer2 amp1*steer1直接把三个方向的导向矢量按幅度叠加当成权。这样得到的方向图只是三个常规波束的叠加不会像 LCMV 那样在约束方向之外的位置自动调零。它是“固定多波束”不是“自适应权”。如果在课程设计里用这个当 Capon 结果交差容易被追问权矢量的推导来源。真正的多点约束 Capon 权应该是第 4.3 节的闭式解而不是导向矢量的线性组合。5.4 批量扫描参数快速验证方向图做参数实验时我一般把信号生成封装成循环批量跑不同配置L_set [128 512 2048]; amp_scale [1 10 100]; for k 1:length(L_set) for m 1:length(amp_scale) L L_set(k); amp1 amp_scale(m)*randn(1); % ... 信号生成与权计算 ... F_dB 20*log10(abs(f)/max(abs(f))); % 检查零陷深度theta1 处是否低于 -30 dB idx find(theta theta1, 1); null_depth(k, m) F_dB(idx); end end验证方向图质量就看三件事主瓣峰值是否落在 theta0、零陷位置是否精确对准 theta1 和 theta2、方向图是否关于 0 度大致对称。零陷深度低于 -30 dB 算合格-50 dB 属于干扰抑制很充分的状态一般受限于快拍数和数值精度。把这组检查写进脚本参数调起来就不用肉眼看图了。本文还有配套的精品资源点击获取
MATLAB波束形成仿真：从均匀线阵方向图到自适应波束形成
简介这份MATLAB波束赋形Beamforming仿真代码文档面向无线通信、雷达、声纳等信号处理方向的学习者和工程技术人员可作为课堂教学、课程设计或工程预研的参考资料。文档以doc格式完整收录了均匀线阵方向图绘制、16/128/1024元阵元数对波束宽度和分辨率的影响、阵元间距大于半波长引发栅瓣等仿真代码并给出了按定义方向图与最优权FFT结果的对比以及最大信噪比准则方向图、功率谱和ASC旁瓣相消MSE准则等示例可直观比较不同天线数量和多类波束成形算法的效果。资源包共1个doc文件大小1.23MB已有335人浏览学习。读者既能直接复制运行各小节MATLAB脚本并观察仿真曲线也能借助分节注释梳理阵列流型、加权系数和干扰抑制等关键知识点是快速上手波束赋形实验的实用参考。1. 波束形成的起点把天线阵列搬进 MATLAB雷达在 0 度方向盯住目标20 度方向却有一部干扰机在持续压制常规接收机只能把增益做高而波束形成要做的是给 8 个阵元各分配一个复数权系数让期望方向同相叠加、干扰方向相位相消。这份 MATLAB 资源把这条链路完整串了起来从均匀线阵方向图、阵元数与波束宽度的定量关系到最大信噪比max-SNR、MSE 旁瓣相消、LCMV 和 CaponMVDR这些自适应准则最后落到协方差矩阵不同估计方法对方向图的实际影响。适合正在做阵列信号处理课程设计、雷达或 5G 物理层仿真以及想把手写方向图程序升级成自适应波束形成器的人。整套代码是纯 MATLAB 脚本不依赖额外工具箱拿过来改参数就能跑。2. 均匀线阵方向图阵列流形、半波长间距与波束宽度2.1 导向矢量与阵列流形均匀线阵ULA里最核心的概念是导向矢量。假设 N 个阵元沿直线等间距排列间距为 d一个从 θ 方向来的窄带平面波到达第 n 个阵元时相对参考阵元的相位延迟是 2πnd·sinθ/λ。把这 N 个相位延迟写成列向量就得到该方向上的导向矢量a exp(1i*2*pi*d_lamda*sin(theta)*[0:element_num-1]);其中d_lamda d/λ是归一化阵元间距theta是来波方向[0:element_num-1]生成从 0 到 N-1 的阵元索引。这个向量描述了阵列对某个方向入射信号的幅度和相位响应而方向图本质上就是扫描所有 θ 时加权矢量 w 与导向矢量 a(θ) 内积的模值。d_lamda1/2是均匀线阵里最常出现的配置。间距取半波长有两个原因一是保证整个可见区-90 度到 90 度内不出现栅瓣二是让阵列孔径在物理尺寸受限时尽量大从而获得更窄的主瓣。间距小于半波长会牺牲分辨率大于半波长则会引入空间模糊后面第 2.4 节专门说这个问题。原代码里用了imagsqrt(-1)这种老式复数写法在新版 MATLAB 里imag是内置函数名建议统一换成1i否则跑完这段再调用imag()会出错。2.2 方向图仿真8 阵元均匀线阵先看资源里第一个完整的可运行程序它生成 8 阵元、来波方向 0 度时的均匀线阵方向图% 8 阵元均匀线阵方向图来波方向 0 度 clc; clear all; close all; imagsqrt(-1); % 复数单位新版建议用 1i element_num8; % 阵元数 d_lamda1/2; % 阵元间距 d / 波长 lamda thetalinspace(-pi/2,pi/2,200); theta00; % 来波方向 wexp(imag*2*pi*d_lamda*sin(theta0)*[0:element_num-1]); for j1:length(theta) aexp(imag*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]); p(j)w*a; end figure; plot(theta,abs(p)), grid on xlabel(theta/radian) ylabel(amplitude) title(8 阵元均匀线阵方向图)w是对准 theta0 方向的导向矢量a是扫描方向的导向矢量p(j)w*a是两者的内积。因为w和a都是列向量w是共轭转置内积结果是一个复数取abs后就是该扫描角上的幅度响应。这里严格说是“常规波束形成器CBF”的方向图也就是均匀加权情况下的阵列响应主瓣在 0 度第一旁瓣大约 -13.3 dB这是均匀线阵均匀加权的经典理论值跑完代码看 y 轴幅度比就能对上。2.3 波束宽度与阵元数的关系资源里第二段代码用 16、128、1024 三个阵元数做对比核心公式藏在asin(sin(theta) - λ/(N·d))这个反解里element_num116; % 阵元数对比组 element_num2128; element_num31024; lamda0.03; % 波长 0.03 米对应 10 GHz d1/2*lamda; % 半波长间距 theta0:0.5:90; for j1: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)); endfai、psi、beta的实际含义是当波束指向 theta 时由于阵元数有限导致的主瓣偏移量单位是弧度。lamda/(element_num*d)这一项就是波束宽度的尺度因子N 越大这个值越小主瓣越窄。在 theta0 处16 阵元的偏移约 7.2 度128 阵元降到不足 1 度1024 阵元则接近 0.1 度量级。把三条曲线画在一张图里能直观看到波束宽度随阵元数增加的收紧过程。实际系统里阵元数决定阵列孔径孔径越大角分辨力越高但成本和通道数也线性上涨这是个基本权衡。2.4 间距超过半波长时栅瓣与空间模糊资源里对栅瓣只给了现象描述“当阵元间距 d λ/2 时会出现栅瓣导致空间模糊”。要复现这个现象改一行参数就行把第 2.2 节代码里的d_lamda从1/2改成0.75或1方向图会在可见区出现第二个与主瓣等高的峰。栅瓣的物理原因是相位差 2πd·sinθ/λ 在 θ 从 -90 度扫到 90 度时变化范围超过 2π阵列无法区分两个不同方向的来波产生空间模糊。这也是为什么均匀线阵的半波长间距是“最大无栅瓣间距”而不是经验值。做稀布阵或稀疏阵设计时栅瓣抑制是排布优化的重要约束这里的原理同样适用。3. 方向图是权矢量的傅里叶变换从常规波束到最大信噪比准则3.1 为什么方向图能跟 FFT 对应对均匀线阵导向矢量的第 n 个元素是 exp(j·2π·d_lamda·sinθ·n)。把 sinθ 看成空间频率轴的变量导向矢量就是一组空间复指数权矢量 w 与 a(θ) 的内积本质上就是 w 的空间傅里叶变换在某个 sinθ 处的值。换句话说方向图不是“算出来的曲线”而是权矢量的频谱。这个视角的意义在于设计方向图可以等价为设计一个空间滤波器时域 FIR 滤波器的窗函数设计、频域采样设计方法在这里都能平移过来用。这段代码里的fft(w,128)就是直接对权矢量做傅里叶变换得到方向图。3.2 定义法与 FFT 法的结果对比element_num32; % 阵元数 source_num1; % 信源数 d_lamda1/2; thetalinspace(-pi/2,pi/2,200); theta00; % 来波方向 wexp(ima*2*pi*d_lamda*sin(theta0)*[0:element_num-1]); for j1:length(theta) aexp(ima*2*pi*d_lamda*sin(theta(j))*[0:element_num-1]); p(j)w*a; end subplot(1,2,1) plot(theta,abs(p)), grid on % 按定义扫描计算的方向图 pfftfftshift(fft(w,128)); subplot(1,2,2) plot(linspace(-pi/2,pi/2,128),abs(pfft)), grid onfft(w,128)先把 32 点权矢量补零到 128 点再做 FFT补零相当于频域插值让方向图曲线更平滑。fftshift把零频分量移到数组中间对应把 sinθ0也就是 0 度方向放到横轴中心。两幅图主瓣位置和形状一致差别只在横轴映射方式左边是 sinθ 的等间隔采样右边是 FFT 频点的等间隔采样二者在空间频率域是同一回事。这里有个经常被忽视的细节FFT 法的横轴是linspace(-pi/2,pi/2,128)它把 128 个频点线性映射到 -90 度到 90 度。这个映射只对 sinθ 是线性的对 θ 本身不是。所以在大角度区域FFT 画出来的方向图和按定义扫描的结果会有轻微形变角度越大差异越明显。3.3 最大信噪比准则广义特征分解最大信噪比准则的目标是让输出信噪比最大化。设阵列输出 y wx其中 x 包含期望信号 s、干扰 j 和噪声 n则输出信噪比可以写成 Rayleigh 商形式SINR (wRs w) / (wRnj w)其中 Rs 是信号自相关矩阵Rnj 是干扰加噪声自相关矩阵。最大化这个比值等价于求解广义特征分解Rs v λ Rnj v取最大广义特征值对应的特征向量。element_num8; % 阵元数 d_lamda1/2; theta-90:0.5:90; theta00; % 来波方向 theta120; % 干扰方向 L512; % 采样单元数 for i1:L amp010*randn(1); % 信号幅度随机保证各快拍不相关 amp1200*randn(1); % 干扰幅度功率远大于信号 ampn1; % 噪声幅度 s(:,i)amp0*exp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]); j(:,i)amp1*exp(ima*2*pi*1/2*sin(theta1*pi/180)*[0:element_num-1]); n(:,i)ampn*(randn(element_num,1)ima*randn(element_num,1)); end Rs1/L*s*s; % 信号自相关矩阵 Rnj1/L*(j*jn*n); % 干扰噪声自相关矩阵 [V,D]eig(Rs,Rnj); % 广义特征分解 [D,I]sort(diag(D)); % 特征值升序排列 WoptV(:,I(8)); % 取最大特征值对应的特征向量amp1200*randn(1)意味着干扰幅度标准差是信号的 20 倍干扰功率约为信号的 400 倍这是故意设置的高干信比场景。eig(Rs,Rnj)做的是广义特征分解不是普通特征分解求解的是 Rs·v λ·Rnj·v。sort(diag(D))默认升序取I(8)就是最后一个、也就是最大的特征值对应的列它在数值上等价于I(end)。这个权矢量会自动在 20 度干扰方向形成深零陷同时保持 0 度方向增益效果比固定波束形成好得多。3.4 幅度配比如何影响方向图原代码里amp0、amp1、amp2的配比直接决定方向图质量。干扰幅度远超信号时max-SNR 权会把主要自由度用于抑制干扰零陷很深反过来如果信号幅度远大于干扰算法会倾向于保持主瓣增益零陷相对变浅。这是自适应波束形成的固有行为——自由度有限时权矢量会优先抑制功率大的分量。做仿真时调整这几个幅度参数比改任何其他参数都更能直观理解自适应算法的取舍逻辑。4. 自适应波束形成的三条主线MSE 旁瓣相消、LCMV 与 Capon4.1 三种准则的选型对比max-SNR 在第 3 章已经讲过它需要分别估计 Rs 和 Rnj实际系统里往往分不开。工程上更常用的是下面这三种做法准则优化目标需要的先验信息典型实现MSE最小均方误差输出与参考信号误差最小参考信号或期望信号维纳解R⁻¹r_xdLCMV线性约束最小方差约束方向增益输出功率最小约束矩阵 C 和响应向量 F闭式解R⁻¹C(CR⁻¹C)⁻¹FCapon / MVDR期望方向增益为 1输出功率最小期望方向导向矢量闭式解R⁻¹a/(aR⁻¹a)MSE 是“有参考信号”的准则LCMV 是“多点约束”的准则Capon 是“单点约束”的 MVDR 特例。下面分别拆开看。4.2 MSE 准则旁瓣相消器旁瓣相消SLC的结构是一个主通道用全向天线接收“信号干扰噪声”M 个辅助通道主要接收干扰通过加权辅助通道把主通道里的干扰分量对消掉。M32; % 辅助天线数目 d_lamda.5; theta0-30; % 来波方向 theta160; % 干扰方向 L512; % 采样单元数 for ii1:L amp01*randn(1); amp1200*randn(1); ampn1; jam(:,ii)amp1*exp(ima*2*pi*0.5*sin(theta1*pi/180)*[0:M-1]) ... ampn*(randn(M,1)ima*randn(M,1)); % 辅助通道干扰噪声 s(ii)amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180)) ... amp1*exp(ima*2*pi*0.5*sin(theta1*pi/180)) ... ampn*(randn(1,1)ima*randn(1,1)); % 主通道信号干扰噪声 s0(ii)amp0*exp(ima*2*pi*0.5*sin(theta0*pi/180)); % 纯信号分量 end Rx1/L*jam*jam; % 辅助通道自相关矩阵 r_xd1/L*jam*s; % 辅助通道与主通道的互相关向量 Woptpinv(Rx)*r_xd; % 维纳解用伪逆避免奇异Rx是辅助通道的协方差矩阵r_xd是辅助通道与主通道的互相关。维纳解 Wopt pinv(Rx)·r_xd 的意义是找到一组辅助通道权让加权后的辅助输出尽量逼近主通道里的干扰分量然后从主通道减去。这里用pinv而不是inv很关键M32、L512 时协方差矩阵可能接近奇异伪逆能给出最小范数解避免数值爆炸。代码最后还计算了delta s0 - (s - Wopt*jam)的方差这是衡量对消残差的指标如果对消理想delta应该逼近纯信号 s0方差接近于噪声功率。这个验证思路可以沿用到任何自适应算法里不只是 SLC。4.3 LCMV 准则多点约束的闭式解LCMV 的思路是对某些方向施加硬约束增益固定为指定值在满足约束的前提下最小化输出功率。约束写成矩阵形式 Cw F其中 C 的每一列是一个约束方向的导向矢量F 是对应方向的期望响应。element_num8; theta00; theta130; theta260; % ... 信号生成部分与前面相同三个方向的幅度分别为 amp0/amp1/amp2 ... C[steer1 steer2 steer3]; % 约束矩阵0/30/60 度导向矢量 F[1 0 1]; % 期望响应0 度和 60 度增益为 130 度为 0 winv(Rx)*C*(inv(C*inv(Rx)*C))*F;F[1 0 1]表示把 0 度和 60 度都当作“要保留的方向”30 度方向约束为零增益。这个配置下方向图会在 0 度和 60 度保持响应在 30 度形成零陷。如果改成F[1 0 0]那就是大作业里的经典做法只保留 0 度来波方向45 度和 60 度两个干扰方向全部约束为零方向图会同时在这两个角度出现深零陷。注意theta-90:0.5:90-0.3这种写法终点取 89.7 而不是 90是为了让扫描角度避开边界防止在 ±90 度附近出现数值异常。扫描角度与约束角度重合时LCMV 方向图会在该点出现尖峰或凹陷这是约束点与扫描点重合后的正常现象。4.4 Capon / MVDR单点约束的最小方差Capon 波束形成是 LCMV 在单点约束下的特例只约束期望方向增益为 1其余全部交给“最小化输出功率”去处理。Rx1/L*x*x; % 接收信号自相关矩阵 Rinv(Rx); steerexp(ima*2*pi*1/2*sin(theta0*pi/180)*[0:element_num-1]); wR*steer/(steer*R*steer); % Capon 最优权矢量 for j1:length(theta) aexp(ima*2*pi*d_lamda*sin(theta(j)*pi/180)*[0:element_num-1]); f(j)w*a; % 方向图 p(j)1/(a*R*a); % Capon 功率谱 end闭式解里的steer*R*steer是标量整个式子保证steer*w 1。方向图由f(j)w*a给出而p(j)1/(a*R*a)是 Capon 谱——它估计的是“从 θ 方向入射的信号功率”而不是阵列输出功率。Capon 谱比常规周期图有更高的角度分辨力代价是每次谱估计都要做一次矩阵求逆。原代码里干扰方向取 20 度和 60 度信号幅度 10、干扰幅度 200方向图上两个干扰位置的谱峰会被自适应权压下去而 Capon 谱则在这两个角度显示干扰功率的估计值。5. 协方差矩阵的工程细节多点约束、快拍数与幅度配比5.1 用 Rx 还是用 Rnj结果差很多资源第 9 段代码做了件很有工程意义的事分别用接收信号自相关矩阵 Rx 和干扰加噪声自相关矩阵 Rnj 计算 Capon 权对比方向图差异。Wopt_Rxinv(Rx)*e/(e*inv(Rx)*e); % 用接收信号自相关矩阵 Wopt_Rnjinv(Rnj)*e/(e*inv(Rnj)*e); % 用干扰噪声自相关矩阵用 Rx 时期望信号分量也被算进了协方差矩阵Capon 会把信号本身当作“要抑制的成分”导致期望方向的方向图出现凹陷信噪比越高凹陷越明显。用 Rnj 时协方差矩阵里只有干扰和噪声期望方向的增益被完整保留。实际系统里拿不到纯 Rnj常见的近似做法是先在没有期望信号的时候估计 Rnj或者在协方差矩阵里对角加载牺牲一点信噪比换方向图保真。5.2 快拍数 L 与对角加载大作业部分对比了 L512 和 L2048 两组参数规律很明确快拍数翻四倍协方差矩阵估计方差减半零陷更深、位置更准、旁瓣更干净。快拍数太少时Rx 的特征值分散小特征值对应的方向会出现随机高旁瓣LCMV 的求逆过程也会放大噪声。工程上我一般会加一层保护Rx_loaded Rx 1e-6 * eye(element_num); % 对角加载 w inv(Rx_loaded) * C * (inv(C*inv(Rx_loaded)*C)) * F;对角加载相当于给协方差矩阵的所有特征值加了一个底噪把病态条件数拉下来代价是零陷深度略微变浅。加载量取1e-3到1e-6倍的主对角线均值具体看干扰强度和数值精度。5.3 直接叠加导向矢量做权是固定波束而非自适应资源第 10 段代码里有个值得警惕的写法w amp0*steer0 amp2*steer2 amp1*steer1直接把三个方向的导向矢量按幅度叠加当成权。这样得到的方向图只是三个常规波束的叠加不会像 LCMV 那样在约束方向之外的位置自动调零。它是“固定多波束”不是“自适应权”。如果在课程设计里用这个当 Capon 结果交差容易被追问权矢量的推导来源。真正的多点约束 Capon 权应该是第 4.3 节的闭式解而不是导向矢量的线性组合。5.4 批量扫描参数快速验证方向图做参数实验时我一般把信号生成封装成循环批量跑不同配置L_set [128 512 2048]; amp_scale [1 10 100]; for k 1:length(L_set) for m 1:length(amp_scale) L L_set(k); amp1 amp_scale(m)*randn(1); % ... 信号生成与权计算 ... F_dB 20*log10(abs(f)/max(abs(f))); % 检查零陷深度theta1 处是否低于 -30 dB idx find(theta theta1, 1); null_depth(k, m) F_dB(idx); end end验证方向图质量就看三件事主瓣峰值是否落在 theta0、零陷位置是否精确对准 theta1 和 theta2、方向图是否关于 0 度大致对称。零陷深度低于 -30 dB 算合格-50 dB 属于干扰抑制很充分的状态一般受限于快拍数和数值精度。把这组检查写进脚本参数调起来就不用肉眼看图了。本文还有配套的精品资源点击获取
