阵列信号处理这块DFT波束形成算是“看着简单、上手就翻车”的典型。我最早接触这套东西时以为就是把阵元数据丢进fft函数出来的峰值角度就是来波方向结果第一次跑仿真波束指向偏了七八度还以为是阵列流形推错了。实际上DFT波束形成背后有一套非常明确的适用条件远场、窄带、均匀线阵、阵元间距与波长成固定比例缺一个条件结果就会歪。这篇就围绕这些条件展开把离散傅里叶变换DFT波束形成的数学来路、Matlab实现细节、以及实际调试中哪些地方最容易出错一条一条拆开讲透。内容面向两类人一类是刚学阵列信号处理、需要把公式跑成图的学生和初级工程师另一类是做雷达、通信、超声、声呐方向需要快速搭一个频域波束形成验证链路的从业者。代码全部给到可运行级别参数取值、归一化方式、角度映射关系都会写清楚照着改频率和阵元数就能直接上自己的场景。1. 从延迟求和到空间傅里叶变换DFT波束形成在算什么1.1 远场窄带让时延退化成一个相位因子先把最原始的做法摆出来。一维均匀线阵M 个阵元间距 d平面波从与法线成 θ 的方向入射阵元 m 相对参考阵元的传播时延是 τ_m m·d·sinθ/c。做延迟求和波束形成就是对每个阵元补偿 -τ_m 的时延再相加让同相分量叠起来。问题在于数字系统里精确的时延是做不到的只能做到采样周期的整数倍。补偿非整数采样周期就得用分数延迟滤波器阶数不够就引入幅相误差。而窄带假设给了另一条路如果信号带宽 B 远小于载频 f_c那么一段时延 τ 在复包络上几乎不改变幅度只改变相位即 s(t-τ) ≈ s(t)·e^{-j2πf_c τ}。把这一步代进去时延补偿变成了相位旋转权重就是一个复指数。这一步是整个 DFT 波束形成的根基。它的代价也很明确——窄带近似成立的范围是有边界的第 5 节会专门算这个边界在哪里以及越界之后波束图会怎么坏掉。1.2 导向矢量长得和 DFT 矩阵的列一模一样把上面的相位写成向量形式。取 d λ/2即半波长布阵则 2πf_c·τ_m 2π·(d/λ)·m·sinθ π·m·sinθ。定义空间频率 u sinθ导向矢量就是a(u) [1, e^{jπu}, e^{j2πu}, …, e^{jπ(M-1)u}]^T再看 M 点离散傅里叶变换的基向量w_k [1, e^{j2πk/M}, e^{j2π2k/M}, …, e^{j2π(M-1)k/M}]^T两式一对比只要令 πu 2πk/M也就是 u 2k/M导向矢量和 DFT 基向量就完全重合。这意味着什么意味着对 M 个阵元的同一时刻快拍做一次 M 点 DFT得到的第 k 个输出等价于用指向 u_k 2k/M 的导向矢量做了一次匹配滤波。一次 FFTM 个波束同时出来。这就是 DFT 波束形成最核心的一句话它是空间域的傅里叶变换把阵元域的采样序列变换到空间频率域每个频率 bin 对应一个固定指向的波束。1.3 一次FFT出全部波束的前提条件这个等价关系成立不是无条件的我把它拆成四条缺一条就得改代码条件具体要求不满足会怎样远场平面波假设波前曲率可忽略导向矢量要用球面波模型FFT 的平移不变性失效窄带B/f_c 足够小时延只体现为相位各频点导向矢量不同需逐频点处理等间距阵元位置严格等间隔阵元域不再是均匀采样FFT 无法直接对应d λ/2空间采样率满足奈奎斯特d λ/2 出现栅瓣d λ/2 波束变宽且可见区填不满第四条容易被忽略。空间采样和时域采样是完全对偶的阵元间距 d 对应采样周期1/d 是空间采样率u 域的一个“周期”长度是 λ/d。要让可见区 u ∈ [-1, 1] 恰好覆盖一个完整周期不产生混叠需要 λ/d ≥ 2即 d ≤ λ/2。d λ/2 时可见区刚好铺满 M 个栅格点d 偏大则可见区装不下一个周期u 域出现周期延拓这就是栅瓣d 偏小则一个周期只填了可见区的一部分M 个栅格点挤在 u ∈ [-2d/λ, 2d/λ] 内剩余的可见区没有栅格需要更多阵元才能覆盖。所以那行d lambda/2不是随便写的经验值而是让 FFT 输出恰好无冗余、无混叠地铺满整个可见角的唯一取值。2. Matlab 里把一条线阵的DFT波束形成跑通2.1 参数取法与数据生成我习惯先把物理量定死再写处理链避免中途改参数导致的量纲混乱。下面这套参数规模小、跑得快但包含了真实链路里所有关键环节。c 3e8; % 光速 fc 1e9; % 载频 1 GHz lambda c/fc; % 波长 0.3 m M 16; % 阵元数 d lambda/2; % 半波长布阵 fs 4e9; % 采样率 Ns 2048; % 每阵元采样点数 theta0 20; % 真实来向度 snr_db 20; t (0:Ns-1)/fs; s exp(1j*2*pi*fc*t); % 复指数窄带信号 u0 sind(theta0); % 注意用 sind输入是角度制 tau (0:M-1). * d * u0 / c; % M×1 相对时延 A exp(-1j*2*pi*fc*tau); % M×1 导向矢量 X A * s; % M×Ns 阵列数据 noise (randn(M,Ns) 1j*randn(M,Ns))/sqrt(2); X X 10^(-snr_db/20) * noise; % 按复噪声功率归一后加噪这里有几个坑要提前说清。第一sind和sin千万别混用sind(20)是 20 度sin(20)是 20 弧度后者算出来的 u0 大于 1asin直接报错很多人第一次跑就是这么挂的。第二导向矢量的指数符号是负的因为阵元 m 收到信号比参考阵元晚相位滞后。第三噪声功率按复高斯归一化时要除以 √2否则实部虚部各一份功率总功率翻倍实际信噪比比设定值高 3 dB。2.2 阵元维FFT五行核心代码数据准备好之后空间 DFT 就是一行B fft(X, M, 1); % 沿第 1 维阵元维做 M 点 DFTfft(X, M, 1)里的三个参数分别是数据矩阵、变换点数、变换维度。第三个参数最容易被漏掉。如果 X 是 M×Ns 的矩阵而你不指定维度Matlab 默认沿第一个非单例维度做——恰好也是阵元维看起来没问题。但一旦你把数据转置成 Ns×M 存不写维度就会变成沿时间做 FFT输出看着有值但完全不是波束。变换之后需要把 bin 索引映射回角度k (0:M-1).; u_grid 2*k/M; % 空间频率栅格 u_grid(u_grid 1) u_grid(u_grid 1) - 2; % 折回 [-1,1) theta_grid asind(u_grid); % 空间频率转角度 % 取信号频点上的空间谱 f_axis (0:Ns-1)*fs/Ns; [~, fidx] min(abs(f_axis - fc)); b B(:, fidx); [pk_val, pk_idx] max(abs(b)); fprintf(粗测来向 %.2f 度\n, theta_grid(pk_idx));u_grid 1那一行是关键。M 点 DFT 输出索引 k 从 0 到 M-1对应 u 0, 2/M, ..., 2(M-1)/M。当 u 超过 1 时实际物理意义是负方向的空间频率因为 u 2k/M 与 u 2(k-M)/M 在复指数上完全等价差了整数个 2π。不做这个折回画出来的波束图角度会从 0 单调增到接近 180 度再突跳看着像乱码。还有一点asind的定义域是 [-1,1]浮点误差可能让 u_grid 出现 ±1.0000000002 这种值会返回 NaN。稳妥的做法是加一行u_grid max(min(u_grid,1),-1);做截断。2.3 用显式导向矢量反验FFT版本的正确性我不太相信“看着对”的图一般会用手算的导向矢量去反验一次 FFT 的结果。做法是取 k 0 这一列它对应的指向是 u 0正侧射方向理论上等于把每个阵元的信号直接相加% 方式一FFT 的 k0 输出 y_fft B(1, fidx); % 方式二显式求和 y_sum sum(X(:,fidx)) / M; % 方式三显式导向矢量内积 a0 ones(M,1); y_mvdr (a0 * X(:,fidx)) / M; fprintf(fft%.6f%.6fj sum%.6f%.6fj\n, ... real(y_fft), imag(y_fft), real(y_sum), imag(y_sum));三个值必须完全一致除了 FFT 归一化引入的 M 倍差异。如果对不上问题一定在维度或者索引上先别急着调别的参数。我见过一次对不上的情况查了半天发现是加窗的时候把窗乘到了时间维而不是阵元维——X.*win里 win 是 M×1 列向量Matlab 会自动广播到每一列这是对的但如果 win 写成 1×M 行向量就会广播到每一行变成对时间序列加窗波束图完全不对还很难从图上看出问题。同样地取任意 k 都能验证k_test 3; a_k exp(-1j*pi*(0:M-1).*u_grid(k_test1)); y_direct (a_k * X(:,fidx)) / M; fprintf(k%d fft%.6f%.6fj direct%.6f%.6fj\n, ... k_test, real(B(k_test1,fidx)), imag(B(k_test1,fidx)), ... real(y_direct), imag(y_direct));这个反验习惯帮我省了很多时间尤其是换阵型、加加权、改采样率的时候每次都验一遍能立刻定位问题出在模型还是实现。3. 波束图上的三个“看起来不对”3.1 栅格在 u 域均匀在角度域就不均匀这是个必须建立起来的直觉。M 点 DFT 给出的栅格是 u sinθ 等间隔间隔 2/M。换算成角度θ_k asin(2k/M) 在 0 度附近很密越靠近 ±90 度越稀。M 16 时u 间隔 0.125正侧射附近的角间隔约 7.2 度而接近端射时相邻两点可以差到 20 度以上。for k 0:M/2 fprintf(k%2d u%.3f theta%7.2f 度\n, k, u_grid(k1), theta_grid(k1)); end跑出来会看到 k0 到 k1 差了 7.18 度k6 到 k7 差了 19.5 度。如果拿这套栅格去测一个位于端射附近的信号分辨率会差得离谱——不是算法不行是栅格本身在那儿就粗。处理办法有两个。要么接受 u 域均匀的栅格在 u 域画图很多论文图都是横轴画 sinθ 而不是 θ就是这个原因要么用非均匀的导向矢量扫描代价是丢掉 FFT 的速度优势。我一般两种都出u 域图看主瓣形状和旁瓣电平角度域图看实际可用角度范围内的表现。3.2 真实来向不在栅格上时的扇贝损失导向矢量落在栅格点之间时FFT 输出会有损失这叫扇贝损失scalloping loss。算一下最坏情况真实方向在相邻两栅格正中间u 1/M。此时阵列因子AF(u) (1/M)·sin(πMu/2) / sin(πu/2)代入 u 1/M得 AF (1/M)·sin(π/2)/sin(π/(2M)) ≈ 2/π ≈ 0.637即 -3.92 dB。这个 3.92 dB 是 M 点 DFT 矩形窗条件下的固定值跟 M 大小基本无关M 越大越逼近 2/π。也就是说运气不好的时候测出来的峰值幅度比真实值低将近 4 dB。做幅度估计的场合必须考虑这一点做角度粗测则影响不大。实测验证err_list []; theta_scan linspace(-60, 60, 241); for ii 1:numel(theta_scan) uu sind(theta_scan(ii)); A_tmp exp(-1j*pi*(0:M-1).*uu); Xt A_tmp * s; Bt fft(Xt, M, 1); [pv, pi_] max(abs(Bt(:,fidx))); err_list(ii) 20*log10(pv/M); end plot(theta_scan, err_list); grid on; xlabel(真实来向 (度)); ylabel(峰值输出 (dB));画出来是一条在 0 dB 和 -3.92 dB 之间周期性起伏的曲线起伏的周期正好对应栅格间隔在角度域的分布——中间宽两头密。3.3 补零只是插值不提升分辨率看到栅格太粗第一反应往往是fft(X, 4*M, 1)。这么写没错输出点数变成 4 倍波束图看着平滑多了。但主瓣宽度一点没变。补零的本质是在频域做 sinc 插值只是把已有的频谱曲线画得更密没有增加任何新的信息。我做个直观对比B_M fft(X, M, 1); B_2M fft(X, 2*M, 1); B_4M fft(X, 4*M, 1);在 u 域把三条曲线叠在一起会看到 4M 点的曲线完整穿过了 M 点那些采样点形状完全一样。真正能提升分辨率的只有三条路增大 M加长阵列、提高信噪比、换用超分辨算法MUSIC、ESPRIT 一类。不过补零有一个实际价值降低扇贝损失。补到 2M 点最坏失配从半个栅格降到四分之一个栅格损失从 3.92 dB 降到约 1 dB 以内。这在工程上是很划算的一步只多花一点计算量。4. 旁瓣压不下去的时候窗函数该加在哪一维4.1 矩形窗的 -13.2 dB 是绕不过去的起点均匀加权的 DFT 波束形成第一旁瓣就是经典的 -13.2 dB。在通信、雷达的很多场景里这个电平不够用——强干扰从旁瓣进来会直接把弱目标淹掉。压旁瓣最直接的办法是阵元域加窗。注意“阵元域”这三个字。窗函数要乘在阵元维上也就是数据的行方向这和时域做频谱分析时乘在时间序列上是对偶的。加窗之后win hamming(M); % M×1 列向量 win win / mean(win); % 归一化保证主瓣峰值不因加窗而掉 Bw fft(X .* win, M, 1); % 自动广播到每一列win必须是列向量原因前面说过。归一化用mean(win)而不是sum(win)因为我们最终除以的是阵元数 Msum(win)/M mean(win)这样主瓣峰值保持 0 dB方便对比。4.2 主瓣展宽换旁瓣抑制这是个硬交易几种常用窗的实测数据M 16d λ/2u 域测量窗类型第一旁瓣 (dB)主瓣展宽倍数备注矩形窗-13.21.00分辨率最好旁瓣最高Hamming-42.7约 1.36通用折中滚降快Hann-31.5约 1.50旁瓣滚降 -18 dB/octBlackman-58.1约 1.65旁瓣极低主瓣很宽Chebyshev(30dB)-30.0等电平约 1.25旁瓣等电平可控Chebyshev 窗在阵列处理里特别常用因为它的旁瓣是等电平的——所有旁瓣都压在同一水平不会出现远端旁瓣突然抬高的情况。Matlab 里用chebwin(M, 30)直接生成。win_cheb chebwin(M, 30); win_cheb win_cheb / mean(win_cheb); Bc fft(X .* win_cheb, M, 1);选窗的原则我总结成一句话看干扰落在哪个角度区间。如果干扰来自已知的固定大角度用 Chebyshev 把那个区域的旁瓣压平如果干扰分布未知Hamming 的综合表现最稳如果只关心主瓣内的目标、旁瓣无所谓直接用矩形窗换取最好的分辨率。4.3 加窗之后必须重新标定角度这一点吃过亏。加窗之后主瓣不对称性会变尤其是用非对称窗或者窗长不等于阵元数的时候峰值位置可能发生亚栅格偏移。更常见的问题是加窗后不同指向的波束增益不再一致靠近端射的波束因为有效孔径缩短增益本来就会掉加窗会加剧这个效应。处理办法是做一个角度相关的增益校正表对每个栅格指向先用无噪单频信号跑一遍记录峰值输出存成校正系数实测时除回去。M 16 的时候这个表只有 16 个数开销可以忽略但能让整个视场内的幅度一致性从 ±1.5 dB 收紧到 ±0.2 dB 以内。gain_cal zeros(M,1); for k 1:M a_k exp(-1j*pi*(0:M-1).*u_grid(k)); Xc a_k * s; Bcal fft(Xc .* win, M, 1); gain_cal(k) abs(Bcal(k, fidx)) / M; end Bw_cal Bw ./ gain_cal; % 逐 bin 校正5. 窄带假设能撑多宽从单频DFT走到子带DFT5.1 算一下你的带宽余量还剩多少窄带近似到底能容忍多宽的带宽判据是孔径渡越时间内的信号包络变化要小。阵列孔径 D (M-1)d严格说用 M·d 也可以量级一致最大时延差 τ_max D·sinθ_max/c。带宽 B 对应的最小时变尺度是 1/B要求 τ_max ≪ 1/B。用具体数字感受一下。M 16d 0.15 mD 2.25 mθ 90 度时 τ_max 7.5 ns。若 B 10 MHz1/B 100 ns比值 0.075窄带假设很稳。若 B 400 MHz1/B 2.5 ns比值 3.0窄带假设彻底失效。D (M-1)*d; tau_max D*sind(90)/c; B_list [10e6, 100e6, 400e6, 1e9]; for B B_list fprintf(B%5.0f MHz, tau_max*B %.2f\n, B/1e6, tau_max*B); end经验阈值tau_max*B 0.1时放心用窄带 DFT0.1 ~ 0.3之间会有可见的波束倾斜但方向估计还凑合超过0.3就必须逐频点处理。5.2 时间维FFT加逐频点空间DFT的处理链宽带场景的标准做法是“时间维分帧 FFT 每个频点单独做空间 DFT”。因为不同频率对应的波长不同导向矢量里的 d/λ 项随之变化如果强行用中心频率的导向矢量波束会往低频或高频方向偏。Nfft 1024; Xf fft(X, Nfft, 2); % 沿时间维做 FFT f_ax (0:Nfft-1)*fs/Nfft; f_use f_ax(1:Nfft/2); % 只取正频率一半 P_beam zeros(M, numel(f_use)); for ii 1:numel(f_use) lam_i c / f_use(ii); % 该频点的导向矢量d/lam_i 不再是 0.5 u_ax (0:M-1).*d/lam_i; % 相位系数注意这里不是固定 pi % 直接对频域快拍做匹配用 FFT 等价形式时需按比例缩放 ... end这里有个实现上的细节要说清楚。当 d/λ 不等于 0.5 时空间栅格不再是 u 2k/M而是 u k·λ/(M·d)。频率越低λ 越大同样 M 点 FFT 覆盖的 u 范围越小——视场变窄了。反过来为了在低频覆盖完整视场就得用更多的变换点数。我的处理方式是放弃“通用 FFT”这条捷径改成显式导向矢量矩阵乘theta_vec -60:0.5:60; for ii 1:numel(f_use) lam_i c / f_use(ii); A_mat exp(-1j*2*pi*(0:M-1).*d*sind(theta_vec)/lam_i); % M×Ntheta Xf_i Xf(:, ii); % M×1 P_beam(:,ii) abs(A_mat * Xf_i).^2 / M^2; % 各角度功率 end这是矩阵乘复杂度 O(M·Nθ·Nf)M16、Nθ241、Nf512 的时候大约 200 万次乘加Matlab 里毫秒级。除非阵元数上千或者要做实时一般没必要为了省这点算力去折腾 FFT 的变体。5.3 波束倾斜宽带下最直观的那个坏现象宽带条件下如果用中心频率的导向矢量去做所有频点的加权低频和高频的波束峰值会指向不同角度这叫波束倾斜beam squint。表现是目标明明在 20 度宽带信号合成之后低频分量测出来 23 度高频分量测出来 17 度各频点的测向结果呈一条斜线。验证方法很直接上面那段循环里把每个频点的峰值角度记下来画 角度-频率 二维图ang_est zeros(size(f_use)); for ii 1:numel(f_use) [~, idx] max(P_beam(:,ii)); ang_est(ii) theta_vec(idx); end plot(f_use/1e9, ang_est); grid on; xlabel(频率 (GHz)); ylabel(估计来向 (度));如果是一条水平直线说明处理链正确补偿了频率依赖如果是斜线说明还在用固定导向矢量。这个图我很喜欢用来做快速自检一眼就能看出代码有没有把 d/λ 这个频率项带进去。6. 我踩过的那些坑按出现频率排序6.1 fftshift 用不用取决于你下一步要干什么B fft(X, M, 1)的输出顺序是 k 0, 1, ..., M-1对应 u 从 0 增到接近 2然后折回负值。画图的时候如果不做处理波束图横轴是断的。这时有两种写法第一种用fftshift把零频移到中间B_shift fftshift(B, 1); u_shift (-M/2:M/2-1). * (2/M); % 仅 M 为偶数时成立第二种像我前面那样直接在 u 域折回保持原始顺序画图时用plot(u_grid, ...)。两种都对但绝不能混用——fftshift之后的数组和原始索引对不上如果后面还有B(k1, :)这种按 k 取数的逻辑就会错位。我一般统一用第二种因为它不需要考虑 M 是奇数还是偶数fftshift在 M 为奇数时零频位置是(M1)/2很容易算错。6.2 复信号和实信号的方向会反过来如果生成的是实信号cos(2πf_c t)通过fft之后正负频率各有一份对称谱。做空间 DFT 的时候负频率那一支对应的导向矢量是共轭的指向也反了。表现出来就是明明信号从 20 度来波束图在 20 和 -20 度各出一个峰。解决办法有两个。要么从一开始就用复指数信号建模我推荐这个物理意义清楚不会混淆要么对实信号先做希尔伯特变换取解析信号再走后续流程。X_ana hilbert(X, [], 2); % 沿时间维取解析信号顺便提一句hilbert是沿第一个非单例维度做变换的X 是 M×Ns 时默认沿阵元维必须显式写第二个参数。这个坑我踩过输出看着有值实际上把空间维当成了时间维。6.3 阵元索引起点差一位指向偏半个栅格Matlab 索引从 1 开始但物理推导里阵元编号从 0 开始。写导向矢量的时候如果用(1:M)而不是(0:M-1)相当于整个阵列平移了一个阵元间距波束指向会偏。算一下偏移量相位多出 πu等价于指向平移了 u 域的一个栅格间隔——对 M 16 来说是 0.125在正侧射附近约 7 度。这个错误特别隐蔽因为波束图看着还是正常的只是峰值位置偏了一点点很容易被当成“栅格不够密”而忽略。我的自检办法是用一个已知来向的信号跑一遍看峰值角度和真实角度的偏差是否小于半个栅格间隔。err abs(theta_grid(pk_idx) - theta0); assert(err 180/pi*2/M/2*1.2, 指向偏差超出半栅格检查阵元索引);6.4 归一化和 dB 换算别搞反波束图的纵轴是 dB转换公式20*log10(abs(B)/max(abs(B)))。注意是 20 倍不是 10 倍因为处理的是复幅度而不是功率。如果中间做过abs(...).^2转成功率才用 10 倍。这两种写法在不同代码里混着出现一旦搞混旁瓣电平会显示成 -26 dB 而不是 -13 dB看着像“我的窗函数效果特别好”实际是算错了。还有个归一化的问题fft本身不除以点数所以 M 点 FFT 的输出幅度是时域和的 M 倍。做波束图对比的时候如果一边除了一边没除两条曲线会差 24 dB。我习惯在处理链的最开始就把归一化定死——B fft(X, M, 1)/M之后所有比较都在同一量纲下做。顺带一个实操建议画波束图的时候把纵轴下限定在 -60 dB 左右就够了再往下是数值噪声看不出有效信息反而让旁瓣细节看不清。这套东西真正跑顺之后你会发现 DFT 波束形成的代码量其实很小难点全在参数取值的物理依据和对输出图的正确解读上。我自己现在的做法是把它当成一个“快速原型”先用 FFT 版本几分钟内出个结果确认信号模型、阵列几何、角度映射都没问题再往上叠自适应算法或者宽带处理。这个顺序比一上来就写复杂的加权矩阵要稳得多出问题的时候也更容易定位到底是模型错了还是算法错了。
/YOLO道路破损坑洼检测)
