DFT波束形成Matlab实现：空域多波束与宽带频域处理

1. 先把概念掰开DFT 波束形成到底指哪一种做法在阵列信号处理这个圈子里DFT 波束形成这个说法其实被两拨人用成了两个不同的东西混着用的时候特别容易出问题。我见过不少人拿着别人的代码跑结果波束图死活对不上最后发现是两个人在讲完全不同的实现。第一拨人说的是空域 DFT 波束形成。均匀线阵ULA的导向矢量里第 n 个阵元的相位是exp(-j2π·n·d·sinθ/λ)这个形式和 DFT 的旋转因子exp(-j2π·k·n/N)长得一模一样。于是把阵元维数据直接送进一个 DFT 矩阵工程上就是fft一次就能同时形成 N 个指向不同方向的波束。这就是所谓的用 FFT 做多波束本质是空域滤波器组或者说波束域变换雷达、声呐、麦克风阵列的快速扫描都用这个套路。第二拨人说的是频域 DFT 波束形成指的是宽带信号处理时先把每个阵元的时域数据分帧做 FFT在每一个频率 bin 上按照该频率对应的波长单独算导向矢量加权再 IFFT 和重叠相加回时域。这里的 DFT 是时间维的作用是让宽带信号在频域分而治之从而规避窄带假设失效的问题。这两件事名字像、代码里都出现fft但一个是空间换角度一个是时间换频率做的东西完全不一样。本文两条线都走一遍因为实际系统里它们经常是串联使用的先做时间维 FFT 把宽带拆成窄带再在窄带内做空间维 DFT 多波束。Matlab 实现这件事上这两条线的坑还各不相同。空域那条主要坑在索引到角度的映射和栅瓣频域那条主要坑在分帧边界和逐频点导向矢量的相位对齐。下面我从参数计算开始一路写到能直接跑起来、能对着波束图验证的实现。2. 阵列模型与参数先算数再写代码2.1 ULA 信号模型与空间频率远场窄带平面波从方向 θ 入射到 N 元均匀线阵阵元间距 d以第 0 个阵元为参考相邻阵元接收信号的时延差是τ d·sinθ / c对应载波相位差2π·d·sinθ/λ。写成导向矢量a(θ) [1, e^{-j2π·d·sinθ/λ}, ..., e^{-j2π·(N-1)·d·sinθ/λ}]^T为了偷懒工程上定义一个无量纲空间频率u d·sinθ/λ导向矢量变成a(u) [1, e^{-j2πu}, ..., e^{-j2π(N-1)u}]^T。这一步变量替换是整个 DFT 波束形成的起点因为它把角度这个非线性的东西变成了一个和阵元下标成线性关系的相位斜率。我记得第一次动手时就是没做这个替换直接拿 θ 去写循环结果波束图在 0° 附近被拉得特别宽端射方向挤成一团。原因很简单sinθ和 θ 不是线性关系等角度步进在 sinθ 域里就是不等间隔采样画出来自然变形。做波束图一定要在 sinθ 或 u 域上画再换算回角度看这两条曲线形状完全不同。2.2 关键参数的手算过程我拿一套具体的数来算后面代码全部围绕这套参数你可以直接照着改。参数符号取值说明阵元数N16决定自由度与孔径载频f03 GHz用于波长换算光速c3e8 m/s载波波长λ00.1 mλ0 c/f0阵元间距d0.05 m取 λ0/2阵列孔径L0.75 mL (N-1)·d等效带宽B300 MHz复基带采样率 400 MHz几个必须手算出来的数半功率波束宽度。均匀加权时约等于0.886·λ/(N·d)弧度。代入0.886 × 0.1 / (16 × 0.05) 0.1108 rad ≈ 6.35°。这个数是阵列的角分辨率量级后面的所有角度分辨能力都不能比它更细想做更细就必须加阵元或者提高频率缩短波长。栅瓣条件。当d λ/2时除了主瓣方向还会有别的方向满足相位差为 2π 的整数倍形成和主瓣一样高的假峰。d λ0/2 是临界值考虑带宽 B 之后最低频率对应的波长最长最高频率对应波长λ_min c/(f0 B/2) 3e8/3.15e9 0.0952 m要求d ≤ λ_min/2 0.0476 m。所以我上面写的 d 0.05 m 在带宽边缘其实略微超标会在大约边缘频率处冒出很小的栅瓣。这个细节很少有人在示例里提但真实宽带系统里必须按λ_min/2来定阵元间距否则带内高频端方向图会翘起来。孔径渡越时间。τ_max (N-1)·d/c 0.75/3e8 2.5 ns。判断窄带近似是否成立的经验判据是B·τ_max 1。这里300e6 × 2.5e-9 0.75已经不算是远小于 1了属于必须做频域处理的区间。如果按窄带那样全带宽只用一个 λ0带宽边缘会出问题具体误差下面第 4 节会量化。2.3 窄带假设的体检什么时候不能直接空域 DFT窄带假设的核心是信号在整个孔径上的相位变化可以用单一频率的相位差来近似。数学上就是包络在阵元之间不发生形变。信号带宽越宽、孔径越长包络在孔径上的变化越剧烈这个近似就越站不住。我把这条判据做成一个简单的心算规则先算 B·τ_max小于 0.1 就放心用窄带0.1 到 1 之间要留神做波束图时会看到主瓣变形大于 1 就必须走频域分 bin 的路子。举个具体的失败案例。假设只用一个导向矢量按 f0 算去处理 300 MHz 带宽的信号最高频率处每个阵元的相位偏差是2π·d·sinθ·Δf/c取 sinθ 1、Δf 150 MHz得到每个阵元偏2π × 0.05 × 1.5e8 / 3e8 0.157 rad。16 个阵元一路累积下来最多差15 × 0.157 2.36 rad超过 π意味着阵列两端阵元的相位已经完全反相。这时候波束不仅没对准增益还可能掉到负值。所以宽带系统里用 f0 算一组权重打天下这种做法基本必挂。3. 空域 DFT 多波束的 Matlab 实现3.1 导向矢量与 DFT 矩阵的那层关系把a(u)和 N 点 DFT 矩阵的第 k 行放一起看DFT 矩阵第 k 行第 n 列元素是e^{-j2πkn/N}导向矢量第 n 个元素是e^{-j2πnu}。两者对上当u k/N时完全相等。也就是说第 k 个 DFT 输出对应的空间频率是 u k/N对应角度满足d·sinθ/λ k/N即sinθ_k kλ/(Nd)。取 d λ/2 代进去sinθ_k 2k/N。k 从 0 到 N-1sinθ 从 0 一直走到接近 2。因为 sinθ 的物理范围只有 [-1, 1]超过 1 的部分自动折返到负角度上这就是 DFT 的周期性在空间域的体现。用fftshift把零频挪到中间后索引 k -N/2 到 N/2-1 就整齐对应sinθ 2k/N从 -1 走到 1覆盖全空域。d λ/2 这个取值特别舒服就是因为 N 点 DFT 的角度网格刚好铺满整个可见空间不多不少。换成 d λ/3DFT 网格只覆盖 ±2/3 的可见空间剩下的角度你根本扫不到得靠补零或者额外算内积。3.2 一段可直接跑的波束扫描代码这段代码我改过好几版最后留下的是最不容易出错的那种写法导向矢量统一用列向量波束图矩阵按波束 × 角度组织画图前先归一化到 dB。% ---- 参数 ---- N 16; % 阵元数 d 0.05; % 阵元间距 (m) c 3e8; f0 3e9; lam0 c/f0; % ---- 波束索引 - 角度 (d lam0/2 时) ---- k (-N/2 : N/2-1).; sin_beam 2*k/N; % 空间频率网格 theta_beam asind(sin_beam); % 对应角度 (deg) % ---- 空间 DFT 矩阵 (酉矩阵, 带 1/sqrt(N) 归一化) ---- n (0:N-1).; Fdft exp(-1j*2*pi*n*k./N) / sqrt(N); % N x N Fdft fftshift(Fdft, 2); % 行顺序与 sin_beam 对齐 % ---- 用一条单频平面波验证 ---- theta0 20; % 真实来向 x exp(-1j*2*pi*d*n*sind(theta0)/lam0); % N x 1 快拍 y Fdft * x; % N x 1 多波束输出 [~, idx] max(abs(y)); fprintf(估计角度 %.2f deg\n, theta_beam(idx));跑出来theta_beam(idx)应该是 20°允许有网格量化误差因为波束网格步进是2/N 0.125的 sinθ对应角度间隔在 20° 附近大约是 5.7°。这个量化误差是 DFT 波束形成的固有代价不是代码写错了。想更细就得在阵元维补零也就是在fft里把点数从 N 提到 4N 甚至 8N相当于对波束图做插值角度网格变密但分辨率不变后文会说明这个变密不等于变好的区别。3.3 波束图、栅瓣与加窗的取舍验证实现对不对最靠谱的是画波束图而不是只看峰值位置。波束图的计算是B(θ) |w_k^H · a(θ)|扫一遍 θ 网格。sin_grid linspace(-1, 1, 4001); a_grid exp(-1j*2*pi*d*n*sin_grid/lam0); % N x L % 取第 k0 个波束的权重Fdft 的一行 k0 8; % 约 0 deg w Fdft(k0, :).; B abs(w * a_grid); % 1 x L B 20*log10(B / max(B)); plot(sin_grid, B); xlabel(sin\theta); ylabel(dB); grid on;三个必须看的地方主瓣宽度是否和手算的 6.35° 对得上。如果明显更宽先检查你是不是把w的共轭弄反了取错共轭会让整个方向图镜像翻转看起来能对准但角度是反的。旁瓣电平。均匀加权也就是纯 DFT第一旁瓣约 -13.2 dB这个数很经典。想要更低就得加窗比如改用 Hamming 或者 Taylor 加权。代价是主瓣展宽Hamming 大概让主瓣宽 1.3 倍旁瓣压到 -43 dB 左右。这是典型的用分辨率换动态范围探测强目标旁边的小目标时压旁瓣比抠那点分辨率更值钱。栅瓣。把 d 改成 0.08 m 再跑一遍你会看到 ±90° 附近冒出几个和主瓣差不多的尖峰。这就是空域混叠没救只能改 d。注意加窗的时候窗函数是加在阵元维上的也就是w w .* win不是加在时间上。很多人会下意识套用一维信号的窗对着时间序列做那就完全跑偏了。另外一个实用经验是波束之间其实并不干净地互不干扰。虽然 N 个 DFT 波束在sinθ 2k/N这些离散点上满足正交性F 是酉矩阵能量守恒波束域变换不损失信息但相邻波束在主瓣区域的重叠非常严重一个目标通常会同时点亮 2 到 3 个相邻波束。做角度估计时用峰值位置能凑合想要精细角度就得上插值或者别的超分辨方法别指望单靠 DFT 网格能抠出亚度级精度。4. 宽带频域 DFT 波束形成STFT 域逐频点加权4.1 分帧、FFT、逐频点加权、重叠相加宽带处理的骨架就是四步我用 300 MHz 带宽、400 MHz 复基带采样率来写每个阵元的复基带数据分帧帧长 Nfft帧间重叠 50%。每帧做 FFT得到X_n[m]n 是阵元m 是频率 bin。对每个 bin 单独算这个频率下的导向矢量加权求和。结果 IFFT 回时域重叠相加WOLA拼回连续输出。第二步和第三步是关键。第 m 个 bin 对应的物理频率是f_m f0 f_bin[m]该频率下的波长λ_m c/f_m导向矢量必须用λ_m而不是λ0。fs 400e6; Nfft 1024; hop Nfft/2; n (0:N-1).; % 频率 bin 转成 [-fs/2, fs/2) fb (0:Nfft-1)./Nfft*fs; fb(fb fs/2) fb(fb fs/2) - fs; f_phys f0 fb; % 实际射频频率 lam c ./ f_phys; % 每个 bin 的波长 % 期望方向的阵元时延 (N x 1) tau0 n * d * sind(theta0) / c; % 逐频点导向矢量: A 是 N x Nfft A exp(-1j*2*pi*(f_phys.) .* tau0); % 这里 (1xNfft).*(Nx1) 自动扩展 % 波束形成权重 (匹配滤波形式, 1/N 归一化) W conj(A) / N; % N x NfftW里每一列对应一个频率 bin 的 N 个阵元权重结构非常清晰同一条代码覆盖了所有子带不需要 for 循环逐 bin 算。这样写不只快更重要的是不容易出维度错。帧处理部分win sqrt(hann(Nfft, periodic)); % 分析合成各用一次, 平方和补偿 frame_idx 1:hop:(size(x,1)-Nfft1); y_out zeros(size(x,1)Nfft, 1); win_sq_sum zeros(size(y_out)); for ii 1:numel(frame_idx) s frame_idx(ii); Xf fft(x(s:sNfft-1, :) .* win, Nfft, 1); % Nfft x N Yf (conj(W). * Xf.).; % Nfft x 1 逐 bin 加权 yt real(ifft(Yf)); % 复基带也可保留虚部 seg (s:sNfft-1).; y_out(seg) y_out(seg) yt .* win; win_sq_sum(seg) win_sq_sum(seg) win.^2; end % 归一化, 消除窗的重叠增益 valid win_sq_sum 1e-8; y_out(valid) y_out(valid) ./ win_sq_sum(valid);注意最后那步归一化。如果不做窗平方和补偿帧重叠区域的幅度会比单帧区域高合成信号会周期性起伏听起来像突突的调幅看起来像一段段方波。这是最常见的翻车点之一我第一次做 WOLA 就被这个折腾了半天一直以为是波束形成的相位算错了。4.2 频域导向矢量的符号约定与对齐这一节讲清楚两个容易乱的符号问题。第一相位符号。复基带接收信号模型是x_n(t) s(t-τ_n)·e^{-j2πf0·τ_n}。傅里叶变换后得到X_n(f) S(f)·e^{-j2π(f0f)τ_n}。想让 N 个阵元相干相加权重必须取conj(a_n) e^{j2π(f0f)τ_n}。也就是权重是导向矢量的共轭因为波束形成输出是w^H x共轭转置里已经含了一个共轭。这一点在复数域代码里非常容易写混尤其当你在不同函数之间传参的时候。我建议统一成权重矩阵 W 直接对应w^H x别在中间到处转置。第二帧对齐。所有阵元必须用同一时间窗切帧。如果为了省事让某个通道独立做 STFT 并取不同起点就会引入额外的帧间相位差。避免方法是把多通道数据组织成[采样点 × 通道]的矩阵一次fft出来物理上就不可能错位。第三块卷积的群延迟。如果你不走 STFT 加权、而是把每个频点的权重 IFFT 成长度 Nfft 的 FIR 滤波器、再去做分块卷积那么每路滤波器会引入Nfft/2的群延迟必须统一补偿否则不同通道的时延就被你人为改掉了波束直接指向错误方向。这是 overlap-save 实现里最阴的一个坑因为你从波形上完全看不出异常只是增益低了几 dB很容易被归咎于噪声。4.3 完整实现与验证方法验证频域实现对不对最直接的办法是拿一个已知方向、已知波形的宽带信号走一遍完整流程看两件事波形是否保真、角度是否对准。% 构造宽带线性调频信号 (复基带) t (0:1/fs:200e-6).; s exp(1j*pi*(B/ (t(end)) )*t.^2); % LFM, 带宽 B % 各阵元时延 tau0 n * d * sind(theta0) / c; x zeros(numel(t), N); for nn 1:N x(:, nn) interp1(t, s, t - tau0(nn), linear, 0) ... .* exp(-1j*2*pi*f0*tau0(nn)); end注意interp1那里的t - tau0(nn)时延是亚采样级的纳秒级的时延在 400 MHz 采样下不足一个采样点直接整数移位会损失精度。用线性或者样条插值补误差可以接受要求更高就上分数延迟滤波器。跑完后对比三组结果处理方式主瓣宽度峰值增益带宽边缘表现全带单权重窄带假设偏宽且畸变明显下降严重失配频域逐 bin 加权接近理论值接近 10log10(N)正常频域 阵元加窗宽 1.3 倍略降正常旁瓣低理论增益10log10(16) 12 dB。实测如果只有 8 dB 左右先查是不是导向矢量用错了频率再查帧对齐。这个 12 dB 是个很好的体检指标白噪声背景下做一次比对一下输入输出信噪比差得多说明实现有问题。还有一点频域处理天然支持每个 bin 单独设计权重这意味着你可以在某些 bin 上做干扰抑制例如把已知的窄带干扰所在 bin 的权重置零或者做投影这在时域 DFT 波束形成里是几乎免费的副产品。做声呐或者雷达抗干扰的时候这个能力比波束本身更有价值。5. 从 DAS 到 MVDRDFT 框架下还能怎么往上加5.1 协方差矩阵估计与对角加载DFT 波束形成本质是数据独立的DAS延时相加权重只看几何不看数据。好处是稳健、不需要训练样本、算得快坏处是旁瓣高、抗干扰能力差一个强干扰落在旁瓣上就能把你的目标淹掉。在频域框架下升级到 MVDRCapon几乎是顺手的每个 bin 上估一个协方差矩阵然后解一个约束优化。% 在某个频点 m 上累积 K 个快拍 Xm Xf(m, :).; % N x 1, 单快拍示例 % 实际用 K 帧累积: Xm_all 是 N x K R (Xm_all * Xm_all) / K; % 对角加载, 提高小样本下的稳健性 dl trace(R)/N * 1e-2; R R dl*eye(N); a_m A(:, m); % 该频点的导向矢量 w_mvdr (R \ a_m) / (a_m * (R \ a_m));几个实打实的经验快拍数至少要是 2N 到 3N否则 R 是病态的MVDR 会自消把目标也一起压掉。N 16 的话至少 32 帧稳妥起见用 64 帧以上。对角加载量不要省。我一般取trace(R)/N的 1% 到 10%具体看信噪比。加载太少小样本下方向图会出现尖锐的凹陷甚至主瓣畸变加载太多退化成 DAS白忙一场。这是一个纯靠试出来的参数建议在 -10 dB 到 20 dB 信噪比区间各跑一遍看主瓣是否稳定。白噪声增益要注意。DAS 的 WNG 就是 N12 dBMVDR 在抑制干扰时 WNG 会掉掉太多意味着它对误差极度敏感。实际系统里通道幅相误差一般有百分之几WNG 掉到 0 dB 以下基本就不能用了。5.2 两种方法的实测对比我在仿真里放了一个目标20°和一个比目标强 20 dB 的干扰-30°分别用 DAS 和 MVDR 处理结果差异非常直观指标DASDFT 波束MVDR对角加载 1%目标方向增益12.0 dB12.0 dB约束保证干扰方向响应-13 dB落在旁瓣-40 dB 以下零陷输出信干噪比差提升约 25 dB需要的快拍数160 左右对通道误差敏感度低高表格里干扰方向响应 -13 dB这一栏很说明问题如果你的目标信噪比比干扰低 20 dB那干扰经过旁瓣进来之后仍然和目标一样强DAS 根本没戏。这也是为什么实际系统里 DAS 通常只作为第一级粗处理后面必须接自适应。5.3 通道幅相误差的现场校正说到误差这个必须单拎出来讲。理论仿真里通道是完全一致的实测里每个射频通道的增益差可能有 ±0.5 dB、相位差 ±5° 甚至更大。相位误差是最致命的因为波束形成全靠相位相干5° 的随机相位误差在 16 元阵上累积起来相干增益能掉 1 到 2 dB旁瓣抬高 3 到 5 dB。校正的常规做法是注入一个已知方向的参考信号近场耦合或者远场辅助源测出每个通道相对参考通道的复增益% 假设 ref 是参考通道数据, x 是所有通道 G mean(x .* conj(ref), 1) ./ mean(abs(ref).^2); x_cal x ./ G; % 逐通道复除, 同时校正幅度和相位这一步看着简单但有两个注意点必须在工作频带内做如果只在单频点校宽带系统的高频端还是错的得逐 bin 校校正系数要定期刷新温度变化会让相位漂移尤其是长时间工作的系统。我在一套设备上就吃过这个亏冷机校正完跑了两小时之后方向图明显变差一开始以为是算法问题查了半天才发现是功放温漂导致的通道相位变化。6. 常见问题速查与踩坑记录6.1 症状 - 原因 - 处置对照表症状最可能的原因处置波束指向和真实来向反号导向矢量/权重共轭搞反了检查w^H x里的共轭位置用单频信号验证波束图出现多个等高尖峰d λ_min/2空域混叠按最高频率重定阵元间距或做稀疏阵设计主瓣比理论值宽很多孔径不足或加窗过重核对0.886λ/(Nd)减小窗的滚降输出幅度周期性起伏WOLA 窗平方和未归一化加win_sq_sum补偿宽带边缘角度偏移用了 f0 的导向矢量改为逐 bin 用λ(f)计算相干增益只有 6 dB通道相位未校准或帧错位注入参考信号校幅相检查分帧索引MVDR 把目标也压掉了快拍不足导致 R 病态增加快拍到 3N 以上加大对角加载端射方向波束特别胖孔径在端射方向投影变短正常物理现象端射附近避开使用补零后角度网格变密但分辨不开补零只是插值不是增加信息需要真分辨率就加阵元或提高频率6.2 我个人踩过的几个坑第一个坑是角度映射。我一开始直接theta linspace(-90, 90, N)和 DFT 输出一一对应跑出来的波束图完全对不上教科书。折腾一晚上才反应过来DFT 网格是等 sinθ 间隔的不是等角度间隔。修正之后还有个新问题d 不取 λ/2 的时候网格覆盖不了全空域我又得手动去补零。现在我的习惯是写死一行注释在代码顶上% 角度网格: sin(theta_k) k*lam/(N*d), 索引需 fftshift再也不犯。第二个坑是复基带和实信号搞混。用实信号做 FFT 的时候正负频率是镜像的只处理正频率半边是对的但权重也得按正频率算不能两边都用同一个。我有一版代码忘了这个导致波束出现对称的双峰一开始还以为是两个目标。第三个坑是对补零能提高分辨率的误解。早年我以为fft(x, 4*N)就能让波束变窄实际上它只是把波束图插值画得更细两个靠得很近的目标该糊还是糊。真正的角分辨率由孔径决定也就是那个0.886λ/(Nd)任何后处理都突破不了除非上超分辨算法那是另一套东西。这个认识转过弯之后我在阵列设计阶段就把孔径算够省掉很多后期折腾。第四个坑是端射方向。用的 d λ/2端射方向±90°的波束会明显变胖、增益下降因为此时阵列的有效投影孔径接近零。如果你的应用需要覆盖大角度范围要么改用平面阵要么老老实实接受端射方向性能衰减别在那儿调参数调不出来的。最后分享一个实用的验证套路。不管是空域 DFT 还是频域 DFT 实现我都会先做一个三点测试放三个不同方向的单频源看波束图峰值位置、看主瓣宽度、看旁瓣电平三个都对上了才算过关。这套流程闭着眼睛执行绝大多数实现错误十分钟内就能定位。