做相控阵或者阵列信号处理的同学肯定天天跟方向图打交道。算方向图本身不难一个for循环遍历角度每个角度把M个阵元贡献加起来完事。但一旦阵元数量过百、角度扫描分辨率要求高或者要在嵌入式系统里实时刷新方向图这个for循环就会变成性能瓶颈。这也正是“基于FFT的阵列方向图快速计算公式推导”这类工作存在的意义——利用阵列方向图与阵元激励之间天然的离散傅里叶变换关系把O(M×Nθ)的逐点计算降到O(NlogN)在FPGA里更是可以直接用现成的FFT IP核一条流水线出结果。这篇文章基于我实际推导和反复调试代码的经验把公式推导、代码实现和几个特别容易踩的坑一次性讲清楚。适合雷达、通信、天线阵仿真方向的工程师也适合准备在Vivado或者STM32F4上落地FFT方向图算法的同学。1. 阵列方向图计算慢在哪先看常规方法的多重循环1.1 方向图在工程中的高频使用场景方向图是描述阵列在不同角度上辐射响应强弱的曲线。相控阵天线设计阶段要看副瓣电平、波束宽度5G波束赋形要给每个用户实时算一组波束方向图MIMO雷达在做发射方向图匹配、接收方向图置零时更是每迭代一次就要重算一轮方向图。这些场景的共同特点是方向图不是一个算完就结束的结果而是要在优化循环、自适应算法、实时波控表格生成中反复调用。在嵌入式设备里情况就更紧张。比如用STM32F4做小型阵列的波束监测或者用Vivado里的FFT IP核做高帧率方向图上报逐点遍历往往撑不住实时性。我见过不少项目方向图计算代码写得没问题但一跑起来CPU时间全花在这个二重循环上整个系统的刷新率被拖到几百毫秒甚至秒级。这种时候FFT方法就不是“锦上添花”了而是能不能落地的关键。1.2 常规逐点计算的复杂度问题常规做法特别直觉对每个角度θ把所有阵元的贡献叠加起来。F[theta] Σ_{m0}^{M-1} w[m] * exp(j * 2π * (d/λ) * m * sinθ)伪代码写出来就是for theta in angle_list: for m in range(M): F[theta] w[m] * exp(j * 2π * (d/λ) * m * sin(theta))假设阵元数M128角度范围-90°到90°步进0.1°那角度数K1801。一次完整的方向图计算需要执行128×1801≈23万次复数指数运算。看着不多但放在自适应算法里如果要做100次迭代就是2300万次。再如果阵列规模到1024元角度步进缩到0.01°单次就是上亿次运算这在任何实时系统里都是不可接受的。更关键的是这个运算量里大量计算是重复的不同角度对应的复指数项之间其实存在严格的周期结构。如果我们能把角度变量和阵元变量同时离散化并找到合适的变量代换就能把这个二重循环变成一次标准的离散傅里叶变换——FFT的强项就在这。2. 核心公式推导方向图采样就是阵元激励的离散傅里叶变换2.1 均匀线阵方向图公式的离散化写法讨论范围先限定在最常用的均匀线阵ULA。设M个阵元沿x轴等间距排列间距为d阵元位置为x_m m·dm 0, 1, ..., M-1。远场方向图忽略幅度锥削和单元方向图可以写成F(θ) Σ_{m0}^{M-1} w_m · exp(j·(2π/λ)·m·d·sinθ)其中w_m是第m个阵元的复激励λ是工作波长。为了书写方便把u sinθ称为方向余弦它在可视角范围内取值[-1, 1]。代入后F(u) Σ_{m0}^{M-1} w_m · exp(j·2π·(d/λ)·m·u)现在关键问题来了我们要在u域上对方向图进行等间隔采样。假设采样点数为N采样间隔为Δu第k个采样点为u_k k·Δuk0,1,...,N-1。代入上式F(u_k) Σ_{m0}^{M-1} w_m · exp(j·2π·(d/λ)·m·k·Δu)如果令Δu λ/(d·N)那么F(u_k) Σ_{m0}^{M-1} w_m · exp(j·2π·m·k/N)这时你再看离散傅里叶变换的定义X[k] Σ_{n0}^{N-1} x[n] · exp(-j·2π·n·k/N)两个式子的指数符号只差一个负号。如果把方向图公式里的指数符号处理一下等间隔采样方向图就是阵元激励序列w_m的离散傅里叶变换或逆变换。换句话说一阵列的激励做一次FFT/IFFT拿到的就是整个方向图在u域的等间隔采样。这个结论在数学上非常干净但实际工程里有两个细节必须掰扯清楚否则代码写出来方向图会左右颠倒。2.2 变量代换从sinθ到DFT的k轴用DFT视角重新审视上面推导。我们选择的u_k k·λ/(d·N)对应sinθ的等间隔采样。这意味着FFT输出的第k个点对应的不是角度θ本身而是u sinθ。角度需要通过θ arcsin(u)换算回来。这里重点说物理含义阵列方向图在u域上的采样间隔完全由d、λ和N决定。Δu λ/(d·N)。以小角度近似来看θ≈u所以角度域的采样间隔大约是Δθ ≈ λ/(d·N)弧度。dλ/2、N256时Δu 1/(0.5×256) 0.0078125对应角度步进约0.45°。这个精度对绝大多数方向图分析都够了。那u域的覆盖范围呢当k从0到N-1且我们只用FFT的前半段或经过fftshift之后覆盖范围是u ∈ [-λ/(2d), λ/(2d)]。当dλ/2时正好覆盖[-1,1]也就是把整个可见区扫完。如果dλ/2覆盖范围会超过[-1,1]多出来的部分对应“不可见区”实际使用中要截掉。如果dλ/2覆盖范围小于[-1,1]不对我想反了——让我重新说明白dλ/2时λ/(2d)1理论上可见区[-1,1]比一个周期还大方向图在可见区内会出现重复的波束峰也就是栅瓣。这一点在4.2节会详细说。2.3 补零和点数N角度网格密度由什么决定原始阵列只有M个阵元激励序列长度是M。但为了做FFT我们要选择一个变换点数NN通常大于M甚至大得多。这时需要在激励序列w后面补N-M个零再做N点FFT。补零的物理意义是什么补零等效于在u域的连续方向图上做更密集的采样。还是拿上面的公式看u_k k·λ/(d·N)N越大Δu越小方向图曲线上的采样点越密。但注意补零不会提高阵列的真实角分辨率——波束宽度由M、d、λ这些物理参数决定补零只是让方向图的采样网格更细让峰值位置、零陷位置看得更准图形更平滑。这一点特别容易被新手误解成“补零能提高分辨率”其实不能。举例8元均匀线阵dλ/2波束宽度大约为0.886λ/(M·d) ≈ 0.22弧度约12.8°。如果N16Δu1/80.125角度采样间隔约7.2°这个网格比波束宽度还粗方向图峰值可能刚好落在两个采样点之间画出来波形会很丑。补到N128后Δu1/640.0156角度间隔约0.9°曲线就平滑多了。2.4 用FFT还是IFFT符号问题的本质这是最容易踩坑的地方。方向图公式里是正指数exp(j·2π·m·k/N)而FFT的定义是负指数exp(-j·2π·m·k/N)IFFT才是正指数。所以严格来说把激励补零后直接做一次IFFT才能得到与u_k k·λ/(d·N)对应的方向图采样。对应关系写成公式就是F(u_k) N · IFFT(w_padded)[k]其中u_k k·λ/(d·N)如果你非要用FFT也不是不行但输出方向和上面的方向图正好镜像——峰值会出现在u轴的另一侧。修正方法有几种对激励序列取共轭再做FFT得到的方向图再取共轭用FFT但把横坐标取负即把k换成-kCLI命令式地记住一句话“要算正指数加和就找正指数变换IFFT优先”。我个人建议代码里默认用IFFT理由就是方向图公式的指数形式和IFFT完全一致少一层脑内翻转。后面第3节的代码也是这样写的。3. Python实战三步实现FFT方向图并验证3.1 第一步构造激励向量并补零先用一个最简单也最能说明问题的例子16元均匀线阵dλ/2激励采用等幅分布但加上一个指向30°的相移让波束偏转。这样做完FFT后主瓣应该出现在u0.5因为sin30°0.5的位置方便验证。import numpy as np import matplotlib.pyplot as plt M 16 # 阵元数 d 0.5 # 阵元间距以波长为单位 lam 1.0 # 波长归一化为1 theta0 np.deg2rad(30) u0 np.sin(theta0) # 阵元索引 m np.arange(M) # 指向u0的激励指数为负使方向图公式里的指数在uu0处叠加为0 w np.exp(-1j * 2 * np.pi * (d / lam) * m * u0)这里激励指数的符号需要注意。方向图公式F(u) Σw_m·exp(j2π(d/λ)m·u)要取得最大值必须让w_m的相位和exp(j2π(d/λ)m·u0)相互抵消所以w_m exp(-j2π(d/λ)m·u0)。接着补零。N取256把16点激励补齐到256点N 256 # FFT点数建议2的幂 w_pad np.pad(w, (0, N - M)) # 尾部补零到N补零位置放尾部即可因为后续做IFFT时补零在尾部对应u域更高频/更偏角度的位置不影响主瓣中心区域。3.2 第二步做IFFT并生成u轴坐标核心就一行F np.fft.fftshift(np.fft.ifft(w_pad, N)) * N乘以N是为了抵消IFFT自带的1/N因子让幅度回到和方向图公式一致的水平幅度归一化时再除以M。fftshift是把零频率移到数组中心这样方向图的主瓣就落在u0附近符合我们从左到右看曲线的习惯。u轴坐标这样生成u np.fft.fftshift(np.fft.fftfreq(N, dd)) * lamfftfreq产生的频率序列对应DFT的k轴乘以λ/d就是u sinθ。注意fftfreq的第二个参数传入的是阵元间距d归一化波长单位这里和2.2节的代换公式是严格对应的。最终画图只取可见区|u|≤1的部分mask np.abs(u) 1 F_db 20 * np.log10(np.abs(F) / M 1e-12) plt.figure(figsize(8, 4)) plt.plot(u[mask], F_db[mask]) plt.xlabel(u sin(theta)) plt.ylabel(Normalized pattern (dB)) plt.ylim(-40, 5) plt.grid(True) plt.show()运行结果显示主瓣出现在u0.5对应30°副瓣约为-13.3dB和均匀分布阵列的理论副瓣电平完全一致。波束宽度也能从曲线上直接读出来和解析公式对得上。3.3 第三步用逐点法验证结果FFT方法算出来的方向图到底准不准最稳妥的验证方法就是和常规逐点法对比。逐点法公式直接用方向图定义def pattern_direct(w, u): # w: 原始激励u: 方向余弦数组 m np.arange(len(w)) return np.array([ np.sum(w * np.exp(1j * 2 * np.pi * (d / lam) * m * uu)) for uu in u ])然后对比两种方法在u u_grid处的幅度差u_grid np.linspace(-1, 1, 1001) F_direct pattern_direct(w, u_grid) F_fft np.fft.fftshift(np.fft.ifft(w_pad, N)) * N u_fft np.fft.fftshift(np.fft.fftfreq(N, dd)) * lam # 在FFT采样点上插值对比 F_direct_interp np.interp(u_fft[mask], u_grid, np.abs(F_direct)) err np.max(np.abs(np.abs(F_fft[mask]) - F_direct_interp)) print(fMax abs error {err:.2e})实测误差在10的负12次方量级纯粹是浮点舍入误差。这说明只要公式方向和坐标映射正确FFT方法不是近似而是精确的离散化计算和逐点法在对应采样点上完全等价。3.4 参数速查表不同d/λ和N下的u轴间隔做工程时经常要快速估算给定阵元间距和FFT点数角度网格到底能到多细。下面这张表可以直接抄。阵元间距d/λFFT点数Nu域覆盖范围u采样间隔Δu对应角度步进(约)0.564[-1, 1)0.0156250.9°0.5256[-1, 1)0.0039060.22°0.51024[-1, 1)0.0009770.056°0.25256[-2, 2)0.003906可见区外还有冗余0.8256[-0.625, 0.625)0.001953可见区内出现栅瓣从表里能直接看出两个规律。第一N越大角度网格越细但覆盖范围不变半周期λ/(2d)固定。第二d/λ大于0.5后可见区[-1,1]超出了u域的一个工作周期方向图会出现多个主瓣也就是栅瓣。4. 工程落地中的常见问题与避坑经验4.1 补零倍数和FFT点数怎么取才合适补零倍数没有绝对标准但有一个实用的经验法则N取原始阵元数的4到16倍且向上取到2的幂。取4倍以上是为了让方向图峰值附近的采样点足够多画出来的曲线不吃亏取2的幂是为了适配FFT的高效实现以及FPGA里FFT IP核的输入点数要求。如果只是做波束扫描的粗看NM也可以用但峰值可能被低估因为峰值不一定正好落在采样网格上。要做精细的副瓣分析或者零点定位N取16M以上也不夸张反正FFT计算量是NlogN从64点到1024点的计算时间增加并不夸张。实测128点FFT在普通PC上是微秒级在STM32F4上用CMSIS-DSP库做256点CFFT也是百微秒级别性能完全不是问题。有一点必须提醒补零不改变方向图的物理分辨率不要指望把N无限拉大就能分出两个原本叠在一起的波束。阵列的角分辨率由阵列口径M·d决定补零只是“插值”不是“分辨”。4.2 栅瓣出现的位置和DFT周期性用FFT算方向图后有个很有意思的现象如果d/λ0.5方向图输出里会在主瓣之外出现一个形状一模一样的“假峰”这其实就是栅瓣。它的出现原因可以从DFT的周期性角度理解。方向图F(u)作为u的函数周期是λ/d。当dλ/2时周期等于2刚好和可见区宽度[-1,1]一样所以可见区内只有一个完整周期。当d0.8λ时周期变成1.25小于可见区宽度2于是可见区内能看到大约1.6个周期重复出现的波束峰就是栅瓣。FFT输出天然是周期序列所以栅瓣位置根本不需要额外计算它自己就会出现在u轴上对应的地方。这个视角特别有用在做阵列设计时直接看FFT方向图输出就能快速判断间距是否会导致栅瓣而不需要再按传统方式算所谓“栅瓣出现条件”。4.3 从一维到二维均匀面阵方向图的fft2实现一维线阵搞明白了二维均匀面阵就是直接套用二维FFT。设Mx×My矩形栅格面阵x方向间距dxy方向间距dy激励矩阵W[m,n]方向图写作F(u,v) Σ_m Σ_n W[m,n] · exp(j·2π·(dx/λ)·m·u) · exp(j·2π·(dy/λ)·n·v)这本质上是一个二维DFT。代码实现也是在两个维度上分别补零然后做ifft2Mx, My 16, 12 dx dy 0.5 Nx Ny 256 w2d np.ones((Mx, My)) # 实际替换成你的二维激励 w2d_pad np.pad(w2d, ((0, Nx - Mx), (0, Ny - My))) F2 np.fft.fftshift(np.fft.ifft2(w2d_pad, s(Nx, Ny))) * (Nx * Ny) u_ax np.fft.fftshift(np.fft.fftfreq(Nx, ddx)) * lam v_ax np.fft.fftshift(np.fft.fftfreq(Ny, ddy)) * lam注意一个特别容易搞反的点二维激励矩阵的第一个维度对应x方向第二个维度对应y方向而meshgrid生成的坐标网格也是第一维是行y方向、第二维是列x方向。如果你用np.meshgrid生成网格再填激励要确认行列顺序和fft2的维度定义一致否则方向图会转置。我的习惯是统一用“第一维是x、第二维是y”的数组形状然后用extent参数画图避免绕晕。4.4 非均匀阵列、嵌入式移植与FFT IP核的衔接FFT方法最大的局限是要求阵元等间距即阵列是均匀线阵或均匀面阵。非均匀阵列比如稀疏阵、圆形阵、共形阵的阵元位置不满足等差关系方向图公式里不会出现标准的傅里叶变换核直接套FFT必然出错。这时候要么用逐点法要么用非均匀FFTNUFFT或者做相位补偿后分段处理。工程上如果只是稀疏阵有一种近似做法是把空缺阵元位置补零凑成均匀网格再做FFT但这样做方向图会引入额外的栅瓣和误差只适合粗看不适合精细设计。嵌入式落地方面Vivado的FFT IP核配置时通常提供forward和inverse选项结合前面2.4节的结论应该选择inverse模式或者把激励取共轭后走forward模式。IP核输出顺序默认是自然序或位反转序需要在后处理里根据配置手动做重排这一步对应PC端的fftshift但意义略有不同fftshift是把零频搬到中心而位反转重排是把FFT输出的乱序恢复成自然序两者别混。STM32F4虽然没有硬件FFT但CMSIS-DSP库的arm_cfft_f32做实数或复数FFT足够快。它的输入是“实部虚部交替排列”的数组输出也是同样格式而且要求点数必须是2的幂。用这个库做方向图时先把补零后的复激励按实虚交替填进数组调用arm_cfft_f32后再手动把k0分量移到数组中央。这个“手动fftshift”在嵌入式上没有现成函数我自己写了一个简单的索引重排踩过的坑是如果忘了把补零放在激励尾部或者把虚部序列填错位置方向图会直接变成噪声。4.5 归一化、对数和边界处理的小技巧方向图通常用dB表示幅度归一化是除以阵元数M均匀激励时峰值幅度正好等于M。代码里我习惯写成20*log10(abs(F)/M 1e-12)加1e-12是为了防止零点处取对数产生-inf绘图时inf会导致曲线断开。这个细节虽然小但在批量跑方向图时能少很多麻烦。还有一个边界问题u轴上的可见区是|u|≤1但FFT输出的u轴范围是[-λ/(2d), λ/(2d)]当d≠λ/2时u轴范围可能比[-1,1]大也可能小。绘图和计算副瓣电平时一定要先做mask截断否则会把不可见区的电平算进统计里得出一个明显偏低的“假副瓣”。这个坑我在做一版波束扫描程序时踩过当时就是把不可见区的毛刺当成主瓣旁边的副瓣分析了半天最后才发现是mask忘了加。这段内容后续可以这样扩展这套FFT方向图方法看着简单但它能延伸出去的地方不少。比如我最近在折腾自适应波束形成每个迭代步都要重算方向图和梯度用FFT版本之后一次1000点的方向图计算从毫秒级降到了微秒级整个迭代链路瞬间就通顺了。再比如用FFT方向图配合波束空间变换可以快速估计到达角或者把阵列响应矩阵的乘法全部变成FFT流水线在FPGA上做实时处理时资源占用也低得多。我个人在实际操作中的体会是凡是遇到“逐点遍历角度”的方向图计算先别急着优化循环停下来想一想能不能把角度采样和阵元响应的关系写成傅里叶变换的形式——很多时候答案是可以。最开始用FFT算方向图时峰值出现在负角度折腾了一晚上才发现是fft和ifft的符号差异。后来我习惯把“先构造一个已知指向的波束做符号标定”写进代码模板每次换平台都先跑这一步这个习惯帮我避免了很多低级错误。希望这篇推导和排坑记录也能让你少走几步弯路。