简介面向无线通信、雷达及声纳阵列信号处理学习者的 MATLAB 波束赋形专题文档系统讲解均匀线阵方向图绘制、波束宽度与波达方向及阵元数的关系、栅瓣产生与抑制、最优权傅里叶变换、最大信噪比准则方向图与功率谱、ASC旁瓣相消MSE准则等核心内容。资源包含8/16/128/1024阵元等多种配置的仿真对比直观展示阵元数对波束宽度、分辨力及旁瓣性能的影响并配有完整可运行的MATLAB代码、仿真图与文字说明适合从入门到进阶的开发者边看边练也可支撑课程设计或科研预研。资源包为单个doc文档大小1.23MB结构紧凑便于快速定位学习。已有335人浏览学习/下载。文档中多段示例代码可直接复制运行例如8阵元均匀线阵方向图、最大信噪比准则下的方向图与功率谱、旁瓣相消MSE准则等能帮助读者掌握波束赋形权值计算、信干噪比优化与空间干扰抑制方法加深对空域滤波技术的理解。1. 从一份 doc 开始的波束形成空域滤波为什么值得自己写一遍手机上收音的降噪、雷达从杂波里捞目标、声呐在浅海分辨潜艇底层是同一件事从多路信号的相位差里反推空间的来波方向再把阵列响应按方向加权。matlab beamforming 这个组合之所以流传广是因为 MATLAB 把矩阵运算和信号处理函数都备好了缺的只是从花哨公式到能跑通仿真的那一层工程封装。很多人手里都流转过一份叫 matlab beamforming.doc 的旧文档概念写得很全但要让它变成真正能给出波束图、能抗干扰、能处理宽带语音的代码中间还得补齐导向矢量、协方差估计、稳健性和宽带处理四块。下面按我平时搭仿真验证的顺序来推进先写相移波束形成再做 MVDR再拆宽带最后讲一套通用的验证和校准套路。2. 波束形成的数学内核导向矢量、波束图与相移波束形成2.1 导向矢量波束形成的“单位冲激响应”在远场平面波假设下阵列各阵元收到的信号只差一个由到达角和阵元位置决定的相位。均匀线阵ULA的导向矢量写成function a steering_vec(theta_deg, M, d_lambda) % 生成均匀线阵(ULA)的阵列流形矩阵 % theta_deg: 来波方向单位度可传行向量 % M : 阵元数量 % d_lambda : 阵元间距单位用波长归一化 theta theta_deg(:).; % 相邻阵元相位差 2*pi*d*sin(theta)/lambda phase 2 * pi * d_lambda * sind(theta); a exp(1i * phase .* (0:M-1).); endd_lambda是这段代码里最需要留意的参数。半波长间距d_lambda 0.5是默认选择比它大扫描范围两端会出现栅瓣比它小阵列物理孔径变小主瓣变宽角度分辨率下降。sind在 MATLAB 里接收的是角度而非弧度这点和sin容易混导致相位算错后波束图完全不对称。实际里我习惯把它和deg2rad换算成套件任何改动都只发生在这一层封装里后面所有算法都调用同一个函数避免一处一处改角度单位。如果把阵列响应类比成 FIR 滤波器的冲激响应来波方向对应频率扫描角度对应频率扫描波束形成本质上就是空域滤波。相移补偿就是旋转因子这和 DFT 里的exp(jwn)是同构的只差一个空间采样位置换时间采样位置。2.2 相移波束形成延迟-求和的最小可运行实现先写一个最朴素的相移波束形成器把阵列指向 30 度方向视角对准均匀加权w a0 / M然后扫-90 ~ 90度画出波束图。M 16; d_lambda 0.5; theta_scan -90:0.1:90; % 目标来波方向 a0 steering_vec(30, M, d_lambda); w a0 / M; % 均匀加权每路幅度相同只做相位补偿 A steering_vec(theta_scan, M, d_lambda); P 20 * log10(abs(w * A)); P P - max(P); % 主瓣方向归一化为 0 dB plot(theta_scan, P, LineWidth, 1.5); grid on; xlabel(方位角 / deg); ylabel(归一化功率 / dB); ylim([-40 5]);代码只做了三件事构造导向矢量、用共轭转置做匹配滤波、归一化后画图。w * A展开后就是每个扫描角度上的匹配响应w要先共轭是因为导向矢量是复数匹配滤波要求相位对齐而不是单纯内积。扫描步长取0.1度对波束图足够平滑如果要定位峰值或做高精度 DOA 估计再把步长缩到0.01度否则没必要增大计算量。波束形成在时域上看就是延迟求和相位补偿等价于把每个阵元的接收信号对齐到同一波前。窄带假设下延迟可以退化成复指数相乘所以 MATLAB 里做窄带波束形成完全不用处理分数延迟这是它能用这么短代码跑起来的原因。2.3 波束图怎么读主瓣、旁瓣和孔径的关系均匀加权 ULA 的第一旁瓣高度固定在-13.26 dB附近这是矩形窗的空域版本和频谱泄漏是同一个数学问题。天线阵列的孔径越大波束越窄半功率波束宽度近似为BW_3dB ≈ 0.886 * lambda / (M * d * cos(theta0))theta0是波束指向角。用这个近似估算 10 元和 20 元阵列的分辨能力比每次跑仿真快得多。阵元数 M间距指向3dB 主瓣宽度近似第一旁瓣80.5 λ0°12.7°-13.3 dB160.5 λ0°6.4°-13.3 dB320.5 λ0°3.2°-13.3 dB旁瓣太高会让强干扰从旁瓣漏进来均匀加权并不是实际工程的第一选择。加 Hamming 窗或 Taylor 窗可以把旁瓣压到-40 dB以下代价是主瓣展宽约一半。MATLAB 里直接用w win * a0 / M就行这里win是和阵元数等长的窗向量相当于把 FIR 设计里的加窗法平移到了空域。真正到了干扰抑制需求明确的场景固定窗函数就不够用了下面进入自适应。3. 自适应波束形成MVDR 的 MATLAB 数值实现与稳健性参数3.1 为什么要自适应固定波束形成扛不住强干扰相移波束形成的权重只由a0决定来波里混进一个功率高 20 dB 的干扰时即便它在旁瓣方向泄漏功率也足以把目标信号淹没。波束形成的优势体现在阵列增益10lg(M)也就是均匀加权下信噪比最多提升M倍但这是对白噪声而言。对有色干扰需要根据接收数据的协方差矩阵调节权重把零陷对准干扰方向。MVDR最小方差无失真响应的约束写成min w * R * w s.t. w * a0 1目标方向增益固定为 1同时最小化输出总功率。闭式解是w inv(R) * a0 / (a0 * inv(R) * a0)R 是干扰加噪声的协方差矩阵。在仿真里R通常由快拍数据估计直接inv(R)不是好习惯矩阵条件数稍大结果就不可信应当用R \ a0解线性方程。MVDR 的方向图会在干扰处自动形成零陷零陷深度不受M直接限制而受 R 估计精度和数据平稳性限制。3.2 MVDR 实现代码与对角加载参数function w mvdr_weights(Rxx, a0, delta) % MVDR 权重求解带对角加载 % Rxx : M x M 协方差矩阵由采样快拍估计 % a0 : 目标方向导向矢量 % delta: 对角加载系数标量 Rl Rxx delta * eye(size(Rxx)); w Rl \ a0; w w / (a0 * w); % 无失真约束归一化 end调用端还有一个容易被忽略的用法对角加载系数的量纲要和 R 的量纲一致否则加载不加载没有工程意义。N 2000; % 快拍数 X ...; % M x N 基带复采样数据 R (X * X) / N; % 采样协方差估计 target_a steering_vec(30, M, d_lambda); delta 1e-3 * trace(R) / M; % 相对功率的对角加载 w mvdr_weights(R, target_a, delta);对角加载的本质是在小特征值方向上注入白噪声压低协方差矩阵的动态范围防止 ”噪声子空间“ 被估计误差抬起来。参数调不好会出两个极端加载太小零陷深但极不稳定加载太大零陷被抹平退化成常规波束形成。加载系数相对 trace(R)/M典型效果适用场景0零陷极深信号自消风险高无限精度、快拍充足1e-6 ~ 1e-4零陷深数值敏感高精度仿真、快拍 10M1e-3 ~ 1e-2零陷 -20 ~ -40 dB稳健实测数据、快拍较少 1e-1接近常规波束形成强失配、低信噪比快拍数少于2M时采样协方差矩阵的特征值开始散开大特征值变大、小特征值变小MVDR 的性能急剧下降。这时候除了对角加载还可以用前后向平滑或 Toeplitz 化处理后者的做法是把 R 沿反对角线取平均写起来只有几行代码效果却能抵御阵元间的相位误差。3.3 导向矢量失配时的救法WNC 和 LCMV 推广MVDR 最大的工程坑不是矩阵求逆而是a0不准。阵元位置标定误差、互耦、通道幅度相位不一致都会让真实导向矢量和理论值偏离MVDR 会把这个失配当成目标信号的一部分进行相消结果就是目标也被抑制输出 SNR 反而比常规波束形成更差。这种现象叫信号自消在低快拍下尤其明显。一个实用的兜底方案是用白噪声增益约束WNC替换对角加载约束权重向量的范数限制白噪声增益不低于某个阈值例如-10 dB。实现时把加载系数放在导向矢量一侧做迭代求解替代直接对 R 加对角。LCMV线性约束最小方差则是把单个无失真约束扩展成多个线性约束适合同时需要零点约束和方向图保形的场景约束矩阵 C 和目标响应向量 f 一进去解就是w R \ C * inv(C * R \ C) * f这一套在 MATLAB 里用矩阵左除几行就能写但参数的物理含义必须先说清楚约束越多自由度消耗越大能抑制的独立干扰数目越少。ULA 的可用自由度是M-1一个约束点吃掉一个自由度这是设计时先要算清楚的账。4. 宽带波束形成从窄带假设失效到频域/时域两套实现4.1 什么时候窄带会失效麦克风阵列的典型带宽相移波束形成建立在窄带假设上整个信号带宽内导向矢量基本不变。判断标准是阵列渡越时间要远小于信号相关时间也就是D c / BD 是阵列孔径B 是信号带宽c 是波速。对 8 kHz 采样带宽的语音c / B 343 / 8000 ≈ 4.3 cm一个 16 元、间距 4 cm 的麦克风阵列孔径就有 60 cm完全不满足窄带条件。雷达脉冲时宽大、相对带宽小窄带假设通常成立语音和声呐就得老老实实走宽带波束形成。宽带处理有两条主流路径频域做法是把数据切帧做 STFT在每个频率子带里独立做窄带波束形成最后合回时域时域做法是给每个通道接一串 FIR 抽头用约束最小二乘或自适应算法迭代出抽头系数。频域方法调试方便代价是块延迟时域方法适合实时系统但滤波器阶数选择和数值稳定性都要花功夫。4.2 频域实现STFT 分帧、子带 MVDR 与合成频域宽带波束形成的骨架分三步分帧加窗做 FFT逐频率点估计协方差并求解权重重叠相加合成时域。% 参数fs16000帧长512hop256M8d0.04m Nfft 512; hop 256; win hann(Nfft, periodic); frames buffer(x, Nfft, Nfft-hop, nodelay); X fft(bsxfun(times, frames, win), Nfft, 1); % 单通道分帧 % 多通道按同样方式分帧得到 X_ch: Nfft x n_frames x M Y zeros(size(X)); % 从第2个bin到Nyquist bin每个频率独立处理 for k 2:Nfft/21 f_k fs * (k-1) / Nfft; d_lambda d / (343 / f_k); % 频率升高等效阵元间距变大 a_k steering_vec(theta0, M, d_lambda); Xk squeeze(X_ch(k, :, :)).; % M x n_frames Rk Xk * Xk / size(Xk, 2); wk Rk \ a_k; wk wk / (a_k * wk); % 每个频率单独保证无失真 Y(k, :) wk * Xk; end % 再对 Y 做重叠相加 IFFT 得到时域输出这段代码里的关键陷阱是d_lambda必须随频率变化。低频段d / lambda很小阵列电尺寸小波束很宽方向性弱高频段电尺寸变大可能出现空间混叠所以要限制处理频带上限。实际参数调试时我一般固定f_low 300 Hz、f_high fs/2低于300 Hz的 bin 直接不处理因为低频段的阵列增益实在有限。参数推荐值说明采样率16 kHz语音常用最高 bin 8 kHz帧长512 (32 ms)频率分辨率约 31 Hz帧移25650% 重叠合成时抑制窗效应窗型周期性 Hann主旁瓣均衡适合重叠相加处理频带300 Hz ~ 7.5 kHz避开低频噪声和 ds混叠频域方法还有个数值上的好处每个 bin 上的快拍数等于帧数16 kHz、10 秒语音约有 624 帧远大于2M 16协方差估计充裕。要注意的是相邻帧的窗重叠会造成数据相关性等效快拍数要打个折扣但工程上影响不大。4.3 时域 FIR 波束形成抽头延迟线和多通道维纳解时域宽带波束形成的思路是把每个通道扩展成T个抽头形成M*T维的输入向量再用 LCMV 或维纳解求权重。抽头延迟线本质上就是一个稀疏的 FIR 滤波器组目标方向的无失真条件和干扰方向零陷都写在约束矩阵里。T 32; % 每通道抽头数 % 构造扩展数据X_ext(m*t t_index, n) % 用 toeplitz 为每个通道生成延迟矩阵再纵向拼接 % 约束目标方向所有抽头之和为1其余线性约束为0 A_ext kron(steering_vec(theta0, M, 1), ones(T, 1)); % 最小二乘解可以直接用 MVDR 的闭式解 % 只是 R 变成 M*T 维的扩展协方差矩阵 R_ext X_ext * X_ext / length(x); w_ext R_ext \ A_ext; w_ext w_ext / (A_ext * w_ext);T的选取和信号带宽成正比抽头太少频率响应在带内起伏抽头太多协方差矩阵维度变大需要更多快拍数据。一般按T ≈ 10 * fs / f_center起步再用仿真对比输出 SNR 收敛曲线。MATLAB 优化工具箱里的fmincon也可以用来做抽头系数的非线性约束优化但多数场景下线性约束闭式解已经够用不必把问题复杂化。频域和时域方案的工作量差异主要体现在调试频域方法可以单独看某个频率 bin 的波束图定位是哪一个频段出了问题时域方法只能看到整体频率响应。但时域方法延迟小得多适合回音消除这种对延迟敏感的场景。做仿真验证时我建议先走频域确认算法逻辑正确后再移植到时域。5. 波束形成仿真验证信号注入、失配校准与空间谱定位5.1 模拟多源阵元数据信噪比、分数延迟与相干源陷阱仿真验证的第一步是用已知信号构造阵列数据。注意分数延迟不能简单取整高频段会有明显相位误差。常见做法是先对信号做 8 倍上采样移动整数个采样点后再抽取回原速率fs 16000; c 343; M 8; d 0.04; t (0:fs*10-1)/fs; s_target sin(2*pi*1000*t); s_interf sqrt(10^(20/10)) * sin(2*pi*1500*t); % 比目标高20dB up 8; s_up resample(s_target, up, 1); delay_up round((0:M-1) * d * sind(30) / c * fs * up); X zeros(M, length(s_target)); for m 1:M X(m, :) resample(s_up, 1, up); % 用上采样序列移动整数延迟 end X X awgn(X, 10, measured); % 加噪SNR10dB要验证 MVDR 的干扰抑制能力用独立噪声源生成干扰就够。但要测空间平滑或去相关算法就必须注意相干源问题当两个信号完全相干同一信号的不同延迟协方差矩阵会出现奇异MVDR 会把它们当成一个源处理。这时候要么在信号里加微小的随机扰动打破相干性要么专门做前后向平滑验证算法。5.2 画空间谱和用 MUSIC 验证 DOA 定位能力验证权重算得对不对除了看波束图还要看空间谱。MVDR 的伪谱直接利用协方差矩阵峰值位置对应来波方向theta -90:0.5:90; A steering_vec(theta, M, 0.5); P_mvdr zeros(size(theta)); for i 1:length(theta) a A(:, i); P_mvdr(i) 1 / abs(a * (R \ a)); end P_mvdr 10 * log10(P_mvdr / max(P_mvdr)); plot(theta, P_mvdr, LineWidth, 1.5); xlabel(方位角 / deg); ylabel(归一化空间谱 / dB); grid on;MVDR 谱峰越尖锐说明 R 估计越稳出现多个毛刺说明快拍不足或加载系数偏小。MUSIC 走的是特征分解路线把协方差矩阵分成信号子空间和噪声子空间用噪声子空间和导向矢量的正交性找峰分辨率比 MVDR 谱更高但要求信源数已知。做阵列信号处理时两者搭配使用MUSIC 确认来波方向MVDR 验证该方向上的无失真增益。图存下来时用exportgraphics(gcf, beam_pattern_16.png, Resolution, 300)比print这个老接口在 MATLAB 2022b 之后的版本里更稳导出论文用的 EPS 也是一样的路径。5.3 失配自检三板斧位置误差、幅度相位误差、校准仿真通过不等于外场能用我会在交付前做三种失配注入测试每个都模拟真实的工程坑。第一板斧是阵元位置误差给每个阵元的位置乘(1 0.01 * randn)也就是 1% 的间距扰动。常规波束形成在这种误差下几乎无感但 MVDR 零陷深度可能损失 10 dB 以上这是判别代码里有没有对角加载的最直接方法。第二板斧是通道幅度相位误差幅度加2%的随机起伏、相位加2°抖动测试方式是画 50 次蒙特卡洛的输出 SINR 分布。理想情况下 MVDR 的 SINR 方差明显大于常规波束形成如果方差太大就要把加载系数往上调一档。第三板斧是校准。用已知方向的扬声器或喇叭在消声环境采集阵列响应把每个通道的复数增益存成calibration_16k.mat实际数据处理时先做X X ./ cal_vector再估计协方差。校准文件按频率点分开存频域波束形成的每个 bin 乘各自校准系数。这套流程从仿真平滑过渡到外场基本不会出现仿真很好、实测全瞎的断崖。本文还有配套的精品资源点击获取