
简介针对不同频率下宽度波束形成的仿真需求这份Matlab资源提供了完整的实现与结果展示适合信号处理、阵列信号处理方向的本科生、研究生及科研人员用于教学与仿真参考。包内共9个文件包含LCMV.m、LFMsource.m、fft_8_1.m等3个可运行脚本覆盖核心算法实现与频域处理流程5张png图直观呈现不同频率下的波束图和仿真结果便于对照分析另有说明txt梳理运行方法和注意事项即使不熟悉Matlab环境也能按说明快速上手。整体包体仅518KB轻量易用。目前已有114人学习浏览内容围绕波束形成、LCMV准则、LFM信号源等典型知识点展开代码结构清晰、注释易读适合在课程设计、毕业设计或课题预研中直接参考和二次开发也可作为理解频域波束形成特性的入门实例。1. 同一阵列在不同频率下的波束宽度为什么会差出好几倍一套 8 元等间距线阵工作频率从 1 kHz 抬到 4 kHz主瓣半功率宽度能从约 19 度缩窄到不到 5 度。这不是阵列坏了而是阵元间距与波长的比值 d/λ 变了阵列在电尺寸上变大了。拿到这个压缩包时我第一反应就是先看 fft_8_1.m 和 LCMV.m 这两个脚本因为它们恰恰对应了波束形成里最常用的两条路频域常规波束形成和自适应约束波束形成。再加上 LFMsource.m 负责生成线性调频源三个文件拼起来就是一条完整的仿真链路。本文按原理 → 频域实现 → 自适应实现 → 频率调参验证的顺序拆解这套代码新手可以直接对照运行熟手可以重点关注频率切换时哪些参数必须同步调整。2. 波束形成的数学基础与 MATLAB 仿真信号链2.1 阵列导向矢量波束形成的“坐标系”波束形成的本质是给每个阵元信号乘一个复权重让目标方向的信号同相叠加其他方向因相位不一致而相互抵消。对均匀线阵假设阵元间距为 d信号从 θ 方向入射第 n 个阵元相对参考阵元的时延为τn n·d·sinθ / c在频域这个时延就变成相移 e^(-j2πf·τn)。把所有阵元的相移排成列向量就是导向矢量 a(f, θ)a(f, θ) [1, e^(-j2πf·d·sinθ/c), …, e^(-j2πf·(N-1)·d·sinθ/c)]^T代码里的频率 f、阵元数 N、间距 d、声速 c 都在这个公式里。需要注意声速取值压缩包说明里如果按水下场景处理c 一般取 1500 m/s按空气声学处理则取 343 m/s。这直接决定波束指向是否正确是仿真结果对不上的第一排查点。2.2 频率、阵元间距与主瓣宽度的定量关系均匀线阵宽边方向的半功率波束宽度近似为BW ≈ 0.886·λ / (N·d)λ c/f 代入后BW ≈ 0.886·c / (N·d·f)。也就是说在阵元数和间距固定的前提下频率越高波束越窄。下表给出了 c1500 m/s、N8、d0.5 m 时的理论值这个表可以拿来做仿真结果的对照基准频率 f (Hz)波长 λ (m)d/λ半功率波束宽度 BW度10001.5000.333约 19.020000.7500.667约 9.540000.3751.333约 4.880000.18752.667约 2.4注意 d/λ 0.5 时会出现栅瓣风险仿真时如果扫到多个等高的峰值先别怀疑代码检查一下是不是阵元间距已经超过半波长。很多 MATLAB 教程不会强调这个边界条件但在实际调参时它比算法本身更容易出问题。2.3 LFMsource.m 生成了什么信号LFMsource.m 生成的是线性调频信号也叫 chirp 信号。它的瞬时频率随时间线性变化数学形式为s(t) A·exp(j·2π·(f0·t K·t²/2))其中 K 是调频斜率。线性调频信号在雷达和声呐里常用是因为它带宽大、经匹配滤波后能获得很高的距离分辨力。下面是 LFMsource.m 这类生成函数的典型实现框架function s LFMsource(N, fs, f0, f1) % N: 采样点数, fs: 采样率, f0: 起始频率, f1: 截止频率 t (0:N-1) / fs; % 时间轴 K (f1 - f0) / (N / fs); % 调频斜率 s exp(1j * 2 * pi * (f0 * t 0.5 * K * t.^2)); end这段代码先构造时间轴再计算调频斜率最后生成复解析信号。用复信号而不是实信号的原因波束形成后续要做相位补偿复信号可以直接在指数项上操作避免实信号频谱正负半轴叠加带来的麻烦。压缩包里的 LFMsource.m 可能还包含加窗、幅度加权之类的内容但核心就这四行。如果你在别的项目里复用记得检查 fs 是否满足奈奎斯特条件否则线性调频信号会混叠仿真出来的波束图主瓣位置都会偏移。2.4 MATLAB 版本差异与脚本兼容性这套代码在 MATLAB 2014a、2019a、2021a 上都能跑原因在于它只用了 FFT、矩阵乘、eig/svd 这类基础函数没有依赖较新的工具箱。常见的运行问题出在路径和脚本名上直接把脚本拖进编辑器按 F5 运行有时会因为当前文件夹不在搜索路径里报未定义函数解决方法是右键文件夹选择添加到路径或者用 cd 命令切到脚本所在目录。另一个跨版本差异在绘图样式上2019a 之后默认的 colororder 变了同一套 plot 代码画出来的线条颜色不同这不影响数据正确性但如果你要把仿真图放到论文里建议统一用 MATLAB 2021a 出图。3. fft_8_1.m频域 FFT 波束形成的算法拆解3.1 数据流从时域阵元信号到频域快拍fft_8_1.m 代表的是频域常规波束形成。它的输入是各阵元的时域采样信号输出是不同扫描角度下的波束输出功率。整体流程分三步先对各阵元信号做 FFT 变换到频域再对每个频点乘上对应的相位补偿因子最后对各阵元做相干求和并统计能量。这比时域延时求和快得多因为一次 FFT 就把整个频带都处理完了后面每个频点是纯复数乘法。接收信号的构造在仿真里通常是这样的假设有一个从 θ0 方向来的线性调频信号 s(t)那么第 n 个阵元收到的信号是 s(t - τn)在频域表示为Xn(f) S(f)·e^(-j2πf·n·d·sinθ0/c) Nn(f)这里 Nn(f) 是噪声的频域表示。代码里通常会用一个循环生成各阵元的数据再按列做 FFT。我习惯把阵元放在行方向、时间放在列方向这样每行是一个阵元的一帧数据fft 之后每行是这个阵元的频谱。% 构造 8 阵元接收数据示意写法与 fft_8_1.m 逻辑对应 N 8; d 0.5; c 1500; fs 10000; theta0 30; % 目标入射角 s LFMsource(1024, fs); % 线性调频源信号 S fft(s); % 源信号频谱 % 生成各阵元频域接收信号 X zeros(N, length(S)); for n 0:N-1 tau n * d * sind(theta0) / c; % 各阵元时延 X(n1, :) S .* exp(-1j * 2 * pi * (0:length(S)-1) * fs/length(S) * tau); end这个循环里最关键的是 exp 参数里的频率项。很多人直接把 FFT 的下标当作频率却忘了乘上频率分辨率 fs/Nfft导致相移完全错乱。注意这里 (0:length(S)-1) * fs/length(S) 才是每个 FFT bin 对应的实际物理频率。3.2 指向性补偿与波束输出功率拿到频域快拍后下一步是对每个扫描角度 θ 构造导向矢量并做相位共轭相乘再求和这就是数字波束形成里的延时-求和操作在频域的等价形式。对第 k 个频点、第 i 个扫描角输出为y(θi, fk) w^H(θi, fk) · X(fk)其中 w(θi, fk) a(θi, fk) 本身就是导向矢量因为常规波束形成的最优权就是导向矢量本身。扫描的粒度一般取 0.1 度或 0.5 度太粗会错过主瓣峰值太细增加计算量但精度不会显著提升。% 波束扫描主循环 scan_theta -90:0.5:90; % 扫描角度范围 freq_bin 50; % 选取感兴趣的频点 bin f_k freq_bin * fs / length(S); % 该 bin 对应的物理频率 P zeros(size(scan_theta)); for i 1:length(scan_theta) theta scan_theta(i); a exp(-1j * 2 * pi * f_k * d * (0:N-1). * sind(theta) / c); P(i) abs(a * X(:, freq_bin))^2; % 该方向输出功率 end P_dB 10 * log10(P / max(P)); % 归一化并转 dB这段代码先把扫描角定义好再取一个固定的频点做功率统计。P(i) 的计算用了共轭转置 a这是频域波束形成的标准写法。最终得到的 P_dB 就是你要画的波束方向图0 dB 处对应主瓣指向。如果不取单一频点而是对频带内多个频点的功率做平均得到的就是宽带波束输出但那要额外处理各频点的指向一致性fft_8_1.m 之所以叫这个名字我理解就是固定频点下的 8 元 FFT 波束形成。3.3 参数速查表与运行步骤参数名在代码中的位置推荐初值调整说明N阵元数脚本开头8增大使波束变窄计算量线性增加d阵元间距脚本开头0.5 m必须小于 λ/2否则出现栅瓣c声速脚本开头1500水下默认空气声学改 343fs采样率生成信号处10000需大于信号最高频率的 2 倍scan_theta扫描循环-90:0.5:90栅瓣观察需要扫到 ±90 度freq_bin主循环50对应频率 488 Hz按需修改运行步骤很简单把压缩包解压到无中文路径的目录打开 fft_8_1.m直接运行看生成的 1.png 或直接看 Figure 窗口。建议先别改任何参数把原始结果跑出来再改 freq_bin 观察波束宽度变化这样能最快建立频率-宽度的直观认识。如果运行报错看错误提示是不是函数未定义优先检查 2.4 节说的路径问题。4. LCMV.m线性约束最小方差波束形成4.1 为什么常规波束形成在强干扰下会失效fft_8_1.m 的问题在于它永远使用固定的导向矢量作为权值不管其他方向有没有强干扰。当干扰功率比目标信号高 20 dB 时常规波束形成的旁瓣泄漏会把目标完全淹没主瓣里看不出任何目标峰值。LCMV 的思路是在保证目标方向增益不变的约束下让输出功率最小化从而自动在干扰方向形成零陷。LCMV 的数学优化问题写成min w^H·R·w约束条件 C^H·w f其中 R 是接收数据的协方差矩阵C 是约束矩阵f 是约束响应向量。目标方向约束为增益 1也就是 a^H(θ0)·w 1。它的闭式解为w_opt R^(-1)·C·(C^H·R^(-1)·C)^(-1)·f这个式子看起来复杂但代码实现只有三行左右难点反而在协方差矩阵估计是否准确。4.2 约束矩阵与最优权向量的求解LCMV.m 里的核心代码框架大致如下% R: 协方差矩阵, C: 约束矩阵, f: 约束响应 K 200; % 快拍数 R zeros(N, N); for k 1:K xk X_noisy(:, k); % 第 k 个快拍 R R xk * xk; % 累加外积 end R R / K; % 求平均得到协方差估计 C a_theta0(:); % 约束矩阵只取目标方向导向矢量 f 1; % 约束增益为 1 w_lcmv (R \ C) / (C * (R \ C)) * f; % 矩阵求解得到最优权 % 用最优权做波束扫描 P_lcmv abs(w_lcmv * a_scan).^2;注意代码里用了反斜杠运算符 R \ C 而不是 inv(R) * C。R \ C 本质上是高斯消元求解数值稳定性远好于显式求逆在矩阵接近奇异时优势尤其明显。LCMV 的波束图会明显区别于 fft_8_1.m 的结果目标方向增益保持为 1但会在干扰方向出现一个很深的凹陷凹陷深度可以到 -40 dB 以下。这里有个容易犯的错约束矩阵 C 不能只放目标方向导向矢量如果你知道干扰的方向并且想强制零陷需要把干扰方向导向矢量也放进 C 里并且把对应的希望增益设为 0。这就是线性约束四个字的含义——你给算法规定了必须满足的增益条件。4.3 快拍数不足与仿真发散的处理LCMV 仿真最常见的失败现象是权向量计算结果巨大、波束图乱振甚至仿真发散。原因通常是快拍数 K 太少协方差矩阵 R 病态甚至奇异。信号维数是 8理论上 K 至少大于 8但实践里 K 100 时波束图都会很不稳定。我一般会在 LCMV.m 源码里加一行对角加载R_loaded R 0.01 * eye(N); % 对角加载加载量取 R 对角线均值的 1% w_lcmv (R_loaded \ C) / (C * (R_loaded \ C)) * f;对角加载的本质是给协方差矩阵的对角线加一个小常数抬高它的最小特征值避免求逆时除以接近零的数。加载量的取值没有黄金标准一般从 R 对角线均值的 0.1% 试到 10%观察波束图是否稳定。加了加载之后零陷深度会稍浅但换来的是波束形状的稳定这对教研场景来说更重要。如果你在仿真结果里看到 LCMV 的输出功率在某些角度跳变剧烈优先怀疑的对象就是 R 条件数而不是约束矩阵写错。下面用一个表对比两个脚本的特性方便你决定什么时候用哪个对比项fft_8_1.mCBFLCMV.m计算量每频点复杂度低约 O(N log N)需要估计协方差并求解约 O(N^3)抗干扰能力无自适应机制旁瓣固定干扰方向自动形成零陷阵元误差敏感度低高幅度相位误差会破坏零陷对快拍数的依赖不依赖依赖快拍少时需对角加载频率变化的影响波束宽度随频率变化约束矩阵需按频率更新否则指向偏实际项目里两者常常配合使用先用 fft_8_1.m 快速扫一个粗测角度再在该角度附近用 LCMV 做精细处理和干扰抑制。5. 用波束图验证频率-宽度关系及调参技巧5.1 半功率波束宽度的实测方法仿真跑完后别只盯着一幅图看用一个脚本把主瓣的半功率宽度自动量出来这才是能复用到其他工程里的本事。测量逻辑是先找峰值再从峰向两侧找第一个降到 -3 dB 的点两点之间的角度差就是半功率波束宽度[peak, idx_peak] max(P_dB); % 找主瓣峰 left find(P_dB(1:idx_peak-1) -3, 1, last); % 左边界 right idx_peak find(P_dB(idx_peak1:end) -3, 1, first); % 右边界 bw scan_theta(right) - scan_theta(left); % 半功率宽度这套逻辑假设波束图只有一个主瓣如果出现栅瓣测出来的宽度会异常大。所以每次测量前先检查主瓣两侧是否存在等高的第二峰也就是栅瓣。栅瓣的角位置可以通过 arcsin(kλ/d) 估算简单做法就是把 d 改小到 λ/2 以下对比波束图的变化。5.2 频率切换后需要同时改哪些参数用这套代码仿真正是在不同频点下对比时有一个高频踩坑点直接改 fft_8_1.m 里的 freq_bin 就去跑把频率从 1 kHz 调到 4 kHz但 d 还停留在 0.5 m这时 d/λ 1.33已经大于 0.5会在 ±48 度附近出现栅瓣看起来像多了一个目标。正确的调参顺序是先根据目标频率算出 λ c/f再检查 d确保 d ≤ λ/2最后再跑仿真。相比之下LCMV 脚本还要多一步重新计算目标方向的导向矢量并更新约束矩阵 C否则约束方向会和真实目标方向错开自适应权会把目标当成干扰抑制掉。5.3 零陷位置的快速校验做 LCMV 仿真验收时我会用一个笨但有效的办法强行在协方差矩阵里注入一个已知方向的干扰信号然后看波束图的零陷是否落在该方向。注入干扰的代码就一行R R sigma_j^2 * (a_j * a_j); % a_j 为干扰方向导向矢量如果零陷位置偏差超过 1 度说明声速值设错或频率参数没对齐优先回到 2.1 节核对导向矢量公式。这个方法也可以反着用扫描零陷位置反推实际频率在有频率偏移的被动探测场景下这本身就是一种频率估计算法。把这两步跑通这份压缩包的仿真结果就不只是几张图而是一套可以继续扩展的频率-波束联合分析工具。本文还有配套的精品资源点击获取