阵列信号处理这块，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 的导向矢量做了一次匹配滤波。一次 FFT，M 个波束同时出来。

这就是 DFT 波束形成最核心的一句话：它是空间域的傅里叶变换，把阵元域的采样序列变换到空间频率域，每个频率 bin 对应一个固定指向的波束。

1.3 一次FFT出全部波束的前提条件

这个等价关系成立不是无条件的，我把它拆成四条，缺一条就得改代码：

第四条容易被忽略。空间采样和时域采样是完全对偶的：阵元间距 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 大于 1，asin直接报错，很多人第一次跑就是这么挂的。第二，导向矢量的指数符号是负的，因为阵元 m 收到信号比参考阵元晚，相位滞后。第三，噪声功率按复高斯归一化时要除以 √2，否则实部虚部各一份功率，总功率翻倍，实际信噪比比设定值高 3 dB。

2.2 阵元维FFT：五行核心代码

数据准备好之后，空间 DFT 就是一行：

B = fft(X, M, 1); % 沿第 1 维（阵元维）做 M 点 DFT

fft(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 的 k=0 输出 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_test+1)); y_direct = (a_k' * X(:,fidx)) / M; fprintf('k=%d fft=%.6f%+.6fj direct=%.6f%+.6fj\n', ... k_test, real(B(k_test+1,fidx)), imag(B(k_test+1,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(k+1), theta_grid(k+1)); end

跑出来会看到 k=0 到 k=1 差了 7.18 度，k=6 到 k=7 差了 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)，因为我们最终除以的是阵元数 M，sum(win)/M = mean(win)，这样主瓣峰值保持 0 dB，方便对比。

4.2 主瓣展宽换旁瓣抑制，这是个硬交易

几种常用窗的实测数据（M = 16，d = λ/2，u 域测量）：

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走到子带DFT

5.1 算一下你的带宽余量还剩多少

窄带近似到底能容忍多宽的带宽？判据是孔径渡越时间内的信号包络变化要小。阵列孔径 D = (M-1)d（严格说用 M·d 也可以，量级一致），最大时延差 τ_max = D·sinθ_max/c。带宽 B 对应的最小时变尺度是 1/B，要求 τ_max ≪ 1/B。

用具体数字感受一下。M = 16，d = 0.15 m，D = 2.25 m，θ = 90 度时 τ_max = 7.5 ns。若 B = 10 MHz，1/B = 100 ns，比值 0.075，窄带假设很稳。若 B = 400 MHz，1/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时放心用窄带 DFT；0.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)，M=16、Nθ=241、Nf=512 的时候大约 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(k+1, :)这种按 k 取数的逻辑，就会错位。

我一般统一用第二种，因为它不需要考虑 M 是奇数还是偶数（fftshift在 M 为奇数时零频位置是(M+1)/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 版本几分钟内出个结果，确认信号模型、阵列几何、角度映射都没问题，再往上叠自适应算法或者宽带处理。这个顺序比一上来就写复杂的加权矩阵要稳得多，出问题的时候也更容易定位到底是模型错了还是算法错了。