简介数字多波束形成DBF是阵列信号处理的核心方向之一广泛用于雷达、通信与声呐系统的空域滤波设计这份仿真代码适合理解多波束扫描与方向图绘制原理的初学者也可供相关课题工程师快速搭建验证环境。资源包共2个文件包含1个MATLAB脚本与1个说明文档压缩包仅2KB脚本是完整可运行的方向图仿真主程序说明文档则梳理了文件构成与使用注意事项结构非常精简便于下载后直接上手。代码支持同时配置四个波束的方位角、俯仰角与载波频率阵元坐标、间距和数量均可自由调整能直观展示波束指向与参数设置一一对应的效果仿真分别绘制方位面与俯仰面方向图参数设置、权值计算、波束扫描等步骤划分清晰关键位置附文字注释方便逐段理解DBF实现流程。目前已有988人学习下载对学习阵列信号处理或准备课程设计的学生可用这份代码替换不同阵元排布与波束参数观察方向图变化加深对数字波束形成原理的理解。1. 数字多波束形成 DBF 仿真在 MATLAB 里到底仿什么做相控阵、雷达或者测向接收机的人早晚会遇到“数字多波束形成”这个需求同一副天线阵列采集到的数据要同时输出多个方向的波束而不是一个一个地扫。MATLAB 的 DBF 仿真代码核心不是把方向图画得漂亮而是把阵列流形、复数权值、快照矩阵这几块数据结构搭对让方向图、多波束输出和测角结果互相验证。本文针对均匀线阵ULA的常规多波束形成给出一套能从零跑通、能改参数、能排查错误的 MATLAB 实现思路。这里先声明一个检索上的坑标题里的 DBF 是 Digital BeamForming数字波束形成跟 Oracle 数据库的 DBF 数据文件没有关系搜资料时别被带偏。2. 单波束复现均匀线阵方向图的最小 MATLAB 代码2.1 阵列流形向量相位差公式到矩阵表达均匀线阵是 DBF 仿真最常用的阵型N 个全向阵元等间距 d 排成直线远场窄带信号以平面波入射。以阵列法线方向为 0°入射角 θ 定义在这个坐标系里第 n 个阵元相对第 0 个阵元的波程差是 n·d·sinθ相位差是 2π·n·d·sinθ/λ。把波长归一化掉令 d_lambda d/λ则阵列流形向量写成a(θ) [1, exp(j2πd_lambda·sinθ), …, exp(j2π(N-1)d_lambda·sinθ)]^T在 MATLAB 里生成流形向量最省事也最不容易错的方式是用“列向量乘行向量”的外积一次性生成扫描角矩阵N 16; % 阵元数 d_lambda 0.5; % 阵元间距/波长半波长是标准配置 theta_deg (-90:0.1:90); % 扫描角范围步进0.1° theta_rad theta_deg * pi / 180; n (0:N-1).; % 阵元索引列向量 N×1 A exp(1j * 2 * pi * d_lambda * n * sin(theta_rad)); % N×length(theta_deg)这里的A就是流形矩阵每一列对应一个扫描角度每一行对应一个阵元。n * sin(theta_rad)利用 MATLAB 的隐式扩展完成外积生成 N×L 的相位指数矩阵。关键不是“会乘”而是明白A(:, k)是什么——它代表第 k 个扫描角度下的理想阵列响应。后续所有方向图、波束权值、数据快照都要和它保持同样的维度约定。新手最容易在这里把 n 写成行向量导致矩阵维度转置后面矩阵乘法全部错乱排查时往往耗费大量时间。要验证流形矩阵是否生成正确可以抽查两个特殊方向。θ0° 时所有阵元同相流形向量元素全是 1θ90° 端射方向时相邻阵元相位差是 jπ符号交替。用两行命令快速自检disp(A(:, find(theta_deg 0, 1))); % 法线方向预期全为1 disp(A(:, find(theta_deg 90, 1))); % 端射方向预期相邻符号交替正常情况输出里第一段是 10i 的向量第二段是 1、-1 交替。和预期不一致基本是角度定义或 sin 参数写错了。这个自检步骤虽然简单但能挡住大量后续返工。2.2 常规波束形成权值计算代码与参数常规波束形成CBF就是按期望方向构造权向量再做共轭加权求和。权向量取指向方向 θ0 的流形向量加权输出为 y(θ) w^H a(θ)其中 w a(θ0)。当入射方向 θ 等于 θ0 时各阵元信号相位对齐、同相叠加输出最大否则相位逐渐失配输出随方向图衰减。MATLAB 实现如下% 波束指向角 theta0_deg 20; theta0_rad theta0_deg * pi / 180; % 指向20°的权值注意这里没取共轭w会做共轭转置 w exp(1j * 2 * pi * d_lambda * n * sin(theta0_rad)); % N×1 % 方向图对每个扫描角做加权求和 pattern w * A; % 1×L 复数 pattern_dB 20 * log10(abs(pattern) / max(abs(pattern))); % 绘图 figure(Color,white); plot(theta_deg, pattern_dB, LineWidth, 1.2); xlabel(入射角 (deg)); ylabel(归一化增益 (dB)); ylim([-60 0]); grid on; title(sprintf(ULA方向图 N%d 指向%g°, N, theta0_deg));代码里有三个点需要说明。第一w是共轭转置Hermitian不是普通转置因为流形向量是复数只有共轭加权才能把相位补偿回来。第二方向图以w * A的复幅度取模、再以最大值归一化后转 dB主瓣顶部永远在 0 dB不同扫描角之间比较相对电平才有意义。第三ylim取 -60 dB 是为了让旁瓣和零点结构完整显示如果只画线性幅度第一旁瓣看起来只是一个小突起很容易误导旁瓣水平的判断。波束宽度可从图中直接读N16、dλ/2、θ020° 时主瓣 -3 dB 宽度约 6.4°~7°因为扫描到 20° 时波束会略微展宽和理论公式 0.886λ/(N·d·cosθ0) 对上才算代码可信。提示方向图输出是匹配滤波后的复值不能先用 abs 再平方来替代 abs(pattern)那会在较低电平上引入错误的刻度。2.3 用主瓣、旁瓣和零点验证单波束代码代码跑通后第一件事不是做多波束而是验证单波束方向图的三个特征量。主瓣峰值位置用max找到峰值角度应该和 theta0_deg 一致在扫描网格精度内[~, idx_max] max(abs(pattern)); fprintf(实际峰值角度: %.2f°\n, theta_deg(idx_max));第一旁瓣电平不加窗时均匀线阵理论第一旁瓣是 -13.26 dB。如果图里明显高于这个值多半是扫描角度步进太粗导致峰值采样不足明显低于这个值可能是无意中加了窗。零点位置方向图零点出现在 sinθ sinθ0 ± mλ/(Nd)m 为整数。我习惯同时画一张以 sinθ 为横轴的方向图周期性结构在 sinθ 坐标里更直观figure(Color,white); plot(sin(theta_rad), pattern_dB, LineWidth, 1.2); xlabel(sin\theta); ylabel(归一化增益 (dB)); ylim([-60 0]); grid on;在 sinθ 坐标系里均匀阵列方向图变成周期函数零点间隔均匀栅瓣位置一目了然。这一步验证习惯建议一直保留它不增加多少代码量却是最快区分“代码写错”和“参数选错”的手段。把单波束验证做扎实后面多波束、FFT 实现、窗函数设计才能在一个可信的基线上推进。3. 多波束形成两种套路权值矩阵法与 FFT 快速法3.1 权值矩阵法一次矩阵乘法得到全部波束输出单波束跑通后多波束就是往多个方向复制权值向量。假设希望在 [-40° -20° 0° 20° 40°] 五个方向同时形成波束常见做法是构造一个 N×M 的权值矩阵 W每个快照向量 xN×1直接做矩阵乘法 y W x得到 M×1 输出。下面是单个快照的完整代码% 多波束指向角集合 theta_list [-40 -20 0 20 40]; M length(theta_list); % 权值矩阵每一列对应一个波束的权值 W exp(1j * 2 * pi * d_lambda * n * (theta_list * pi / 180)); % N×M % 模拟一个从20°方向入射的远场信号快照 theta_src 20; s exp(1j * 2 * pi * d_lambda * n * sin(theta_src * pi / 180)); % N×1 % 各波束输出复数 y_beam W * s; % M×1 % 打印各波束输出功率相对值含阵列增益 for ii 1:M fprintf(指向 %5.1f° : %7.2f dB\n, theta_list(ii), ... 10*log10(abs(y_beam(ii))^2)); end逻辑说明W的每一列是某个波束的权值W做共轭转置后第 k 行输出本质是“第 k 个波束权值”与“入射信号快照”的内积。没有加噪声和阵元方向图时输出幅度只由信号来向相对各波束指向的方向图增益决定20° 入射时指向 20° 的波束输出最大其余波束依次衰减衰减量正好等于方向图上对应角度的相对增益。这里有一个容易误解的地方打印的 dB 值是未归一化的阵列输出功率包含 20·log10(N) 的阵列处理增益所以不要拿它和方向图里的 0 dB 主瓣直接对比要对比的是不同波束之间的相对值。多快拍场景把快照排成矩阵 X [x1, x2, …, xT]N×T输出 Y WX一个波束占一行时间序列沿列方向。这样比 for 循环逐快拍计算快得多而且代码更接近工程中“先收集数据块、后做波束形成”的处理方式。要检查某一路波束的时间波形直接取Y(k, :)即可。3.2 FFT 快速波束形成的映射关系与边界当阵元间距 d λ/2 时权值矩阵其实和 DFT 矩阵有相同的结构第 k 个波束权值的相位是 2π·n·k/N 的线性递增。因此 N 个波束可以用一次 N 点 FFT 全部算出来这就是 FFT 波束形成。映射关系为 sinθ_k 2k/Nθ_k asin(2k/N)k 取 -N/2 到 N/2-1对应 fftshift 后的索引% 前提d_lambda 0.5即 d λ/2 X_fft fftshift(fft(s, N)); % N点FFT后循环移位 k (-N/2 : N/2-1).; % 每个频点对应的波束指向角 theta_fft asin(2 * k / N) * 180 / pi; % 只保留物理有效方向|sinθ| 1 valid abs(2 * k / N) 1; theta_fft_valid theta_fft(valid); X_fft_valid X_fft(valid); fprintf(FFT波束形成有效指向角:\n); disp(theta_fft_valid);参数说明关键在fftshift。FFT 原始输出 k 从 0 到 N-1对应 sinθ 从 0 先增大到接近 2 再绕回不经过移位的话角度映射完全错位fftshift后 k 变成对称区间才能和 [-1,1] 的 sinθ 一一对应。有效波束数量略小于 N因为部分 k 对应的 sinθ 绝对值超过 1在实物空间没有对应方向。还要注意角度间隔并不均匀法线附近 sinθ 间隔均匀端射区由于 asin 的非线性角度间隔明显变大。这不是代码 bug而是 FFT 波束形成内部固有的映射特性。3.3 两种实现怎么选波束自由度与计算量权值矩阵法灵活度最高可以自由指定任意角度集合包括非均匀间隔、局部加密、指向角不落在 FFT 网格上代价是计算量随波束数 M 线性增长。FFT 方法快适合全空域覆盖且波束等 sinθ 间隔的场景但实际指向角被 FFT 网格限定端射区角分辨率变差。工程上常见做法是前级用 FFT 快速扫描粗搜检测到目标后切换权值矩阵法做细测角两者是互补关系而不是取代关系。从自由度角度理解N 元阵列最多形成 N 个独立正交波束超过 N 个波束并不增加信息量只是对同一空域做不同加权组合FFT 方法输出的 N 个波束恰好就是一组正交基权值矩阵法可以在这组正交基里自由旋转。所以在仿真代码里建议两种都实现互相验证对同一份快照数据用 FFT 方法和权值矩阵法在各自主瓣方向上的输出幅度应完全一致差异超过 0.1 dB 就要检查映射表或共轭方向。这个一致性检查我每次跑新项目都会做一次是防止方向图代码“看起来对但实际错”的后悔药。4. DBF 仿真参数怎么定阵元数、间距、窗函数与测角精度4.1 阵元间距与栅瓣dλ/2 方向图会冒出的假主瓣方向图在 sinθ 变量下是周期函数周期为 λ/d。当 d λ/2 时周期小于 2在 [-1,1] 物理区间内出现多个周期这就是栅瓣。栅瓣位置满足 sinθ_g sinθ0 ± mλ/dm 为整数。在 MATLAB 里验证栅瓣最快的方法是把 d_lambda 改成 0.8 或 1.0重新画方向图马上能看到第二个主瓣冒出来。所以无论仿真还是工程阵元间距硬约束是 d ≤ λ/2。注意这是通带不模糊条件宽带系统还要考虑最高频率分量对应的最短波长设计时要按最短波长定间距。栅瓣不只在仿真里影响美观更直接影响多波束测角如果主瓣和栅瓣同时覆盖到目标测角结果会出现模糊。仿真中排查栅瓣的标准动作是画 sinθ 横轴方向图看全区间是否有多个等高峰。还有一个容易忽略的点扫描角只在 [-90°, 90°] 内画图时部分栅瓣可能落在区间外没显示换个扫描范围立刻现形。建议第一版代码就把扫描范围扩大到 [-90°, 90°] 全区间不要自欺欺人。4.2 阵元数决定波束宽度和独立波束数上限阵元数 N 决定半功率波束宽度法线方向近似为 BW_3dB ≈ 0.886λ/(N·d)弧度。N16、dλ/2 时约 6.3°N32 时约 3.2°。波束越窄测角精度越高但覆盖同样空域需要更多波束。这里有一个常见的工程权衡表阵元数 N法线波束宽度(°)约可布独立波束数典型用途812.86~7简易测向、车载雷达166.414~16相控阵雷达、5G阵列323.228~30高精度测角、电子对抗641.655~60大型相控阵、射电阵列实际阵列能布多少个波束不是简单等于 N因为相邻波束之间要留交叠保证无盲区覆盖。多波束 DBF 一般按 1.5~2 倍波束宽度为间隔布波束所以全空域覆盖的波束数大致是总覆盖角度除以波束间隔。仿真里想要“每个波束都有效”注意波束越靠近端射方向展宽越明显间隔可以放宽法线附近加密这是均匀阵列在 sinθ 空间的固有属性。4.3 锥削窗压低旁瓣的代价是主瓣展宽和测角曲线变缓均匀权值旁瓣 -13.26 dB很多场景不够用。加窗锥削可以在权值上点乘窗函数% 在既有波束权值上叠加Hamming窗 w_win w .* hamming(N); pattern_win w_win * A; figure(Color,white); plot(theta_deg, 20*log10(abs(pattern_win)/max(abs(pattern_win))), LineWidth, 1.2); ylim([-60 0]); grid on; title(Hamming窗后的方向图);Hamming 窗第一旁瓣约 -43 dB主瓣宽度相对均匀窗展宽约 1.5 倍Blackman 窗旁瓣更低但主瓣展宽更多。如果你做的是多波束比幅测角旁瓣太低不一定是好事比幅测角依赖相邻波束方向图的交叠区斜率主瓣展宽会改变交叠点位置和斜率测角曲线需要重新标定。我在实际项目中吃过这个亏换了窗函数后测角误差从 0.5° 漂到 1.5°重新标定后才拉回来。所以窗函数设计要和测角算法一起调不能单独看着方向图好看就换。仿真里务必把加窗和未加窗两组方向图对比着看确认主瓣展宽可以接受。注意加窗后的阵列增益比均匀窗低因为部分阵元被加权削弱等效阵元数下降。这个损失在低信噪比场景下可能比旁瓣干扰更致命仿真时要同时评估输出 SNR 而不是只看旁瓣抑制。5. 多波束 DBF 仿真避坑5 个常见问题与排查思路5.1 方向图峰值不在预期指向sin/cos 约定不统一现象设定波束指向 20°方向图峰值出现在约 70° 或根本找不到明显主瓣。原因绝大多数是角度坐标系约定混用有的资料以阵列端射方向为 0° 用余弦函数 a(θ) exp(j2πd·n·cosθ)有的以法线为 0° 用正弦函数两套代码拼接时没统一。解决在工程代码开头注释写明“本工程统一以阵列法线为 0°正角度在法线右侧”所有函数只用 sin(θ) 计算相位差。如果来源公式不得不与约定不一致宁可重写也不在原有式子上做临时角度变换那种补丁最容易在后期引入隐蔽误差。仿真验证时用find(theta_deg theta0_deg)精确检查峰值位置是否落在 0.1° 以内偏差大直接查公式不要急着调参数。5.2 扫描范围内出现等高峰栅瓣作怪现象方向图在预期主瓣之外出现一个幅度几乎相等、间隔规律的峰值尤其是在大角度扫描时。原因d_lambda 大于 0.5主瓣和栅瓣同时进入可视范围。解决先把 d_lambda 改回 0.5 确认栅瓣消失如果系统设计要求更大的实际间距就要接受栅瓣存在并改用非均匀阵列稀疏阵、非线性阵打破周期性。仿真阶段的原则是先保证均匀线阵模型完全正确再引入非均匀阵结构不要一步到位写稀疏阵代码否则无法区分是稀疏设计问题还是基础代码问题。画 sinθ 横轴图后栅瓣会表现为等间隔重复峰这个特征可以快速和旁瓣区分。5.3 多波束输出幅度不一致无法直接比幅测角现象同一信号从 10° 入射指向 10° 的波束输出比指向 0° 的波束高但两者的差值不是方向图上直接读出的相对增益甚至出现相邻波束输出反转。原因仿真里没有计入阵元方向图比如实际阵元不是全向而是余弦方向图 cosθ或者波束数目少、交叠区增益跌得太快又或者你在比较时用了未归一化的输出功率。解决仿真中如果假设全向阵元各波束峰值增益理论上一致输出比值可以直接用方向图解释如果加入阵元方向图要在每个波束输出上乘以阵元增益因子并在比幅测角时用“比值-角度”查表曲线而不是理论公式。我一般会在仿真里同时保存方向图矩阵和波束输出向量逐个角度核对出现不一致先看是不是扫描步进太粗导致峰值采样偏差。5.4 FFT 波束形成在端射区角度错位现象用 FFT 方法做多波束信号从 60° 入射时输出最大点对应的索引换算出来的角度不是 60° 而是 50° 或 70°。原因一是漏了fftshift导致 k 与角度映射错位二是用线性映射 θ k·λ/(Nd) 替代了正确的 asin 映射。FFT 波束形成的指向角必须用 asin(2k/N) 计算尤其在端射区 asin 非线性明显线性近似误差可达数度。解决打印theta_fft_valid与信号真实来向对比确认有效波束指向角表和预期一致。注意 FFT 波束形成只能输出等 sinθ 间隔的波束端射区角度间隔大信号来向落在两个波束之间时能量会分散到相邻两个波束这是物理规律不是代码问题。要获得端射区更密的覆盖应改用权值矩阵法。5.5 快照矩阵维度方向或共轭方向弄反现象执行 Y W * X 时报维度错误或者不报错但输出全是乱码级小数值。原因X 的维度约定不统一。我统一约定“每列是一个快照每行是一个阵元”这样 X 是 N×TW 是 M×NY 是 M×T。如果你把 X 写成 T×N矩阵乘法要么不匹配要么结果语义完全相反。解决在构建快照矩阵后立即用size打印检查注释写明“X: N×T, W: N×M, Y: M×T”。还有一类问题是用w.普通转置替代w共轭转置导致相位补偿方向反了方向图变成陷零而不是峰值。遇到输出幅度异常偏小优先检查这两个地方。6. 把 DBF 仿真往工程推进自适应零陷和测角验证6.1 在常规方向图上叠加一个自适应零陷常规多波束只能固化为指定权值遇到强干扰时容易在干扰方向保留高旁瓣。工程里常见做法是在权值基础上叠加一个自适应修正项把零陷压到干扰方向。最小均方LMS类方法虽然简单但在仿真里用来验证“权值可自适应更新”的链路很有价值% 在20°目标附近放一个-50dB干扰用LMS微调权值 rng(2024); target exp(1j*2*pi*d_lambda*n*sin(20*pi/180)); interf exp(1j*2*pi*d_lambda*n*sin(-30*pi/180)); x_noisy target interf 0.01*(randn(N,1)1j*randn(N,1)); w_adapt w; % 以20°常规权值初始化 mu 0.01; % 步长过大发散过小收敛慢 for iter 1:200 y_out w_adapt * x_noisy; error conj(y_out) * x_noisy; % 用输出共轭乘输入近似梯度方向 w_adapt w_adapt - mu * error; end pattern_adapt w_adapt * A;逻辑说明这里用的是简化的 LMS 变体真正工程上更常用 MVDR 或对角加载后的采样协方差求逆但 LMS 代码短、便于观察零陷形成过程。步长 mu 是最敏感的参数太大会发散太小要几千次迭代。跑完画出pattern_adapt的方向图-30° 方向应出现明显零陷20° 主瓣基本保持。这个练习做完再换 MVDR 就很容易理解为什么协方差矩阵需要对角加载来保证数值稳定。6.2 用蒙特卡洛跑一次测角精度验证自适应零陷验证的是单一场景而 DBF 仿真最终要回答的问题是“这套参数测角能到什么精度”。我习惯最后跑一组蒙特卡洛固定信噪比随机生成 500 次带噪声快照用多波束比幅测角统计误差均值和标准差。以 ±30° 内布 5 个波束为例测角误差约在 1/20 波束宽度量级。这个数字会随窗函数、波束间隔和 SNR 明显变化建议每次改参数后都重新统计。做完这步你的 DBF 仿真就不只是画图好看而是能支撑方案论证和参数选型的定量依据了。我自己的教训是单次方向图验证通过的项目蒙特卡洛往往能暴露比幅测角在高 SNR 下的曲线拐点问题早跑早省心。希望帮到你。本文还有配套的精品资源点击获取