简介这份资源是一套基于FDTD时域有限差分方法的电磁仿真MATLAB程序聚焦于将仿真得到的近场数据转换为远场数据面向学习电磁计算、天线辐射特性分析的学生与工程人员。压缩包内共1个文件为MATLAB脚本.m整体约4KB体量轻便便于直接阅读与二次修改。程序围绕FDTD近远场转换这一关键环节展开通常涉及网格与时间步长定义、电磁源与边界条件设置、时间迭代更新、近场数据存储以及借助傅里叶变换完成近场到远场的转换与方向图分析可帮助读者理解FDTD算法的基本流程与远场特性提取思路。目前已有828人学习下载适合作为天线设计、雷达与无线通信方向入门电磁仿真的参考脚本也可用于对照调试与算法验证。1. near_to_far_EM.zip 里的近远场变换FDTD 仿真做完远场方向图到底怎么拿到手做 FDTD 仿真的人多半遇到过这个场景微带天线在时域里跑完了S 参数看着还行可一旦要出远场方向图、增益、轴比软件菜单里翻来翻去就是找不到一个能直接用的按钮。近场数据明明就在监视器里躺着远场却像隔了一层毛玻璃。near_to_far_EM.zip这个包名其实已经把答案写在脸上了——它干的就是近远场变换near-to-far-field transformation把 FDTD 计算域内一个闭合面上的近场切向分量通过等效原理和格林函数积分外推到你想要的那个远场观察方向上去。这不是某个软件的专属功能而是一套可以独立实现的电磁后处理流程。适合谁看适合已经能用 FDTD 跑通天线或散射体、但卡在后处理出图这一步的人也适合想自己写脚本、不想被 GUI 绑死的人。下面从原理到代码把这条路走一遍。2. 近远场变换的物理底子等效原理怎么把近场“搬”到远场2.1 为什么近场不能直接当远场用FDTD 是在有限计算域里解 Maxwell 方程你得到的近场包含 evanescent 波和球面波前幅度和相位都随距离剧烈变化。远场定义在 ( kr \to \infty ) 的区域那里只剩横向的辐射场幅度按 ( 1/r ) 衰减方向图形状不再随距离改变。两者之间差了一个积分变换。如果你直接把近场某个点的值当成远场方向图会完全失真——主瓣可能偏移旁瓣电平能差十几 dB。常见做法是选一个包围源的闭合面惠更斯面在这个面上记录切向的 E 和 H然后用表面等效电流 ( \mathbf{J}_s \hat{n} \times \mathbf{H} ) 和等效磁流 ( \mathbf{M}_s -\hat{n} \times \mathbf{E} ) 替代真实源对外部区域做辐射积分。2.2 等效原理与远场积分公式设闭合面 ( S ) 包围所有源外部区域无源。根据等效原理S 外的场可以由 S 上的等效面电流和面磁流产生。远场区 ( \mathbf{r} ) 方向的电场可以写成[ \mathbf{E}_{far}(\hat{r}) -j k \frac{e^{-jkr}}{4\pi r} \int_S \left[ \eta \hat{r} \times (\hat{r} \times \mathbf{J}_s) \hat{r} \times \mathbf{M}_s \right] e^{j k \hat{r} \cdot \mathbf{r}} dS ]其中 ( \eta ) 是波阻抗( k ) 是波数( \mathbf{r} ) 是 S 上的源点位置。这个积分就是近远场变换的核心。实际实现时你不需要在每个远场方向上都做一次面积分——那样太慢。标准做法是先对 S 上的等效流做二维傅里叶变换得到角谱再通过驻定相位法近似到远场。但如果你只是要几个主平面上的方向图直接数值积分也完全可行尤其是用 Python 或 MATLAB 写脚本时。2.3 近场监视器怎么摆位置、尺寸和采样密度惠更斯面的位置很讲究。太靠近源近场变化剧烈采样不够细就会混叠太靠近 PML吸收边界会污染近场数据。我一般把面放在距离辐射体约 ( \lambda/4 ) 到 ( \lambda/2 ) 的位置并且确保它完全包围源不穿过任何金属或介质结构。采样密度方面每个波长至少 15 到 20 个网格点如果做的是宽带仿真按最高频率的波长来算。面的形状可以是长方体或圆柱长方体实现简单圆柱更省内存。在 FDTD 软件里这通常对应一组“近场监视器”或“DFT 监视器”记录面上每个网格点的 E 和 H 的频域复数值。提示如果近场面上有金属或介质穿过等效原理的前提就被破坏了变换结果会完全不可信。摆面之前先确认它悬在自由空间里。3. 从 near_to_far_EM.zip 到可运行脚本近远场变换的最小实现3.1 数据准备从 FDTD 导出近场切向分量假设你已经用 FDTD 跑完了一个微带天线的仿真并且在某个频点上导出了惠更斯面上的近场数据。常见的数据格式是每个网格点上有 Ex、Ey、Ez、Hx、Hy、Hz 六个复数值外加坐标。不同软件导出的排列方式不一样但核心信息就这些。下面是一个用 Python 读取和整理近场数据的示例假设数据存成了简单的文本或 CSV每行是 x, y, z, Ex_real, Ex_imag, ... 这样排列。import numpy as np def load_near_field(filename): 读取近场数据文件返回坐标和复场分量。 假设文件每行格式 x y z Ex_re Ex_im Ey_re Ey_im Ez_re Ez_im Hx_re Hx_im Hy_re Hy_im Hz_re Hz_im data np.loadtxt(filename) coords data[:, 0:3] # 形状 (N, 3) E data[:, 3:9:2] 1j * data[:, 4:9:2] # Ex, Ey, Ez 复数 H data[:, 9:15:2] 1j * data[:, 10:15:2] # Hx, Hy, Hz 复数 return coords, E, H # 使用示例 coords, E_near, H_near load_near_field(huygens_surface.csv) print(f读取到 {len(coords)} 个采样点)这段代码的关键在于把实部和虚部分开存储的复数重新组合。参数说明data[:, 3:9:2]取的是 Ex_re, Ey_re, Ez_re 三列data[:, 4:9:2]取的是对应的虚部。如果你的数据格式不同只需要改列索引。注意坐标单位要和后续波数 ( k ) 的单位一致一般用米。3.2 计算等效面电流和面磁流有了近场 E 和 H下一步是算等效流。你需要知道惠更斯面每个点上的外法向 ( \hat{n} )。对于长方体面法向只有六个方向可以按坐标判断。下面这个函数根据坐标自动判断法向然后计算 ( \mathbf{J}_s ) 和 ( \mathbf{M}_s )。def compute_equivalent_currents(coords, E, H, surface_bounds): 根据坐标判断外法向计算等效面电流 Js 和面磁流 Ms。 surface_bounds: dict包含 xmin, xmax, ymin, ymax, zmin, zmax 返回 Js (N,3) 和 Ms (N,3) 复数数组。 xmin, xmax surface_bounds[xmin], surface_bounds[xmax] ymin, ymax surface_bounds[ymin], surface_bounds[ymax] zmin, zmax surface_bounds[zmin], surface_bounds[zmax] tol 1e-6 # 容差根据网格步长调整 N coords.shape[0] Js np.zeros((N, 3), dtypecomplex) Ms np.zeros((N, 3), dtypecomplex) for i in range(N): x, y, z coords[i] # 判断在哪个面上确定外法向 if abs(x - xmin) tol: n np.array([-1, 0, 0]) elif abs(x - xmax) tol: n np.array([1, 0, 0]) elif abs(y - ymin) tol: n np.array([0, -1, 0]) elif abs(y - ymax) tol: n np.array([0, 1, 0]) elif abs(z - zmin) tol: n np.array([0, 0, -1]) elif abs(z - zmax) tol: n np.array([0, 0, 1]) else: raise ValueError(f点 ({x},{y},{z}) 不在指定的惠更斯面上) Js[i] np.cross(n, H[i]) Ms[i] -np.cross(n, E[i]) return Js, Ms逻辑说明外法向必须指向外部区域也就是远离源的方向。np.cross(n, H)计算面电流-np.cross(n, E)计算面磁流。容差tol要根据你的网格步长设置如果步长是 1 mmtol 取 1e-6 米可能太小因为浮点误差可能更大建议取步长的十分之一左右。这个循环在点数多的时候会慢可以用向量化改写但为了可读性先这样写。3.3 远场积分对每个观察方向做面积分现在到了核心步骤。给定一个远场方向 ( \hat{r} (\sin\theta\cos\phi, \sin\theta\sin\phi, \cos\theta) )计算该方向的远场。公式里的面积分在离散数据上就是求和每个点贡献的相位因子是 ( e^{j k \hat{r} \cdot \mathbf{r}} )。下面是一个直接求和的实现。def compute_far_field(coords, Js, Ms, theta, phi, freq): 计算给定方向 (theta, phi) 的远场电场。 theta, phi 单位弧度。freq 单位Hz。 返回远场电场矢量 (3,) 复数。 c 299792458.0 # 光速 mu0 4 * np.pi * 1e-7 eps0 8.8541878128e-12 eta np.sqrt(mu0 / eps0) # 波阻抗约 376.73 欧姆 k 2 * np.pi * freq / c # 波数 r_hat np.array([ np.sin(theta) * np.cos(phi), np.sin(theta) * np.sin(phi), np.cos(theta) ]) # 预计算相位因子 phase np.exp(1j * k * (coords r_hat)) # 形状 (N,) # 积分项eta * r_hat × (r_hat × Js) r_hat × Ms term_J np.cross(r_hat, np.cross(r_hat, Js)) # (N,3) term_M np.cross(r_hat, Ms) # (N,3) integrand eta * term_J term_M # (N,3) # 对每个分量做加权求和 E_far np.sum(integrand * phase[:, np.newaxis], axis0) # 乘以常数因子 -j k e^{-jkr}/(4 pi r)这里 r 取 1 米归一化 const -1j * k / (4 * np.pi) E_far * const return E_far参数说明theta和phi是球坐标角度freq是仿真频率。coords r_hat是矩阵乘法得到每个点的投影距离。phase是相位延迟。integrand是每个点的贡献。最后乘的常数里省略了 ( e^{-jkr}/r )因为方向图只关心相对幅度和相位距离因子对所有方向相同可以归一化掉。如果你要绝对增益需要把 ( r ) 和输入功率代进去。3.4 扫出方向图主平面和三维方向图有了单方向的远场计算扫一圈就得到方向图。下面代码计算 E 面和 H 面。def compute_pattern(coords, Js, Ms, freq, planeE, n_points181): 计算主平面方向图。 planeE 时 phi0theta 从 0 到 2pi planeH 时 thetapi/2phi 从 0 到 2pi。 返回角度数组和归一化功率方向图dB。 angles np.linspace(0, 2 * np.pi, n_points) pattern np.zeros(n_points) for i, ang in enumerate(angles): if plane E: theta, phi ang, 0.0 else: theta, phi np.pi / 2, ang E_far compute_far_field(coords, Js, Ms, theta, phi, freq) pattern[i] np.sum(np.abs(E_far)**2) # 功率 pattern_dB 10 * np.log10(pattern / np.max(pattern) 1e-30) return np.degrees(angles), pattern_dB逻辑说明功率方向图取电场模平方归一化到最大值再转 dB。加1e-30是防止 log10(0)。n_points181对应 2 度分辨率够看主瓣和第一旁瓣。如果要三维方向图把 theta 和 phi 都扫一遍存成二维数组再用 matplotlib 的plot_surface画出来。4. 避坑与排查近远场变换翻车的五个血泪现场4.1 方向图完全对称或明显错误现象算出来的方向图跟预期完全不符比如全向性太强或者主瓣指向明显偏了。原因最常见的是惠更斯面没有完全包围源或者面穿过了金属结构。另一个可能是法向判断错了比如把内法向当成了外法向导致等效流符号反了。解决先检查面的几何位置确保所有辐射结构都在面内。然后在代码里打印几个点的法向确认指向外部。如果符号反了把 Js 和 Ms 都取负号再试。4.2 高频时方向图出现栅瓣或噪声现象频率升高后方向图出现很多不该有的尖峰或剧烈抖动。原因近场采样密度不够。近远场变换要求每个波长至少 15 到 20 个采样点如果按低频设计的网格在高频下就欠采样了。解决要么加密近场监视器的采样要么在变换前对近场数据做插值。插值可以用 scipy 的RegularGridInterpolator但注意插值会引入误差最好还是从 FDTD 里直接导出更密的网格。4.3 远场幅度随频率剧烈波动现象扫频时远场某个方向的幅度出现不合理的谐振峰或深零点。原因近场数据里包含了 PML 的反射。如果惠更斯面离 PML 太近吸收边界没有完全吸收反射波会污染近场。解决把面往内移至少离 PML 有 10 个网格以上。另外检查 PML 层数是否足够通常 8 层以上且反射系数设到 1e-6 以下。4.4 计算速度慢到无法接受现象扫一个三维方向图要几个小时。原因直接对每个方向做面积分复杂度是 O(N_angles * N_points)。如果 N_points 是几万N_angles 是几千计算量爆炸。解决用角谱法。先对等效流做二维 FFT得到 ( k_x, k_y ) 谱然后通过驻定相位法映射到角度。这样一次 FFT 就能得到所有方向速度提升几个数量级。或者至少用向量化把角度循环也去掉用矩阵乘法一次算完。4.5 与商业软件结果对不上现象自己算的方向图和 HFSS 或 CST 导出的远场数据有差异旁瓣差几个 dB。原因可能是坐标系定义不同或者归一化方式不同。HFSS 导出的远场数据通常是增益dB而你的脚本算的是归一化功率。另外HFSS 可能用了不同的积分面或不同的远场近似。解决先对比主瓣宽度和指向如果这些对上了旁瓣差异可能来自数值精度。检查你的近场数据是否包含了所有极化分量以及是否在同一个频点上。如果差异很大回到 4.1 检查面。注意不同软件对 ( e^{j\omega t} ) 和 ( e^{-j\omega t} ) 的约定不同相位因子的符号会差一个共轭。如果你的结果相位完全反了试试把相位因子取共轭。5. 进阶技巧用角谱法把变换速度提上去以及一个验证习惯直接积分法适合验证和调试但真要扫三维方向图角谱法才是正路。思路是等效流 ( \mathbf{J}_s ) 和 ( \mathbf{M}_s ) 在惠更斯面上是二维分布的对它们做二维傅里叶变换得到谱域 ( \tilde{\mathbf{J}}(k_x, k_y) ) 和 ( \tilde{\mathbf{M}}(k_x, k_y) )。远场方向 ( \hat{r} ) 对应谱域里的 ( k_x k \sin\theta \cos\phi )( k_y k \sin\theta \sin\phi )。然后远场电场可以写成谱域量的代数组合不需要再做面积分。具体实现时注意 FFT 的采样间隔要和空间采样间隔匹配否则会出现混叠。下面是一个简化的角谱法框架。def compute_far_field_spectral(Js_grid, Ms_grid, dx, dy, freq): 角谱法输入 Js 和 Ms 在规则网格上的二维数组形状 (Ny, Nx, 3)。 返回 kx, ky 网格和对应的远场电场谱。 Ny, Nx, _ Js_grid.shape k 2 * np.pi * freq / 299792458.0 # 二维 FFT Js_hat np.fft.fft2(Js_grid, axes(0, 1)) Ms_hat np.fft.fft2(Ms_grid, axes(0, 1)) # 频率轴 kx 2 * np.pi * np.fft.fftfreq(Nx, ddx) ky 2 * np.pi * np.fft.fftfreq(Ny, ddy) KX, KY np.meshgrid(kx, ky) # 只取传播波区域 kx^2 ky^2 k^2 mask (KX**2 KY**2) k**2 kz np.sqrt(np.maximum(k**2 - KX**2 - KY**2, 0)) # 这里省略了从谱到远场的具体代数映射实际实现需要根据极化分解 # 返回谱域数据供后续映射到角度 return KX, KY, kz, Js_hat, Ms_hat, mask这个框架只是起点完整的角谱法需要处理极化分解和驻定相位近似。但核心思想是把面积分换成 FFT复杂度从 O(N_angles * N_points) 降到 O(N_points log N_points)。对于大尺寸问题这是唯一可行的办法。验证方面我养成了一个习惯每次写完近远场变换脚本先拿一个已知解析解的源来测——比如一个短偶极子。短偶极子的远场方向图是 ( \sin\theta )你可以在 FDTD 里放一个点源导出近场跑一遍变换看结果是不是 ( \sin\theta )。如果对不上说明代码里有 bug而不是数据问题。这个习惯帮我省了无数 debugging 时间。另外主瓣宽度和第一旁瓣电平是两个最敏感的指标跟商业软件对比时先看这两个。如果它们对上了基本可以放心。希望帮到你。本文还有配套的精品资源点击获取