常规波束形成MATLAB实现：从原理、代码到MVDR升级指南

简介这是一份用于常规聚焦波束形成的MATLAB仿真程序面向水声信号处理、阵列信号处理等方向的学习者与研究人员可帮助理解波束形成的基本原理与实现流程。程序以M11元均匀线阵为例通过设置声速、采样频率、阵元间距及声源位置等参数依次完成阵列坐标生成、声程差计算、信号与噪声模拟、接收数据自协方差矩阵求解并在扫描范围内完成聚焦波束扫描成像最终输出归一化后的三维mesh声图和二维伪彩图同时提供峰值定位与主瓣宽度估计。压缩包共1个doc文件约29KB文档中包含完整MATLAB脚本及关键步骤注释结构清晰便于直接运行和二次修改适合作为课程设计、毕业设计或论文仿真的参考。资源已有1678人学习适合正在学习常规波束形成、需要结合仿真理解空间谱成像的入门及进阶用户。1. 常规波束形成没那么玄十几行MATLAB就能跑通常规波束形成CBFConventional Beamforming是阵列信号处理里最不值得恐惧的一块它不做任何统计自适应只靠“把阵元信号对齐到某个方向”这一件事就能换来整排天线阵的阵列增益。配合MATLAB十几行代码就能画出波束图、扫出信号来向整个过程几乎没有黑匣子。可现实中很多人一上来就跳去自适应波束形成被协方差矩阵求逆和特征分解搞得一头雾水回头才发现CBF的空间谱才是所有高级算法的白盒基准。这篇笔记面向两类人正在做DOA估计或阵列接收的工程师以及课程作业里被“波束形成”困住的学生。目标很简单——用一段可复现的MATLAB程序把常规波束形成的原理、参数、坑和下一步升级一次讲透。2. 常规波束形成原理与最小可运行代码从对齐时延到相位旋转2.1 从时延对齐到相位旋转CBF的数学骨架阵列信号处理里有个朴素到容易被忽略的事实把多根天线收到的信号对齐后再合并得到的信噪比提升是单天线无法获得的。常规波束形成要做的正是这件事。考虑一个等间距排列的N元均匀线阵阵元间距为d一个远场窄带信号从θ方向入射。以第一个阵元为参考第n个阵元相对它的传播时延是τ_n nd sinθ / c。如果信号是窄带的这个时延在频域里等价于一个相位旋转φ_n -2π nd sinθ / λ。CBF做的就是在假设来向为θ_s时把每个阵元的相位反向补偿回去让N路信号同相叠加。把对某个方向θ_s的补偿写成向量a(θ_s) [1, exp(-j2πd sinθ_s / λ), …, exp(-j2π(N-1)d sinθ_s / λ)]^T这就是阵列导向矢量。波束输出为y a^H(θ_s) x当θ_s和真实来向θ一致时各路信号对齐输出模值最大理想情况下输出幅度是单阵元的N倍功率是N^2倍。当θ_s偏离真实来向时相位补偿不匹配输出功率下降。把θ_s扫过整个空域得到的空间谱就是CBF的波束图。这套思路在教科书里叫延迟-求和波束形成在工程里直接叫常规波束形成或时延相加波束形成。CBF之所以值得先跑通是因为它不依赖信号统计特性只依赖阵列几何和来向假设结果可以手算验证。我接触过的自适应算法项目几乎都是从一版CBF空间谱开始的——先确认数据链路和阵列流型没问题再上复杂度高的算法否则翻车时根本分不清是数据问题还是算法问题。两个默认前提要记牢。第一是窄带假设CBF的相位补偿只在单一频率上严格成立宽带信号必须分频带处理或用FIR时延滤波器代替移相器否则波束会发散。第二是阵元间距均匀线阵默认d λ/2这是空间采样定理的结论。超过半波长会出现栅瓣小于半波长则孔径利用率下降、分辨率变差。这两个前提在后面的参数章节里会反复遇到。2.2 12行MATLAB跑通CBF最小可运行版本下面是最小可运行版本不追求效率只求把CBF的骨架讲清楚。信号用复基带随机序列模拟主程序加一个扫描循环总共不到20行。% 最小可运行的常规波束形成CBF示例 % 均匀线阵 远场窄带信号8阵元半波长为间距 c 3e8; % 光速m/s fc 2.4e9; % 工作频率Hz lambda c / fc; % 波长m d lambda / 2; % 阵元间距半波长 N 8; % 阵元数 theta_true 30; % 信号真实来向度 n_snap 1024; % 快拍数 rng(0); % 复基带随机信号N路阵元共享同一个包络 s (randn(1, n_snap) 1j*randn(1, n_snap)) / sqrt(2); % 真实来向的导向矢量 a_true exp(-1j * 2 * pi * d * (0:N-1). * sind(theta_true) / lambda); % 模拟接收数据N x n_snap x a_true * s; % 扫描角度 theta_scan -90:0.5:90; p zeros(size(theta_scan)); for k 1:numel(theta_scan) a_scan exp(-1j * 2 * pi * d * (0:N-1). * sind(theta_scan(k)) / lambda); y a_scan * x; % 波束加权求和 p(k) mean(abs(y).^2); % 快拍平均功率 end p_dB 10 * log10(p / max(p)); plot(theta_scan, p_dB, LineWidth, 1.5); xlabel(扫描角度度); ylabel(归一化功率dB); grid on; ylim([-40 5]);这段代码里(0:N-1).生成N×1的阵元序号列向量是MATLAB写导向矢量最常用的技巧后续写MVDR、ESPRIT都会复用。sind()保证角度单位是度如果换成sin()且角度按度传入峰值位置会偏到完全对不上。a_scan * x是共轭转置加权求和本质是空间匹配滤波mean(abs(y).^2)对快拍维取平均功率最后归一化成dB方便看旁瓣相对电平。实际调试时值得动三个参数N从8改成4或16直接观察主瓣变宽还是变窄theta_scan的步进从0.5°改成1°或0.1°观察峰值是否抖动d从λ/2改成0.7λ观察第二个峰何时冒出来。每改一次重新看图比背公式有用得多。提示a_true和a_scan的相位符号只要收发一致峰值位置就正确如果两处符号同时取反空间谱会在-30°出现峰值这是正常现象不是bug。想观察噪声影响把x替换成x_noisy x 10^(-snr_db/20) * (randn(N, n_snap) 1j*randn(N, n_snap)) / sqrt(2)SNR越低波束图底噪越高但峰值位置基本不动这就是CBF稳健性最直观的体现。3. 把CBF封装成MATLAB函数签名设计、三个调用场景与输出约定3.1 封装成cbf_scan函数签名、输入约定与参数选择脚本能跑通只是第一步。实际工程里DOA扫描、空域滤波、多波束接收三个场景的代码主体完全一样差别只在加权向量和输出维度。我一般会把CBF封装成独立函数避免每个脚本复制粘贴一长串扫描循环。另一个更现实的理由是后续切换到自适应波束形成时函数里只改一行加权向量CBF和MVDR的实测对比成本会低很多。function [p_dB, theta_out] cbf_scan(x, d, lambda, theta_scan) % 常规波束形成空间谱扫描 % 输入: % x N x snap 复数阵列数据N为阵元数 % d 阵元间距米 % lambda 波长米 % theta_scan 1 x M 扫描角度度 % 输出: % p_dB 1 x M 归一化空间谱dB % theta_out 1 x M 扫描角度度 N size(x, 1); theta_out theta_scan; p zeros(size(theta_out)); for k 1:numel(theta_out) a exp(-1j * 2 * pi * d * (0:N-1). * sind(theta_out(k)) / lambda); a a / N; % 均匀加权归一化 p(k) mean(abs(a * x).^2); end p_dB 10 * log10(p / max(p)); end这里的关键设计是a a / N把加权向量归一化。归一化后输出量纲是阵元平均功率峰值位置不变方便后续和MVDR等算法横向对比如果不归一化峰值功率是N^2倍dB图只差一个常数但在对比输出信干噪比时会非常别扭。参数选择上theta_scan的步进一般先用0.5°或1°做粗扫锁定峰值区间后再用0.05°细扫。x必须是复数IQ结构如果数据文件里只有实部不要直接喂给这个函数先做希尔伯特变换恢复解析信号。3.2 调用场景A全空域扫描与DOA估计最常见的调用方式是用CBF做全空域扫描直接读空间谱峰值位置作为DOA估计结果。假设有两个等功率信号分别从-20°和10°入射封装之后的调用代码很干净n_snap 1024; theta1 -20; theta2 10; d 0.0625; % 2.4GHz下半波长米也可用lambda/2计算 lambda 0.125; a1 exp(-1j * 2 * pi * d * (0:N-1). * sind(theta1) / lambda); a2 exp(-1j * 2 * pi * d * (0:N-1). * sind(theta2) / lambda); s2 (randn(2, n_snap) 1j*randn(2, n_snap)) / sqrt(2); x2 a1 * s2(1,:) a2 * s2(2,:); [p_dB, theta] cbf_scan(x2, d, lambda, -90:0.5:90); plot(theta, p_dB);这个场景的边界条件很直观两个信号角度间隔40°CBF能清晰分辨出两个峰。如果间隔缩小到5°8元阵的波束宽度大约12.7°两个峰会合并成一个宽包络峰值位置落在中间这就是CBF的分辨率极限。遇到这种情况要么增加阵元数要么换高分辨率算法靠CBF硬扛是扛不过去的。3.3 调用场景B固定方向滤波与多波束接收DOA扫描输出的是功率值但很多场景需要直接输出某个方向的信号时间序列。固定方向空域滤波的封装更简单function y_beam cbf_filter(x, d, lambda, theta) % 固定方向波束输出返回对齐后的时间序列 N size(x, 1); a exp(-1j * 2 * pi * d * (0:N-1). * sind(theta) / lambda); w a / N; y_beam w * x; % 1 x snap end使用时如果只想提取10°方向的信号调用y_beam cbf_filter(x2, d, lambda, 10)即可。得到的y_beam是加权求和后的复时间序列可以直接做解调、频谱分析或后续的编码解码而不是只有一个功率值。这里要注意CBF对非期望方向没有主动抑制能力干扰信号只是被波束方向图的固定旁瓣衰减典型旁瓣电平约-13dB如果干扰很强泄漏仍然可观这个场景正是第4章加窗和第6章MVDR的意义所在。多波束接收是雷达和声呐里最常见的形态多个固定波束并行输出thetas [-20, 10, 40]; Y zeros(length(thetas), size(x2, 2)); for i 1:length(thetas) a_i exp(-1j * 2 * pi * d * (0:N-1). * sind(thetas(i)) / lambda); Y(i, :) (a_i / N) * x2; endY的每一行就是一个方向的波束输出所有波束并行工作。需要提醒的是这里每个波束都没有做旁瓣抑制如果想压低旁瓣就得把a_i / N替换成导向矢量乘以窗向量再归一化这正是下一章的内容。4. 波束图三个必调参数阵元数、窗函数、阵元间距4.1 阵元数N决定角度分辨力一张表摸清孔径上限均匀线阵、均匀加权的半功率波束宽度近似为BW_3dB ≈ 0.886 * λ / (N * d)单位弧度。当d λ/2时λ/(N*d) 2/N于是得到一张非常实用的对照表| 阵元数 N | 半功率波束宽度 | 大致可分辨的最小角度间隔 | | 4 | 约 25.4° | 约 ±12.7° | | 8 | 约 12.7° | 约 ±6.3° | | 16 | 约 6.4° | 约 ±3.2° |这张表的含义很朴素物理孔径N*d决定了角度分辨力的天花板加窗还会让主瓣变宽约1.2到1.6倍实际分辨力要再打折。比如想分辨间隔5°的两个目标8元阵勉强够16元阵才比较舒服。工程上常见的错误是只加阵元数却不考虑波长同样的N和d工作频率升高后λ变小波束宽度会显著变窄这既是好事也是坑——高频段阵列的角度校准误差会更敏感。4.2 窗函数控制旁瓣加窗代码与主瓣代价均匀加权的第一旁瓣电平是-13.2dB在很多场合高得不能接受。加窗的本质是对各阵元幅度再加一组权压低阵元域上的副瓣泄漏。工程上最常用的是汉宁窗和切比雪夫窗对比关系如下| 窗类型 | 第一旁瓣电平 | 主瓣宽度倍数相对均匀加权 | | 均匀 | -13.2 dB | 1.00 | | 汉宁周期 | -31.5 dB | 约 1.62 | | 切比雪夫 30 dB | -30 dB | 约 1.18 | | 切比雪夫 40 dB | -40 dB | 约 1.28 |切比雪夫窗的最大优势是旁瓣电平可精确设置N固定时想让旁瓣是-30dB还是-50dB都能做到。代价是主瓣宽度随旁瓣要求提高而变宽而且N较小时边缘阵元权重会很小等效孔径缩水。加窗的MATLAB代码非常短% 构造一个含30°信号的接收数据x s30 (randn(1, n_snap) 1j*randn(1, n_snap)) / sqrt(2); a30 exp(-1j * 2 * pi * d * (0:N-1). * sind(30) / lambda); x a30 * s30; % 均匀窗加权 w_uniform a30 / N; % 切比雪夫窗40dB旁瓣 win_cheb chebwin(N, 40).; w_cheb (a30 .* win_cheb(:)) / sum(win_cheb); % 对比两种加权的输出功率 y_uniform abs(w_uniform * x); y_cheb abs(w_cheb * x);窗向量是实数导向矢量是复数用点乘把窗应用到每个阵元上再做求和归一化。从频域角度看加窗相当于对空间谱做周期卷积效果是旁瓣压低、主瓣展宽。N小于8时切比雪夫窗的边缘阵元权重会掉得非常厉害波束宽度比理论值宽不少这是模板里最容易被忽略的边界条件。4.3 阵元间距与栅瓣半波长为什么是红线阵元间距d大于λ/2时除了真实来向外某些扫描角会得到同样高的输出这些位置叫栅瓣。栅瓣位置由sin(θ_g) sin(θ_true) mλ/d决定m为整数。用d 1.5λ、θ 30°代入λ/d 0.667m -1时sin(θ_g) 0.5 - 0.667 -0.167对应约-9.6°。也就是说波束图上除了30°的峰还会在-9.6°附近冒出几乎等高的假峰。栅瓣比旁瓣危险得多旁瓣可以用窗压低栅瓣却和主瓣差不多高加窗根本压不掉只能改小阵元间距或限制扫描范围。代码里最常见的错误是图省事直接把d设成一个整数波长而不是λ/2结果扫描图上目标变成两个。把d改回lambda/2问题立刻消失。实际操作中我一般先按半波长起算再根据安装空间做权衡如果间距确实要拉大就先把可能出现的栅瓣位置算一遍确认它落在扫描范围之外再定方案。注意阵列间距是空间采样频率的物化表现半波长是空间奈奎斯特条件不是经验值。任何超过半波长的间距改动都要先做栅瓣位置计算。5. 常规波束形成MATLAB实现避坑指南5个高频翻车点CBF本身不复杂但实测翻车率极高。下面5条是MATLAB实现里最常踩的坑按现象、原因、解决三段式整理基本覆盖从纯仿真到半实物过渡的全部高频问题。5.1 翻车点1角度单位混用峰值偏移1°起步现象扫描角写成-90:0.5:90信号真实来向30°但波束图峰值出现在29°或31°偏差忽大忽小角度越大偏得越离谱。原因MATLAB的sin()和cos()默认输入弧度。代码里写成sin(theta)但theta是度数或者从其他语言迁移代码时把sind()误改成了sin()都会导致角度映射非线性失真。解决统一用sind()和cosd()。验证方法很简单——把峰值对应的扫描角打印出来和真实来向比对。半实物联调时我第一个检查项永远是角度单位。5.2 翻车点2复数IQ信号被取成实数方向性消失现象IQ数据进数组时写了x real(data)波束图变得几乎全向输出功率对所有扫描角都差不多完全没有方向性。原因阵列增益靠复数相位对齐实现取实部等于丢掉了一半的相位信息加权求和从幅度叠加退化成部分抵消。这不是算法问题是数据已经不再满足窄带复数模型。解决接收链路保持复数结构x complex(I, Q)。如果数据文件里只有实部先用希尔伯特变换或正交混频恢复解析信号再进CBF。这个坑在自采数据里出现率最高因为很多人习惯先把ADC数据存成real数组。5.3 翻车点3快拍数太少波束图旁瓣乱跳现象空间谱整体形状对但旁瓣起伏严重换一批数据运行结果明显不一样时好时坏。原因CBF输出功率是对快拍维做样本平均样本数太少时统计方差大。快拍数低于几十时旁瓣起伏会完全掩盖真实旁瓣电平。解决快拍数加到512或1024这是最直接的办法。如果数据本身很短就只读峰值位置不要读旁瓣电平。注意CBF不涉及协方差矩阵求逆不需要做对角加载或矩阵平滑直接对输出功率做时间平均即可。5.4 翻车点4近场信号当远场处理相位补偿全错现象把发射源放在离阵列很近的位置比如几米仿真或实测的波束指向总是偏峰值变钝位置还随快拍移动。原因CBF导向矢量建立在平面波假设上。近场信号是球面波各阵元的时延差不再等于nd sinθ/c远场公式失效。解决先用远场条件判断源到阵列距离R必须满足R 2D²/λ其中D (N-1)d是阵列孔径。不满足时改用近场聚焦逐元计算真实距离差做相位补偿。工程速算的经验是先把D算出来如果2D²/λ和实际距离在一个数量级甚至更大直接按近场处理别想着用远场模型硬凑。5.5 翻车点5扫描网格太粗峰值存在量化误差现象真实来向21°扫描步进1°空间谱峰值在21°或22°定位误差0.5°真实来向正好落在两个网格点中间时峰值能量分散在相邻格点看起来像两个小峰。原因扫描网格本身就是一组离散假设。CBF不做插值网格多粗输出就有多粗。解决先粗扫找峰再在峰值附近±2°范围内用0.05°或0.1°细扫。如果对精度还有要求在细扫结果上做抛物线插值运算量几乎可以忽略。这个习惯比直接换高分辨率算法成本低得多实测里能省掉大量无意义的精度排查。6. 从CBF到自适应波束形成MVDR的最小改动路径6.1 只换一行加权向量MVDR的MATLAB实现CBF的白盒特性是它最大的价值先跑出一版空间谱确认阵列没接错、角度没偏、数据格式没问题再往自适应算法上走。MVDR最小方差无失真响应是升级成本最低的一步它和CBF的差别只体现在加权向量上。% 沿用前文的8元均匀线阵信号从30°入射 a_s exp(-1j * 2 * pi * d * (0:N-1). * sind(30) / lambda); % 接收数据x包含30°期望信号和-10°强干扰 Rxx x * x / size(x, 2); % 对角加载避免协方差矩阵奇异 Rxx_dl Rxx 1e-3 * trace(Rxx) / N * eye(N); % 常规波束形成加权 w_cbf a_s / N; % MVDR加权期望方向增益为1同时最小化输出功率 w_mvdr (Rxx_dl \ a_s) / (a_s * (Rxx_dl \ a_s)); % 全空域扫描对比两种算法的空间谱 p_mvdr zeros(size(theta_scan)); for k 1:numel(theta_scan) a_k exp(-1j * 2 * pi * d * (0:N-1). * sind(theta_scan(k)) / lambda); w_k (Rxx_dl \ a_k) / (a_k * (Rxx_dl \ a_k)); p_mvdr(k) mean(abs(w_k * x).^2); endMVDR的思路是在期望方向增益等于1的约束下最小化输出功率w^H Rxx w解就是这个式子。它和CBF的本质差距体现在干扰方向上CBF的波束方向图对干扰无能为力固定旁瓣能衰减多少算多少MVDR会自动在干扰来向形成深零陷把干扰功率压到噪声底。代价是它对阵列流型误差和协方差估计误差敏感所以对角加载系数是必加的实测中1e-3不够稳时可以直接加到1e-1。6.2 验证方法从空间谱对比到信干噪比对比验证MVDR正确性的方法非常直接。构造一个30°期望信号加一个-10°强干扰的接收数据分别用CBF和MVDR扫描全空域。CBF的空间谱里-10°位置通常会出现一个峰如果干扰功率比期望信号高这个峰甚至超过期望峰MVDR的空间谱里-10°位置会出现一个深达-40dB以下的凹坑而不是峰。这一步能跑通就说明自适应置零机制真正生效了。我现在的习惯是任何阵列实验第一步永远先跑一版CBF画波束图、做全空域扫描确认数据链路是通的第二步才上MVDR或其他高分辨率算法而且一定把CBF的空间谱留在同一个坐标系里做对照。这个过程的代码改动通常不超过20行却能把波束形成从黑匣子变成可视化工具排查问题时的收益远超成本。CBF是起点不是终点但跳过他直接去追自适应算法往往只会收获一屏幕无从下手的复数矩阵。希望帮到你。本文还有配套的精品资源点击获取