简介程序实现了常规波束形成算法的完整MATLAB仿真流程面向水声信号处理、阵列信号处理方向的科研人员与初学者。代码依次完成声源信号模拟、阵元接收、噪声叠加、自协方差矩阵计算、聚焦扫描及声图绘制并以三维mesh与二维伪彩图展示结果便于直观理解波束形成原理和参数影响。程序内部预设声速1500、采样频率10000、信号频率500等关键参数同时阵元间距、阵元个数、扫描步长均可调整便于开展定位精度、旁瓣水平等对比试验。资源包仅1个文件为doc格式大小29KB内含可运行代码与关键步骤说明适合作为课程设计、毕业设计或算法入门的练习模板也可用于算法验证与二次开发。目前已有1678人学习下载是同类MATLAB波束形成资源中较受欢迎的参考范例。1. 常规波束形成这份 MATLAB 程序到底帮你干了什么拿到这段程序的第一印象是它短整个流程只有几十行却把常规波束形成从参数设置、阵列几何、信号仿真、协方差估计到空间扫描、声图绘制、主瓣参数提取全串起来了。c1500、fs10000、M11、d3再加上一个 (10,20,10) 的单声源程序会生成一张 x-z 平面的聚焦波束声图并自动算出峰值位置和半功率宽度。对刚接触波束形成的人来说它能帮你把教科书里的 R 矩阵和方向矢量变成一张看得见的图对已经在做声源定位的工程师来说它是一份可以改参数、加窗、扩展成三维扫描的基础骨架。真正上手之后你会发现这段代码最大的价值不是输出那个漂亮的亮斑而是把阵元间距、参考点、扫描步长这些容易被忽略的参数全部暴露在了你面前。2. 从信号模型到协方差矩阵为什么先要把数据算对2.1 参数设置声速、采样频率、阵元间距为什么这么选程序开头这一段是全部仿真的地基close all clear all clc c1500; % 水声环境中的声速1500 m/s fs10000; % 采样率 10 kHz远高于信号频率 500 Hz T 0.1; % 信号时长 0.1 秒 t 0:1/fs:T; % 时间轴0.1 秒内 1001 个点最后一个点会落在 t0.1 Llength(t); % 快拍长度 L1001 f500; % 窄带信号频率 w2*pi*f; % 角频率 kw/c; % 波数 k单位为弧度每米 M11; % 阵元个数 Nmid1; % 参考点选在第 1 个阵元 d3; % 阵元间距 3 m m[0:1:M-1]; % 阵元编号 0 到 10c1500 是水中声速的常用近似值。fs 取 10000 对 500Hz 信号来说足够奈奎斯特频率是 5000Hz不会出现采样混叠。真正值得多想的是 d3。频率 500Hz 对应波长 λc/f3m这里阵元间距恰好一个波长已经超出了常规线阵要求的小于等于半波长条件。很多初学程序为了凑整数喜欢把间距设成和波长一样结果后面声图会出现栅瓣。这不是程序 bug而是参数本身埋了一个可以观察的现象。我建议你先按 d3 跑一遍再把 d 改成 1.5 跑一遍对比两次声图空间采样定理的印象会刻进脑子里。f500 和 T0.1 决定了快拍数 L1001。后面 Rp*p/L 要用 L 做平均所以 L 不能太小。这里 1001 个点对协方差估计来说足够了。2.2 阵元坐标与声源位置声程差是波束形成的物理核心接着看几何部分yizeros(M,1); % 所有阵元的 y 坐标都是 0保留 M*1 列向量格式 zizeros(M,1); % z 坐标也是 0 xim*d; % x 坐标从 0 到 30间隔 3 xixi.; % 转成列向量得到 M*1 坐标矩阵 y120; x110; z110; % 声源位置 Ric1sqrt((x1-xi).^2(y1-yi).^2(z1-zi).^2); % 声源到每个阵元的斜距M*1 Rn1Ric1-Ric1(Nmid); % 各阵元相对第一个阵元的声程差这里构建的阵列是沿 x 轴方向的 11 元均匀线阵阵元位置分别为 (0,0,0)、(3,0,0)、(6,0,0) 一直到 (30,0,0)。声源放在 (10,20,10)距离参考阵元的斜距约为 24.5m距离最后一个阵元约为 30m。计算 Ric1 和 Rn1 的目的是得到每个阵元相对参考点的相位差。程序里 Nmid1 表示参考点是第一个阵元这是整段代码的相位基准。需要注意参考点取在阵元端点和取在阵列中心理论输出功率的绝对值不变因为共轭相乘后公共相位因子会抵消但中间计算的数值范围不同。后面避坑章我会专门讲参考点选择。2.3 信号、噪声和接收矩阵这个 SNR 不是严格意义上的分贝值信号构造部分是最容易让人误解的s1cos(w*t); % 参考阵元接收到的信号1*L snr 20; Am 10^(-snr/20); % 噪声幅值系数0.1 n1Am*(randn(M,L)j*randn(M,L)); % 复高斯噪声M*L p1zeros(M,L); for k11:M p1(k1,:)Ric1(Nmid)/Ric1(k1)*s1.*exp(-j*w*Rn1(k1)/c); end pp1n1; Rp*p/L;这段代码模拟的是每个阵元收到的信号幅度按距离比衰减相位按声程差延迟。Ric1(Nmid)/Ric1(k1)就是参考距离与实际距离之比物理上对应球面波扩散衰减。exp(-j*w*Rn1(k1)/c)是窄带相移它等价于时域中的延迟但在频域或窄带假设下直接用复指数表达更简洁。s1 在参考阵元上就是 cos(wt)所以其它阵元是它的幅度缩放和相位旋转版本。snr20 这一行要特别留意。Am10^(-snr/20)得到 0.1意思是噪声的幅度比例是 0.1。但随机数 randn 生成的是标准正态分布样本复噪声的实部虚部功率加起来并不是固定的 -20dB。严格定义信噪比应该用信号功率和噪声功率的比值noise_power signal_power*10^(-snr/10)。程序里这种写法在演示场景中不会出问题但如果你要做定量比较比如改变 snr 看主瓣变化就要改用功率反推噪声方差的方法否则结果会有随机起伏。Rpp/L 是样本协方差矩阵维度 MM。L 在这里是快拍数等于 1001。常规波束形成的输出是 bRb所以 R 的质量直接决定声图质量。如果 L 太小R 会出现秩亏缺旁瓣会乱跳。后面避坑章会讲一个 L1 的极端反例。2.4 为什么用协方差矩阵而不是直接对快拍求和波束形成的经典表达式是 P(θ)w^H R w其中 w 是聚焦方向矢量。有人会问为什么不用mean(abs(b*p).^2)其实两者数学等价。用 R 的好处是一旦你算出 R后面换 MVDR、MUSIC 这类基于特征分解的算法时不用重新组织数据。程序里对每一扫描点计算 bRb本质上是对 L 个快拍的输出功率取平均。我复盘这段代码时发现一个值得注意的点扫描方向矢量 b 里只包含了相位项没有距离衰减补偿项。生成信号时 p1 里面乘了Ric1(Nmid)/Ric1(k1)但扫描时 b 用的是exp(-j*w*Rn/c)没有乘类似的比例因子。这意味着波束形成模型默认是平面波正态但生成数据时是球面波。在近场仿真里这个不一致会让靠近阵列的网格产生偏高响应。如果你想把模型弄得自洽可以在 b 上乘一个Ri(Nmid)./Ri也就是对距离衰减做反向补偿。做完之后峰值位置不变但旁瓣分布会变这是加深理解的推荐实验。3. 扫描成像聚焦方向矢量与 201×201 网格的逐点计算3.1 扫描网格步长 0.1 米的代价和收益step_x0.1; step_z0.1; yy1; % 扫描平面 y 坐标固定为声源所在平面 x[0:step_x:20]; % x 方向 201 个点 z[0:step_z:20]; % z 方向 201 个点这段代码把扫描区域限定在 y20 的垂直平面内。x 和 z 各从 0 到 20步长 0.1所以每个方向上有 201 个点总网格数 201*20140401。对每个网格点程序要计算一次 M11 维的方向矢量并做一次二次型运算总循环量超过四万次。对现代电脑来说这是几秒钟的事情。但如果你把步长改成 0.01网格点就变成 2001*2001 约四百万个内存中光 Pcbf 矩阵就要占约 64MBdouble 类型计算时间会涨到十几分钟。所以步长是个需要权衡的参数粗扫用 0.1 或 0.2锁定峰值后局部细化到 0.01。不要一上来就追求高分辨率那会让自己在调试参数时浪费大量时间。3.2 内层循环到底做了什么核心扫描循环只有几行for k11:length(z) for k21:length(x) Risqrt((x(k2)-xi).^2(y-yi).^2(z(k1)-zi).^2); RnRi-Ri(Nmid); bexp(-j*w*Rn/c); Pcbf(k1,k2)abs(b*R*b); end end对每个网格点 (x,z)先计算该点到 11 个阵元的斜距向量 Ri再减去参考阵元的斜距得到相对程差 Rn。然后构造聚焦方向矢量 b最后用 bRb 求出该点的波束输出功率。这里 b 的维度是 M1R 是 MM所以 bRb 是一个标量取绝对值是为了去掉数值计算中可能残留的微小虚部。这段代码内层循环在 MATLAB 里不算快。如果你追求效率可以改成矩阵化版本预先把所有网格点的坐标向量化一次性计算距离矩阵再批量乘 R。但原版的逐点循环更易读也更容易在断点处检查单个网格点的数值。我一般建议第一次运行时保留循环理解后再做矩阵化优化。3.3 一个网格点的输出功率可以手动验证为了确认 Pcbf 计算没有出错你可以单独检查声源位置 (10,20,10) 的输出。理论上当扫描点正好等于声源位置时方向矢量 b 的相位应该与该阵元接收信号中的相位完全对齐因此 bRb 会达到最大值。如果扫描网格的分辨率足够峰值就会出现在最接近声源的那个网格点。在程序中声源位置 (10,20,10) 恰好落在网格上因为 x10 有对应的索引 101z10 有索引 101。所以峰值位置应该精确对应第 101 行、第 101 列。但由于坐标换算公式用了b/length(x)*20输出显示的坐标会偏移约 0.05m这点在避坑章会展开。3.4 为什么要归一化扫描得到的 Pcbf 是功率值幅度大小取决于信号强度、阵元数、快拍数。为了在不同参数下对比归一化是必要的。程序的处理是for k11:length(z) pp(k1)max(Pcbf(k1,:)); end PcbfPcbf/max(pp);先求每一行最大值再在行最大值中找全局最大。这个写法有点绕但效果等价于Pcbf Pcbf/max(Pcbf(:))。归一化后主瓣最大值是 1其他值在 0 到 1 之间这样就能用 0.5 作为半功率阈值来测量主瓣宽度。如果要看旁瓣归一化到分贝域更直观。4. 绘图与主瓣参数提取mesh、pcolor 和半功率宽度怎么算才靠谱4.1 两种声图三维曲面和二维伪彩色图绘图部分非常标准figure; mesh(x,z,Pcbf); xlabel(x(m)),ylabel(z(m)) title(常规聚焦波束形成单亮点声图) figure; pcolor(x,z,Pcbf); shading interp; xlabel(x(m)),ylabel(z(m)) title(常规聚焦波束形成单亮点声图); colorbarmesh 的优点是显示幅度的高度信息你能直观看到主瓣有多尖、旁瓣有多高。pcolor 加 shading interp 产生平滑的伪彩色图适合在报告里标注峰值位置。需要注意pcolor 的默认颜色映射是线性的。当主瓣很尖、旁瓣相对低时线性色标会让所有小幅值区域都变成同一个深色旁瓣细节全被吞掉。这种情况下我习惯用分贝显示Pcbf_dB 20*log10(Pcbf eps); figure; pcolor(x,z,Pcbf_dB); shading interp; colorbar; caxis([-40 0]);这样能清楚看到主瓣以下 40dB 范围内的旁瓣结构。你会在 d3 时看到除主瓣外的栅瓣分贝图比线性图明显得多。4.2 峰值坐标和半功率宽度这段代码有隐藏偏差程序结尾这段是很多人会直接复制的工具代码但里面藏着一个索引换算问题[a,b] find(Pcbfmax(max(Pcbf))); xm1b/length(x)*20, zm1a/length(z)*20, [a b]find(Pcbf 0.5*max(max(Pcbf))); deta_x(max(b)-min(b))/length(x)*20 deta_z(max(a)-min(a))/length(z)*20第一步找最大值位置。如果 Pcbf 中最大值只出现在一个网格点a 和 b 就是两个标量xm1 和 zm1 也是标量。但问题在于换算b 是列索引范围从 1 到 201。b1 对应 x0b201 对应 x20。正确的距离是(b-1)*step_x (b-1)*0.1。程序写成b/length(x)*20等于(b/201)*20当 b101 时得到 10.05m而正确值是 10m。这就是前面提到过的 0.05m 偏差。换成正确公式之后峰值位置就是精确的 10m。半功率宽度 deta_x 用的是最大最小索引差再乘以 20/201。因为步长 0.1 恰好等于 20/201 的近似所以对半功率宽度影响不大但为了代码自洽建议也改成(max(b)-min(b))*step_x。4.3 多个最大值会让 find 出错Pcbf 是离散网格上的数值理论上最大值点唯一。但在某些对称场景或者 Pcbf 出现平顶时可能有多个网格点同时取到最大值。这时find(Pcbfmax(max(Pcbf)))返回的 a 和 b 是向量xm1 就会输出一串数而不是一个坐标。这是新手最容易蒙圈的地方。更稳健的写法是[max_val, idx] max(Pcbf(:)); [peak_z_idx, peak_x_idx] ind2sub(size(Pcbf), idx); xm1 (peak_x_idx - 1) * 0.1; zm1 (peak_z_idx - 1) * 0.1;max 只返回第一个最大值的位置ind2sub 把线性索引换成行列索引。这样无论网格上是否有多个相等峰值输出都不会报错。如果你确实想统计所有峰值点可以用find(Pcbf max_val - 1e-6)。4.4 半功率宽度和理论值怎么对照理论上一维均匀线阵在扫描方向上的主瓣半功率宽度大约正比于 λ/D其中 D(M-1)*d 是阵列孔径。在程序默认参数下D30mλ3mλ/D0.1也就是约 0.1 弧度的角度分辨率。但程序输出的是投影到 x-z 平面的空间宽度不是一个角度不能直接把 λ/D 转换为米。我建议用仿真结果做相对验证固定其他参数把 M 从 11 改成 21观察 deta_x 是否接近减半。这个实验比背公式更有说服力。5. 避坑指南常规波束形成仿真里六个高频翻车点5.1 阵元间距等于波长声图出现栅瓣现象跑完程序pcolor 图上除了 (10,20,10) 附近的主瓣x 方向对称位置还出现了几乎等高的亮斑。第一次看到时我还以为是双声源定位成功。原因d3m 等于一个波长违反了均匀线阵空间采样条件 d≤λ/2。当阵元间距过大时不同方向的导向矢量会重复栅瓣在空间谱中形成虚假峰。解决把 d 改成 1.5m即半个波长再跑一次虚假亮点消失。如果不修改 d至少在解释结果时要把栅瓣标注出来别把它当成真实目标。5.2 峰值坐标总是偏差约 0.05 米现象声源在 (10,20,10)程序输出 xm110.05、zm110.05。原因坐标换算用了b/length(x)*20索引 1 被映射到 0.1m 而不是 0m所以所有坐标都偏移了约半个步长。解决改成(b-1)*step_x或(a-1)*step_z。这个坑不影响声图视觉但当你拿峰值坐标和真实声源位置比对时会被追问一句“你这里为什么差 0.05”非常尴尬。5.3 峰值位置变量突然变成一长串数组现象代码运行后 xm1 后面跟着一串数脚本也不报错。原因find 匹配到多个等于最大值的网格点。网格分辨率、插值方式或归一化精度都可能造成一对多的匹配。解决用 max(Pcbf(:)) 配合 ind2sub 提取第一个最大位置如果需要统计所有峰值用阈值代替严格相等。5.4 半功率宽度随噪声随机变化无法稳定比较现象同样的参数只改了噪声deta_x 可能从 3.2m 跳到 4m。原因randn 每次生成的噪声不同导致主瓣形状在阈值附近抖动离散网格落在 0.5 边界上的位置也随之变化。解决在代码开头加rng(0)固定随机种子。做参数对比实验时固定同一噪声矩阵保证差异只来源于你改的参数。这是研究工作中很重要的习惯能避免很多让人挠头的偶然现象。5.5 快拍数 L 等于 1 时声图爆炸现象把 T 改成 0.0001或者只取一个时间点声图会出现大量随机尖峰峰值位置也乱跳。原因L1 时 R 是秩 1 矩阵样本协方差等于单个快拍的贡献噪声没有被平均掉波束输出会叠加大量随机相位。解决保证 L 至少是阵元数的若干倍。常见做法是 L10M程序默认 L1001远大于 10*11所以默认参数没有这个问题。如果你要分析快拍数影响建议从 L110 到 L1000 之间取几个值观察收敛趋势。5.6 参考点选在阵列端点导致相位差很大但没有明显输出差异现象Nmid 从 1 改成 6峰值位置不变但数值上 bRb 与原来的差很小。原因参考点改变只是给方向矢量乘了一个公共相位因子二次型 bRb 在复数乘法中会抵消这个因子所以输出功率理论上不变。但在数值计算中端点参考会导致 Rn 的范围更大exp 的参数绝对值更大对浮点精度略不友好。解决仍然建议把参考点设在阵列中心比如 Nmid6这样声程差的正负分布更对称代码逻辑也更接近“阵列中心聚焦”的物理直觉。虽然结果一样但调试中间变量时看起来更干净。6. 进阶玩法把二维声图改成三维定位与多声源分辨到这里你已经有了一套能跑的二维切片波束形成工具。接下来我建议做三件事把程序的价值发挥出来。第一件事是把扫描扩展成三维。把 y 也变成扫描轴比如yscan [0:0.5:20]; % 粗扫描先取 0.5m 步长 Pcbf3D zeros(length(z), length(x), length(yscan)); for k11:length(z) for k21:length(x) for k31:length(yscan) Ri sqrt((x(k2)-xi).^2 (yscan(k3)-yi).^2 (z(k1)-zi).^2); Rn Ri - Ri(Nmid); b exp(-j*w*Rn/c); Pcbf3D(k1,k2,k3) abs(b*R*b); end end end三维网格是三维数组用 isosurface 画等值面能看到声源在空间中的分布。但计算量是二维的三次方量级所以步长要先取大。你会发现 y 方向分辨率其实很差因为线阵只在 x 方向有孔径对 y 轴的定位能力非常弱这就是为什么实际系统中要采用十字阵或面阵。第二件事是把单声源改成双声源。在生成 p1 的循环里叠加第二个声源的贡献同时保持采样频率和阵元坐标不变。你会看到两个主瓣出现在声图里但距离较近时会产生干涉条纹。这时候把方向矢量乘上窗函数比如b b .* hamming(M)旁瓣会下降两个声源的分辨也会变得更清晰。这就是窗函数在波束形成中最直观的作用比看理论曲线有意思得多。第三件事是固定随机种子做一组参数扫描记录。变量分别是 M5、11、21d1.5、3记录峰值误差、半功率宽度和最大旁瓣级。你会发现自己能总结出一张完整的规律表比如 M 增大让主瓣变窄d 超过半波长引入栅瓣旁瓣只随加权变化而不随 M 明显变化。这个表可以做进实验报告也能帮你判断某种场景下换算法值不值得。我从接触这段程序开始一直保留一个习惯每次跑波束形成仿真前先检查 d 和 λ 的比值再用 rng 固定随机种子最后用 ind2sub 而不是 find 提取峰值。这三个动作基本避开了所有让我翻过车的坑。这段程序虽然简单但它是理解全部自适应波束形成算法的基础改一改 R 的构造方式就是 MVDR改一改扫描几何就是近场聚焦改一改窗就是低旁瓣设计。希望这份常规波束形成的 MATLAB 程序能成为你继续往下走的起点也希望我踩过的这些坑能帮你少走几步弯路。本文还有配套的精品资源点击获取