波束形成本质是空间滤波：线性积分视角与MATLAB实现

做阵列信号处理这些年，经常被问到同一个问题：波束形成到底该怎么入门。我的答案通常比较直接——先别急着翻自适应算法的书，先把“空间滤波”这个本质吃透。波束形成输出的数学本质，就是对阵列各阵元的接收信号做一次线性积分加权求和，这个权向量怎么设计，决定了你的“空间滤波器”长成什么样。这套思路一旦打通，再回头看延迟求和、MVDR、超指向性波束形成，会发现它们都是同一个积分架构下换了一套加权策略而已。这篇文章就从线性积分这个底层视角出发，把波束形成的原理完整过一遍，再给出可以直接运行的MATLAB实现，覆盖阵列接收信号建模、导向矢量计算、方向图绘制、目标信号提取，以及宽带信号的频域处理方法。内容适配正在学习阵列信号处理的学生，也适配需要把麦克风阵列技术落到实际项目里的工程师。代码只用标准MATLAB，不需要任何商业工具箱。

1. 从“线性积分”理解波束形成的底层逻辑

1.1 麦克风阵列为什么能做到空间指向

先把一个反直觉的事实摆出来：单个麦克风是全向的，它只能记录声音的强弱，不知道声音从哪里来。但当你把多个麦克风按固定间隔排成一排，情况就完全不同了——同一个声源发出的波，到达不同阵元的时刻会有微小差异。这个时间差在频域里就表现为相位差。远场平面波从角度θ入射时，相邻阵元之间的距离d在波传播方向上投影为d·sinθ，所以相邻阵元之间的时延是d·sinθ/c，对应的相位差是2πf·d·sinθ/c。

波束形成干的事，就是对每个阵元的信号乘上一个复数权重，再叠加起来。某个方向来的信号，如果它的相位差刚好和权重里的相位补偿匹配，就会同相叠加，幅度变大；不匹配的方向，各路信号相位错乱，相互抵消。整个过程就相当于在设计一个“只让特定方向声音通过”的空间滤波器。注意，这里的“滤波器”三个字很关键：时间域滤波器是对时间信号做卷积积分，空间域波束形成是对空间各点的信号做加权积分，两者在数学上完全同构。

1.2 连续孔径的积分表达

把问题先放到连续域来看，数学就更干净了。假设阵列从x = −L/2到x = L/2连续分布，孔径函数为w(x)，它表示在位置x处的加权系数。远场平面波入射时，空间中位置x处的复声压在频域可以写成p(x, f)，那么阵列输出对全孔径做一个积分：

Y(f) = ∫ w(x) · p(x, f) dx

这就是整个波束形成的“总表达式”。如果入射方向为θ，参考位置在原点，那么p(x, f) ∝ exp(−j·k·x·sinθ)，其中k = 2πf/c是波数。现在假设我们想让波束指向θ_s方向，也就是对θ_s方向入射的信号做相位补偿，那么孔径函数可以写成w(x) · exp(j·k·x·sinθ_s)，输出响应变成：

B(θ) = ∫ w(x) · exp(−j·k·x·sinθ) · exp(j·k·x·sinθ_s) dx = ∫ w(x) · exp[−j·k·x·(sinθ − sinθ_s)] dx

这个式子说明，波束方向图本质上就是孔径加权函数w(x)的傅里叶变换，只不过自变量从空间频率换成了方向角。这就是为什么说波束形成是线性积分——它的核心操作就是在空间维度上做积分变换。所有看似高级的波束形成算法，最后都落脚在“选择什么样的w(x)”这个问题上。

1.3 从连续积分到离散求和

实际系统里没有连续孔径，只有M个离散阵元。数学上可以认为在阵元位置x_m处放一个狄拉克采样函数，把连续积分退化成加权求和：

Y(f) = Σ_{m=1}^{M} w_m · p(x_m, f)

这就是代码里最常见的那一行：y = w' * X。权向量w就是连续孔径函数w(x)在离散点上的采样值。理解这一点很有用——比如延迟求和波束形成里那个“除以M”，本质上就是连续积分中1/L做归一化的离散版本，对应矩形积分公式的数值积分系数。

离散采样引出一个重要约束：空间采样定理。连续波场在波数域是有限带宽的，最高空间频率是k = 2π/λ。为了保证不同方向的入射波在阵元采样后还能区分开，阵元间距必须满足d ≤ λ/2。这和时域采样率必须大于两倍信号最高频率是同一个道理，只是把“时间”换成了“空间”。违背这个条件时，方向图会周期性重复，出现与主瓣等高的“栅瓣”，这一点后面专门展开讲。

2. 波束形成的核心原理和几种关键加权策略

2.1 远场平面波假设与适用范围

前面推导全部建立在远场平面波假设上。实际中麦克风阵列放在房间里，声源往往离得很近，什么时候能用平面波近似，要有个数。工程上常用瑞利距离来判断：r_far > 2·D²/λ，其中D是阵列总孔径。举个例子：M = 8个阵元，阵元间距d = 4.3cm（4kHz信号的半波长），阵列总孔径D = (M−1)·d ≈ 0.3m，波长λ ≈ 0.086m。算一下临界距离：2 × 0.09 / 0.086 ≈ 2.1米。也就是说，测试距离在2米以上时，平面波假设成立；距离近了，波前是弯曲的，同一套导向矢量就会出现明显的相位误差，方向图也会变形。做真实声音采集之前，先算这个数，能省掉很多排查问题的时间。

在近场条件下，有两种补救方案：一是改用近场模型，导向矢量里加入距离因子，把阵列流形从“方向”扩展成“方向+距离”；二是增大阵列尺寸时同步拉远测试距离。前者适合麦克风贴脸拾音的场景，后者适合会议室、车载等中远场应用。

2.2 导向矢量的推导和代码里的符号约定

先明确符号约定：假设均匀直线阵沿x轴正方向摆放，第一个阵元位置为参考点，远场平面波从θ方向入射。相邻阵元接收信号的相位差为：

φ = 2πf · d · sinθ / c

于是第m个阵元的相位相对于参考点是 −m·φ（入射波先到参考点，后到第m个阵元），导向矢量写成：

a(θ) = [1, e^{−jφ}, e^{−j2φ}, …, e^{−j(M−1)φ}]^T

在MATLAB里用矩阵运算生成，m从0取到M−1：

function a = steering_vector(theta_deg, f, d, M, c) % theta_deg: 入射角度（度） % f: 频率（Hz） % d: 阵元间距（m） % M: 阵元数 % c: 声速（m/s） idx = (0:M-1)'; phase = 2*pi*f*d/c * idx * sind(theta_deg); a = exp(-1j * phase); % 注意符号与文章约定一致 end

这里最容易出错的是sind和sin的混用。角度统一用度做单位，就老老实实用sind和cosd；如果改成弧度，必须用sin。我自己见过很多次因为弧度换算错误导致方向图整体偏移的代码，排查半天最后发现是“度没转弧度”。建议在代码开头写清楚：“所有角度量使用度数，所有三角函数使用sind/cosd”。

2.3 延迟求和波束形成：最简单也是最重要的方案

延迟求和波束形成（Delay-and-Sum，简称DAS），权向量直接取导向矢量的共轭再归一化：

w_DAS = a(θ_s) / M

它的物理意义非常直接：对θ_s方向来的信号，每个阵元做相位补偿后同相叠加，相当于把阵列的“注意力”集中到该方向。这个方案就是均匀加权w(x) = 1/L的离散实现，所以方向图的旁瓣水平也继承了矩形窗的特性，第一旁瓣约−13.3dB。

波束宽度可以用这个近似公式快速估算：Δθ_3dB ≈ 0.886·λ / (D·cosθ_s)。代入M=8、d=λ/2、θ_s=0°，D=4λ，得到约12.7度。这意味着想分辨两个相隔5度的声源，这个阵列是做不到的。想要更窄的波束，只能增大孔径D——要么加阵元，要么拉大间距，但间距受半波长约束，在固定频率下基本是二选一的博弈。这个公式在做阵列方案选型时非常实用，先算清楚分辨率够不够，再去调算法。

2.4 加窗和自适应加权：同一积分架构下的不同系数

既然方向图就是孔径加权的傅里叶变换，那控制旁瓣就可以直接对w(x)加窗。比如用切比雪夫窗把旁瓣压到−40dB，但代价是主瓣明显变宽。麦克风阵列里有一个经典工程妥协：主瓣宽度和旁瓣高度不可能同时最优，窗函数就是在这个权衡中滑动。实际项目里，如果目标方向附近有强干扰源，选择高旁瓣抑制的窗；如果空间内声源密集、需要角度分辨率，宁可牺牲旁瓣也要保住窄主瓣。

自适应波束形成走得更远：不从“固定方向图”出发，而是从“数据驱动”出发。以MVDR为代表，它的约束条件很简洁——

min w^H R w, s.t. w^H a(θ_s) = 1

其中R是接收信号协方差矩阵，约束保证期望方向响应无失真，目标是最小化输出总功率。闭式解为：

w_MVDR = R⁻¹a(θ_s) / (a^H(θ_s)·R⁻¹·a(θ_s))

这个解就是在同一个线性积分框架下，根据数据统计特性自适应选出的最优权向量。理论上它能在干扰方向自动形成深陷零陷，旁瓣水平远低于13dB。但实际工程里MVDR对阵列流形误差非常敏感，阵元增益失配、相位扰动、方向角估计偏差都会让性能断崖式下跌，远不如固定DAS稳健。所以我的建议是：做算法对比基线时用DAS，做学术演示时用MVDR，做真实硬件系统时先别急着上自适应，把DAS的底子打好再说。

3. MATLAB仿真实现：建模、方向图和信号提取

3.1 仿真参数设计

先把整个仿真平台的参数定下来。我选择M=8阵元、中心频率4kHz、阵元间距恰好等于半波长。这个配置既满足空间采样定理，又不会让波束过窄导致方向图难以观察；采样率取16kHz，满足4kHz实信号采样的奈奎斯特条件。参数总表如下：

干扰信号方向放在−20°，距离主瓣约50度，DAS在这个方向上有一定的空间衰减，但不会特别深——正好能直观展示“固定权向量抑制干扰的能力有限”，为后面讲MVDR留个引子。

3.2 阵列接收信号建模

生成接收信号，核心就一句话：每个阵元收到的信号 = 目标信号 × 该方向导向矢量 + 干扰信号 × 干扰方向导向矢量 + 噪声。注意维度矩阵是M行N列，每一列对应一个采样时刻。

% 基础参数 clear; clc; close all; c = 343; f0 = 4000; lambda = c / f0; M = 8; d = lambda / 2; fs = 16000; T = 0.1; N = round(fs * T); % 1600 个采样点 t = (0:N-1) / fs; theta_s = 30; % 目标方向 theta_j = -20; % 干扰方向 % 目标导向矢量和干扰导向矢量 idx = (0:M-1)'; phase_s = 2*pi*f0*d/c * idx * sind(theta_s); a_s = exp(-1j * phase_s); phase_j = 2*pi*f0*d/c * idx * sind(theta_j); a_j = exp(-1j * phase_j); % 源信号（复指数形式，便于观察相位） s_t = exp(1j*2*pi*f0*t).'; % 目标信号，幅度1 j_t = 0.8 * exp(1j*2*pi*f0*t + 1j*pi/3).'; % 干扰信号，幅度0.8 noise = 0.05 * (randn(M,N) + 1j*randn(M,N)); % 复高斯噪声 % 接收数据矩阵 X，维度 M x N X = a_s * s_t.' + a_j * j_t.' + noise;

这里有个经验之谈：先用单频复指数信号做验证，是最快的跑通流程的方法。你不用关心滤波器设计、包络解调这些细节，只管盯住相位和幅度，方向图、加权、信干比这些核心概念全部一目了然。等单频模型验证通过，再换宽带信号做后续扩展。

3.3 波束方向图计算

方向图就是固定权向量w对扫描方向θ逐一计算的响应幅度。给定权向量w = a_s / M，扫描−90°到90°，每个角度计算该方向的导向矢量，然后取内积的模。

% 延迟求和波束形成权向量 w_das = a_s / M; % 扫描角度方向图 theta_scan = -90:0.5:90; B = zeros(size(theta_scan)); for k = 1:length(theta_scan) phase_k = 2*pi*f0*d/c * idx * sind(theta_scan(k)); a_scan = exp(-1j * phase_k); B(k) = w_das' * a_scan; end % 绘制方向图 figure; plot(theta_scan, 20*log10(abs(B) + eps)); xlabel('角度 (deg)'); ylabel('归一化幅度 (dB)'); grid on; ylim([-40, 5]); title('DAS 波束方向图，期望方向 30°');

跑完这个脚本你会在图上看到：主瓣峰值出现在30°位置，这是方向图正确性的第一个自检点。主瓣宽度从峰值下降约6dB的两个零点之间约十几度，两侧第一旁瓣大约在−13.3dB附近。在整个方向图上，−20°位置对应幅值大概在−10dB到−20dB之间，取决于具体阵列流形，这正是固定波束形成对干扰的“空间衰减”效果。

3.4 目标信号提取和干扰抑制验证

有了权向量，输出信号就是一次矩阵乘法：

y = w_das' * X; % 1 x N 输出信号

为了验证它确实把目标方向信号提取出来了，分别看两个数值：一是输出信号与目标信号s_t的相关系数，二是干扰方向在方向图上的衰减值。

% 相关系数 corr_val = abs(corrcoef(real(y), real(s_t.'))); fprintf('输出与目标信号的相关系数: %.4f\n', corr_val(1,2)); % 干扰方向衰减 phase_j_resp = 2*pi*f0*d/c * idx * sind(theta_j); a_j_resp = exp(-1j * phase_j_resp); resp_j = abs(w_das' * a_j_resp); fprintf('干扰方向空间衰减: %.2f dB\n', 20*log10(resp_j + eps));

实测下来，相关系数会非常接近1，说明期望方向的信号完整保留了；干扰方向衰减通常只有10到20dB。单靠DAS做到这个水平已经不错，毕竟它不需要矩阵求逆、不需要估计背景噪声统计，属于最廉价可靠的方案。但注意：如果干扰落到主瓣内，比如目标30度、干扰40度，固定DAS基本无能为力，这种场景才是自适应波束形成出场的时机。

4. 宽带信号场景：FFT分频积分实现

4.1 为什么单频模型不够用

麦克风阵列处理语音、引擎噪声这类真实信号时，频率是不断变化的。如果还用中心频率f0对应的导向矢量给所有频点加权，会出大问题：波数k = 2πf/c随频率变化，高频信号对应的阵元相位差更大，低频信号相位差更小；同一个权向量扫出去，主瓣在高频偏向一侧，低频又宽得离谱。说白了，固定单频导向矢量处理宽带信号，就像一个用单一焦距的透镜同时去看近景和远景，必然糊。

正确的做法是：把接收信号做FFT，得到每个频点上的窄带信号，然后对每一个频点单独计算该频率下的导向矢量，再做线性积分加权。这正好对应1.2节那个频域积分公式，只是从一次积分变成了每个频点的多次积分。每一种频率的导向矢量都不同，整个波束形成器就演化成一个频率相关的滤波器组。

4.2 逐频点导向矢量加权

用chirp信号做演示，频率从1kHz扫到4kHz。FFT点数选2048，频率分辨率约7.8Hz，足够平滑。代码的核心是一个循环，注意不要偷懒把循环向量化过度，逐频点操作其实更好理解，而且便于后续加频率依赖权重。

% 宽带chirp信号 f_start = 1000; f_end = 4000; s_chirp = chirp(t, f_start, T, f_end, 'linear').'; % 接收数据 X_chirp = a_s * s_chirp.' + a_j * (0.8*chirp(t, f_start, T, f_end, 'linear')).' + noise; % FFT 分频 Nfft = 2048; Xf = fft(X_chirp, Nfft, 2); f_axis = (0:Nfft-1) * fs / Nfft; Yf = zeros(1, Nfft); for k = 2:Nfft/2 fk = f_axis(k); if fk < f_start || fk > f_end continue; end lambda_k = c / fk; % 该频点下的导向矢量 phase_k = 2*pi*fk*d/c * idx * sind(theta_s); a_k = exp(-1j * phase_k); w_k = a_k / M; Yf(k) = w_k' * Xf(:, k); end % 共轭对称重建时域信号 Yf(Nfft/2+2:Nfft) = conj(fliplr(Yf(2:Nfft/2))); y_wide = real(ifft(Yf, 'symmetric')); y_wide = y_wide(1:N);

注意这里我把chirp信号同时用作“目标”和“干扰”，干扰幅度仍然取0.8。ifft(..., 'symmetric')是MATLAB里强制让逆变换结果实数化的小技巧，省得手动抽取实部再担心频谱不对称。

4.3 宽带结果怎么看

先看输出波形：把y_wide和s_chirp叠加画在同一张图上，你会发现输出波形和原chirp几乎重合，但幅度略小——这是加权归一化的正常损耗。再看频谱：用pwelch分别求原信号和输出信号的功率谱，输出里1k~4kHz频段的能量基本被保留，而目标频段外的噪声被压低。

最后评估信干比（SIR）改善情况。单独把目标信号成分和干扰成分分别通过波束形成器：

y_s_only = w_das' * (a_s * s_chirp.'); y_j_only = w_das' * (a_j * (0.8*chirp(t, f_start, T, f_end, 'linear')).'); SIR_in = mean(abs(s_chirp).^2) / mean(abs(0.8*chirp(t, f_start, T, f_end, 'linear')).^2); SIR_out = mean(abs(y_s_only).^2) / mean(abs(y_j_only).^2); fprintf('SIR改善: %.2f dB\n', 10*log10(SIR_out / SIR_in));

宽带场景下，因为每个频点都单独做了积分加权，SIR改善通常比单频模型更稳定。这也是为什么实际麦克风阵列处理语音时，几乎都走“FFT + 分频波束形成”这条路线。

5. 常见问题与排查技巧实录

5.1 空间混叠和栅瓣

空间混叠是阵列信号处理里的头号大坑。当一个方向θ_1的入射波和另一个方向θ_2的入射波在阵元采样后产生完全相同的相位差，方向图上就会在θ_2处出现一个和主瓣等高的虚假响应，这就是栅瓣。判断条件是：

d > λ / (1 + |sinθ_1|)

当d增大或者频率升高（波长变短）时，条件变得更容易触发。一个典型的错误是把阵元间距设成10cm，又用8kHz信号做指向——λ/2是2.1cm，10cm远大于它，方向图上在±60°附近很快冒出栅瓣。处理办法很简单：设计阶段保证最高工作频率下仍然满足d ≤ λ_min/2；如果硬件间距已经固定，就只能限制工作频率上限和扫描角度范围。注意并不是所有栅瓣都完全等高于主瓣，很多时候它以“次栅瓣”的形式只比主瓣低一两个dB，极难发现，所以排查时最好把方向图整个画出来看。

5.2 阵列相位中心不统一

很多人写导向矢量时习惯把第一个阵元作为参考点，代码简单直观。但做真实系统时我强烈建议把坐标系原点放在阵列中心，理由有二：其一，阵列中心参考的导向矢量在数学上具有对称性，方向图响应关于法线对称，便于调试；其二，后续如果要做阵列分割、子阵拼接、甚至是多个阵列间的协同处理，参考点不统一会引入频率无关但角度依赖的相位偏移，极难校准。改进方法是把阵元坐标改成 (−(M−1)d/2, …, (M−1)d/2)：

pos = ((0:M-1) - (M-1)/2)' * d; phase = 2*pi*f/c * pos * sind(theta_deg); a = exp(-1j * phase);

这样做之后，方向图的幅度响应不变，但相位响应变得居中对称。幅度方向图不受影响，所以很多人忽略这一点，但一旦进入自适应处理或多通道级联，麻烦就来了。

5.3 近场条件失效

当你把麦克风阵列放在桌面上，测试用的音箱就在阵列前方0.5米处，瑞利距离算出来却需要2米以上，那平面波假设就已经碎了。这时候方向图上表现出来的现象往往是：实际指向偏离预设角度、主瓣变宽、甚至旁瓣抬升。排查时建议先做一道算术题：根据阵列孔径D和工作频率，用r_far = 2D²/λ算一下临界距离，凡是实测距离小于这个值的场景，都不要强行用远场平面波导向矢量；要么改用近场模型，要么增大测试距离。这里特别提醒：不要因为“我的阵列只有8cm长，0.5米肯定算远场”而想当然，孔径小确实让瑞利距离变小，但高频工作波长更短，两者要一起算。

5.4 阵元增益和相位失配

真实硬件里，每个麦克风通道的增益不可能完全一致，ADC相位也会有偏差。这种失配等效于给权向量w乘上了一个随机扰动，后果是方向图旁瓣被抬高，自适应算法尤其受影响。做一个简单仿真就能看到严重程度：给第二个阵元设置1.2倍增益、叠加15度的相位偏差，再把方向图画出来，原本−13.3dB的均匀旁瓣会被抬高到−6dB甚至更高。工程上应对有两招：一是开机跑一段环境噪声做通道互相关校准，估算增益和相位偏差并补偿；二是在算法层面给权向量加对角加载（diagonal loading），让自适应解对误差不那么敏感。对于DAS这类固定波束形成，幅相失配的影响相对可控，但也别完全忽略——特别是在要求旁瓣低于−20dB的高规格应用里。

最后分享一点实际做下来的体会

这套线性积分视角下的波束形成框架，我用了好几年，最直观的感受是：几乎所有高级算法都能映射回那个积分表达式，区别只是w(x)怎么选。MVDR是数据驱动的统计最优积分核，超指向性是追求极限波束宽度的积分核，加窗则是在主瓣和旁瓣之间做权衡的工程积分核。理解了这一层，再抓任何新算法都很快，因为你会习惯性地先问一句：它到底改了什么加权函数？

做仿真的话，再提醒一个自检小动作：每次跑完方向图，第一件事就是检查主瓣峰值是不是落在预设的θ_s上。这个检查能同时抓住坐标系符号错误、角度弧度混用、权向量共轭问题等一大类低级错误。我自己写波束形成代码时，哪怕跑过一百遍，这个自检也从没省过。这套代码的扩展方向还很多——把第4章的逐频点加权换成GSC自适应结构，就能做语音增强；把方向图扫描扩展到二维平面阵，就能做声源定位成像。骨架就摆在这里，剩下的细节，交给实际应用去打磨。