简介一份面向CST超材料仿真与参数反演的MATLAB源码包适用于电磁仿真、射频微波领域的学习者与研究者核心目标是解决从CST仿真模型中提取S参数并反演超材料几何参数的问题。压缩包内仅含1个m脚本文件整包大小仅1KB脚本中涵盖S参数读取、数据预处理及反演算法调用等环节可配合自建的CST模型和导出的S参数数据直接运行学习。项目知识点覆盖S参数的物理意义与计算方法、超材料电磁响应与负折射率等奇异特性、CST三维建模仿真流程以及通过优化算法反演单元尺寸、填充率等结构参数。已有299人学习下载这份源码可作为了解MATLAB与CST联合仿真的入门参考帮助快速理解参数反演的逆问题求解思路并在实际项目中迁移使用。1. 从 CST 到 MATLAB超材料 S 参数反演的最小闭环CST 里把超材料单胞模型跑完S21 曲线光滑漂亮可一旦想从 S 参数反推等效介电常数和磁导率立刻就会出现折射率虚部为正、阻抗实部为负这类“看起来就不物理”的结果。问题往往不在仿真设置而是许多教程直接套反余弦公式后没有处理多值分支也没有把 CST 导出的复 S 参数做相位连续化处理。下面我拆一个常见源码包里的get_S_Parameter.m它正是用来干这件事的读取 CST 导出的 S 参数用 NRW 方法反演出超材料的等效 n、ε 和 μ。适合做超表面、FSS、吸波体、天线罩材料仿真的工程师和研究生也适合想把 CST 和 MATLAB 联合仿真流程整合的人。2. 均匀性假设与 NRW 反演公式为什么 S21 的相位必须展开2.1 适用前提单层均匀等效介质超材料反演的第一步不是打开 MATLAB而是确认模型能不能被当作“均匀平板”。CST 中通常建立周期边界条件下的单胞模型入射波方向上的厚度 d 通常取结构层高度或者包含基板的整体厚度。只有在这个厚度远小于波长的前提下把单胞等效成均匀介质层才有意义。实际中当单元尺寸达到工作波长的三分之一以上时反演结果会明显出现频散异常这时别说 ε、μ连 S 参数本身都需要重新检查边界条件。还有一个常被忽略的前提单胞在入射方向上应该是左右对称的。对称结构满足 S11 S22S21 S12NRW 的简化公式才成立。如果你的结构是单面金属图案加介质基板正反入射的反射并不完全相同源码里如果只读取 S11 和 S21 就会把非对称信息丢掉。我一般会在 CST 里设置两个波端口导出完整的四组 S 参数后再做对称性判别只有偏差小于 1e-3 时才走简化流程。另外本源文件名的MáS我按 MATLAB 脚本的简写理解它只负责后处理不参与 CST 建模所以流程上要先在 CST 里把 S 参数导出来再交给这个脚本。2.2 NRW 方法的两个关键表达式假设等效介质厚度为 d自由空间波数 k0 2πf/c0。Nicolson-Ross-Weir 方法的核心是两个式子z sqrt( ((1 S11)^2 - S21^2) / ((1 - S11)^2 - S21^2) ) n acos( (1 - S11^2 S21^2) / (2 * S21) ) / (k0 * d)z 是相对自由空间的归一化阻抗n 是等效折射率。得到 z 和 n 后再用ε n/z、μ n * z还原出等效介电常数和磁导率。很多人直接拿这两个公式去算结果要么 n 的虚部为负要么 ε 的实部在某个频点突然跳一个大台阶。原因很简单反余弦函数的值域限制在 [0, π]真实传输相位超过 π 之后会折叠回这个区间直接除以 k0d 得到的不是真实的 n而是 n 加上一个整数倍的2π/(k0d)。这个多值性必须处理。常见做法是先对 S21 的相位做 unwrap然后选择一个整数 m让 n 在频带上连续。下面这段是源码里最核心的分支修正逻辑也是整个反演过程最容易出错的地方% 未修正的折射率n0 已经包含折叠 n0 acos((1 - S11.^2 S21.^2) ./ (2 * S21)) ./ k0d; % 用相位导出的粗略折射率作为基准 phase_S21 unwrap(angle(S21)); n_phase -phase_S21 ./ k0d; % 对每个频点把 n0 修正到与 n_phase 差最小 m round((real(n_phase) - real(n0)) ./ (2 * pi ./ k0d)); n n0 m .* (2 * pi ./ k0d);这里unwrap给出随频率连续变化的 S21 相位n_phase是基于连续相位的折射率参考值每个频点上真实 n 与n_phase之间只可能相差整数倍的2π/(k0d)。用round求出整数 m然后把折叠的n0拉回正确分支。需要说明的是round的结果在高损耗频点可能不稳定因为此时 S21 模值很小相位噪声会被放大。稳妥的做法是把n_phase先做平滑或者只在 |S21| 大于 0.05 的频段应用分支修正。2.3 反演参数的物理约束待求量物理约束反演中如何处理z无源材料要求 Re(z) 0取两个候选符号中实部为正的那个n无源材料要求 Im(n) 0若虚部为负取共轭m整数用连续性约束和相位基线求 roundε, μ极高频段要趋近于 1远离工作频段时候选分支更难选建议只反演局域频带另外需要注意当目标的等效折射率为负也就是左手材料常见情况时上面的 Im(n) 0 判据仍然成立但 Re(z) 0、Re(ε) 0、Re(μ) 0 同时存在会让 z 的正号判据失效。此时如果建模目标确实是负折射率我会在分支选择后额外检查 ε 和 μ 实部是否同号若同号且都为负再把 n 的实部取反。大多数超材料反演源码不会处理这一步所以会出现“仿真能用、反演得到负介电常数但折射率却为正”的怪异结果。3. get_S_Parameter.m 首个难关把 CST 导出的 S 参数读进 MATLAB3.1 CST 导出前的端口设置与参考面在 CST 里得到 S 参数的路径很多但如果你最终要交给 MATLAB 做参数反演必须先把参考面定准。波端口默认安装在端口面上而端口面距离超材料表面通常有一段空气层用于避免高次模这段空气层会引入额外相位延迟。正确做法是在波端口属性里勾选 Deembed把参考面平移到介质表面或者在导出后在 MATLAB 里用传输线理论做 phase de-embedding。get_S_Parameter.m源码里一般不包含 de-embedding 计算它默认你已经处理好参考面这一步最容易忽略。另一个设置是频率范围。CST 默认扫描频点可能不够密导致 S21 相位在相邻频点之间跳变超过一个周期MATLAB 的unwrap很可能无法正确还原。我一般把频率步长设置为工作频带的百分之一以下例如 2-18 GHz 范围取 0.02 GHz 步长。零点不需要参与反演因为 0 Hz 处 S 参数对超材料等效参数没有意义反而可能让后面acos产生除零问题。3.2 导出 Touchstone 格式CST 支持导出 Touchstone.s2p文件也支持导出 ASCII 文本。建议导出.s2pMATLAB 处理这个格式相对成熟。导出时最关键的是数据格式选择S-Parameter Format 有 Magnitude/Phase、Real/Imaginary 等选项。我通常选 Real/Imaginary也就是 RI 格式原因有两点一是能直接拿到复平面上的实部和虚部构建复数一步到位二是幅度和相位格式在相位跨越 ±180° 时容易产生人为不连续虽然数据本身没问题但后续unwrap处理多一道风险。一个完整的.s2p文件头部通常长这样! Created by CST Studio Suite # GHz S RI R 50 10.000000 -0.345620 0.012345 0.876512 -0.123456 ... 10.020000 -0.345512 0.012333 0.876498 -0.123455 ...第一行#之后依次是频率单位、S 参数格式、电阻/阻抗类型、参考阻抗。R 50表示参考电阻为 50 Ω。读取时需要跳过!和#开头的注释行然后按每行 9 个实数的格式解析。CST 有时会导出包含[Two-Port Data]的分段格式但标准单频扫描的.s2p通常是这种行式布局。CST 导出选项推荐值原因文件格式Touchstone.s2pMATLAB 解析简单信息完整数据格式Real/Imaginary避免相位折叠引入额外处理参考阻抗50 Ω和 NRW 公式默认归一化一致频率步长工作频带 / 100 以下保证unwrap没有跨周期跳变3.3 一个不依赖 RF 工具箱的读取函数实验室里并不是每台机器都有 MATLAB RF Toolbox所以get_S_Parameter.m这类源码更常见的做法是自己写解析器。下面这段代码可以直接放到脚本里读取 CST 导出的.s2p或普通文本function data read_s2p_ri(filename) fid fopen(filename, r); if fid -1 error(无法打开文件: %s, filename); end vals []; while ~feof(fid) line strtrim(fgetl(fid)); if isempty(line) || line(1) ! || line(1) # continue; end if line(1) [ continue; % 分段格式这里先跳过 end nums sscanf(line, %f); if numel(nums) 9 vals(end1, :) nums(1:9); %#ok elseif numel(nums) 5 vals(end1, :) [nums(1:5), NaN, NaN, NaN, NaN]; %#ok end end fclose(fid); if isempty(vals) error(文件中没有有效S参数数据); end data.f vals(:,1); data.S11 vals(:,2) 1i*vals(:,3); data.S21 vals(:,4) 1i*vals(:,5); if size(vals,2) 9 data.S12 vals(:,6) 1i*vals(:,7); data.S22 vals(:,8) 1i*vals(:,9); end end这段代码做了三件重要的事第一用strtrim去掉行首尾空格避免sscanf把空行读成空数组第二遇到[开头的分段格式行直接跳过防止解析器崩溃第三对于只有 S11 和 S21 的非完整数据用 NaN 补齐 S12、S22保证主程序里不会出现维度不匹配。之后在主脚本里直接调用d 0.0016; % 结构厚度单位米 sp read_s2p_ri(unit_cell.s2p); f sp.f; S11 sp.S11; S21 sp.S21; [eps_r, mu_r, n, z] invert_nrw(S11, S21, f, d);这里d是唯一需要手动确认的量它必须和 CST 模型在电磁波入射方向上的实际厚度一致。很多初学者把晶格周期当作d导致算出来的 ε 整段向上平移。对于有介质基板的超表面d就是基板厚度加金属贴片厚度如果金属可视为零厚度则直接取基板厚度。建议在 CST 里量一下模型包围盒的 Z 方向尺寸不要凭设计参数猜测。4. 参数反演的两步法从复 S 参数到 ε、μ、n4.1 完整反演函数读取到 S11、S21 后核心反演逻辑可以封装成一个独立函数。下面这段代码就是我整理的 NRW 实现主体也是get_S_Parameter.m在反演阶段最常见的写法function [eps_r, mu_r, n, z] invert_nrw(S11, S21, f, d) c0 299792458; k0 2 * pi * f / c0; k0d k0 * d; % 第一步阻抗 z z1 sqrt(((1 S11).^2 - S21.^2) ./ ((1 - S11).^2 - S21.^2)); z2 -z1; z z1; z(real(z) 0) z2(real(z) 0); % 第二步反余弦折射率 cosnk (1 - S11.^2 S21.^2) ./ (2 * S21); cosnk min(max(cosnk, -1), 1); % 数值截断防止 acos 出现 NaN n0 acos(cosnk) ./ k0d; % 第三步用 S21 相位展开作为分支参考 phase_S21 unwrap(angle(S21)); n_phase -phase_S21 ./ k0d; m round((real(n_phase) - real(n0)) ./ (2 * pi ./ k0d)); n n0 m .* (2 * pi ./ k0d); % 第四步无源约束 n(imag(n) 0) conj(n(imag(n) 0)); eps_r n ./ z; mu_r n .* z; end这个函数把上一章的分支修正直接放到了代码里。代码里z(real(z) 0) z2(real(z) 0)用的是逻辑索引比 for 循环简洁但要求 S11、S21 都是列向量。如果读取函数返回的是行向量运行时会报错所以建议在整个流程开头统一f f(:); S11 S11(:); S21 S21(:);。cosnk的截断范围是 [-1, 1]因为反余弦定义域就在这个区间CST 算出的 S 参数在高损耗频点可能让表达式模值略大于 1不截断的话acos会返回复数后面的unwrap和round全乱套。4.2 分支修正与相位展开的关系unwrap(angle(S21))解决的是“角度折叠”n0中的反余弦还带来“空间周期折叠”两者不能混为一谈。角度折叠是单位圆上的 2π 跳变MATLAB 的unwrap只能在相邻采样间隔不超过 π 时可靠工作空间周期折叠则是传播路径在长度 d 上多累积了 2π 相位它随频率变化。所以即便 S21 相位被unwrap完全展开n0仍可能整体偏移若干个2π/(k0d)。现象来源处理方法S21 相位在 ±π 间跳变反正切函数的周期性MATLABunwrapn 的整体偏移反余弦多值性用n_phase作为基线round求整数倍n 在个别频点突变|S21| 太小噪声放大平滑n_phase或限制频段ε、μ 的实部为负且 n 为正负折射率材料分支错误检查 ε、μ 实部同号必要时将 n 实部取反表格给出四类典型反演异常。反演完成后我建议先画 n 的实部和虚部曲线。理想情况下实部应当平滑且随频率变化符合物理直觉虚部应当非负。如果虚部出现锯齿多半是unwrap在某个频点失去相位连续性处理办法是减小 CST 频率步长后重新导出.s2p而不是在 MATLAB 里做暴力平滑那只会掩盖问题。4.3 负折射率材料的分支处理上面代码里的z(real(z) 0) z2(real(z) 0)假设了无源材料 Re(z) 0这在普通介质和大多数超材料中都成立。但对于负折射率设计Re(n) 0 和 Re(ε) 0、Re(μ) 0 同时存在z 本身仍可保持正值。此时分支修正的 m 计算可能出错n_phase来自传输相位它的实部天然为负而n0来自反余弦主值在除以 k0d 后通常是正数两者相差的不再是整数倍周期。这种情况我一般先按常规流程反演一次得到 n 后检查real(eps).*real(mu)的符号。如果乘积为负说明 n 的符号和解不一致直接把 m 整体减 1重新计算 n、ε、μ。更严谨的做法是回到麦克斯韦方程组的因果性约束但工程上这个“先反演、再判号、再平移一个周期”的流程足够处理大多数左手材料单胞也更容易在脚本里自动实现。5. 最后一公里批量扫描、结果验证与常见返工5.1 用 CST 宏和 MATLAB 循环做批量扫描超材料设计几乎都要扫参数方环宽度、开口间隙、介质厚度。每次手动在 CST 里改参数、导出.s2p、再跑 MATLAB会让人崩溃。CST 的 VBA 宏可以把整个流程录制下来我通常会录制一次仿真和导出操作然后把参数改成一个循环变量。下面是一个典型骨架不同 CST 版本方法名略有差异最好先宏录制一次再替换循环体Sub BatchScan() Dim app As Object Dim i As Integer Set app CSTStudio.Application For i 1 To 8 app.StoreParameter gap, 0.1 * i app.Rebuild app.Solver.Start app.ExportSParameterFile D:\sim\gap_ i .s2p Next i End Sub在 CST 中运行宏之前先确保模型参数名是gap并且已经在参数列表里定义。ExportSParameterFile的路径必须提前存在否则会静默失败。宏跑完后MATLAB 侧用dir批量读取并反演files dir(D:\sim\*.s2p); d 0.0016; for k 1:numel(files) sp read_s2p_ri(fullfile(files(k).folder, files(k).name)); [eps_r, mu_r, n, z] invert_nrw(sp.S11, sp.S21, sp.f, d); [~, tag] fileparts(files(k).name); save(sprintf(%s_invert.mat, tag), eps_r, mu_r, n, z, sp); end这里每轮循环会把原始 S 参数和反演结果一起存进.mat文件后续画损耗角正切、色散曲线都不用重新读.s2p。做参数扫描时我习惯把gap值写进文件名这样从文件名就能重建扫描条件和频点省掉一份额外记录表。5.2 反演结果自检拿 ε、μ 回填 CST最有效的验证不是看曲线是否平滑而是把反演得到的 ε_r、μ_r 作为均匀介质块的参数填回 CST再算一次 S 参数然后与原始超材料单胞的 S 参数对比。具体做法是在 CST 中新建一个尺寸相同的矩形块材料类型设为 Dispersive Material把某个频点反演出的复 ε、复 μ 填进去仿真同一频点如果 S21 幅度差超过 0.1、相位差超过 5°说明这个超材料单元在目标频带内不能等效为均匀介质需要缩小单元尺寸或重新检查亚波长条件。这是我在所有反演脚本里都会保留的一步。很多源码包不提供验证脚本导致你反演出的 ε 看起来再漂亮也可能只是分支选择“恰好连续”而已。如果回填后 S21 对不上优先检查的不是 MATLAB 代码而是 CST 里的 reference plane 有没有校准以及超材料单元是否真的满足周期边界下电尺寸小于 λ/4 的要求。最后提醒一个非常实际的操作批量扫描时不要让 MATLAB 和 CST 同时打开同一批.s2p文件。CST 导出文件时可能还没写完缓冲区MATLAB 读到的数据会缺行我吃过这个亏后会在 CST 宏里对每个导出文件做一次文件大小检查或者在 MATLAB 侧用pause(2)延迟读取确认文件大小不再变化后再进反演流程。本文还有配套的精品资源点击获取 超材料S参数反演详解：CST导出到MATLAB的NRW相位展开与分支修正
简介一份面向CST超材料仿真与参数反演的MATLAB源码包适用于电磁仿真、射频微波领域的学习者与研究者核心目标是解决从CST仿真模型中提取S参数并反演超材料几何参数的问题。压缩包内仅含1个m脚本文件整包大小仅1KB脚本中涵盖S参数读取、数据预处理及反演算法调用等环节可配合自建的CST模型和导出的S参数数据直接运行学习。项目知识点覆盖S参数的物理意义与计算方法、超材料电磁响应与负折射率等奇异特性、CST三维建模仿真流程以及通过优化算法反演单元尺寸、填充率等结构参数。已有299人学习下载这份源码可作为了解MATLAB与CST联合仿真的入门参考帮助快速理解参数反演的逆问题求解思路并在实际项目中迁移使用。1. 从 CST 到 MATLAB超材料 S 参数反演的最小闭环CST 里把超材料单胞模型跑完S21 曲线光滑漂亮可一旦想从 S 参数反推等效介电常数和磁导率立刻就会出现折射率虚部为正、阻抗实部为负这类“看起来就不物理”的结果。问题往往不在仿真设置而是许多教程直接套反余弦公式后没有处理多值分支也没有把 CST 导出的复 S 参数做相位连续化处理。下面我拆一个常见源码包里的get_S_Parameter.m它正是用来干这件事的读取 CST 导出的 S 参数用 NRW 方法反演出超材料的等效 n、ε 和 μ。适合做超表面、FSS、吸波体、天线罩材料仿真的工程师和研究生也适合想把 CST 和 MATLAB 联合仿真流程整合的人。2. 均匀性假设与 NRW 反演公式为什么 S21 的相位必须展开2.1 适用前提单层均匀等效介质超材料反演的第一步不是打开 MATLAB而是确认模型能不能被当作“均匀平板”。CST 中通常建立周期边界条件下的单胞模型入射波方向上的厚度 d 通常取结构层高度或者包含基板的整体厚度。只有在这个厚度远小于波长的前提下把单胞等效成均匀介质层才有意义。实际中当单元尺寸达到工作波长的三分之一以上时反演结果会明显出现频散异常这时别说 ε、μ连 S 参数本身都需要重新检查边界条件。还有一个常被忽略的前提单胞在入射方向上应该是左右对称的。对称结构满足 S11 S22S21 S12NRW 的简化公式才成立。如果你的结构是单面金属图案加介质基板正反入射的反射并不完全相同源码里如果只读取 S11 和 S21 就会把非对称信息丢掉。我一般会在 CST 里设置两个波端口导出完整的四组 S 参数后再做对称性判别只有偏差小于 1e-3 时才走简化流程。另外本源文件名的MáS我按 MATLAB 脚本的简写理解它只负责后处理不参与 CST 建模所以流程上要先在 CST 里把 S 参数导出来再交给这个脚本。2.2 NRW 方法的两个关键表达式假设等效介质厚度为 d自由空间波数 k0 2πf/c0。Nicolson-Ross-Weir 方法的核心是两个式子z sqrt( ((1 S11)^2 - S21^2) / ((1 - S11)^2 - S21^2) ) n acos( (1 - S11^2 S21^2) / (2 * S21) ) / (k0 * d)z 是相对自由空间的归一化阻抗n 是等效折射率。得到 z 和 n 后再用ε n/z、μ n * z还原出等效介电常数和磁导率。很多人直接拿这两个公式去算结果要么 n 的虚部为负要么 ε 的实部在某个频点突然跳一个大台阶。原因很简单反余弦函数的值域限制在 [0, π]真实传输相位超过 π 之后会折叠回这个区间直接除以 k0d 得到的不是真实的 n而是 n 加上一个整数倍的2π/(k0d)。这个多值性必须处理。常见做法是先对 S21 的相位做 unwrap然后选择一个整数 m让 n 在频带上连续。下面这段是源码里最核心的分支修正逻辑也是整个反演过程最容易出错的地方% 未修正的折射率n0 已经包含折叠 n0 acos((1 - S11.^2 S21.^2) ./ (2 * S21)) ./ k0d; % 用相位导出的粗略折射率作为基准 phase_S21 unwrap(angle(S21)); n_phase -phase_S21 ./ k0d; % 对每个频点把 n0 修正到与 n_phase 差最小 m round((real(n_phase) - real(n0)) ./ (2 * pi ./ k0d)); n n0 m .* (2 * pi ./ k0d);这里unwrap给出随频率连续变化的 S21 相位n_phase是基于连续相位的折射率参考值每个频点上真实 n 与n_phase之间只可能相差整数倍的2π/(k0d)。用round求出整数 m然后把折叠的n0拉回正确分支。需要说明的是round的结果在高损耗频点可能不稳定因为此时 S21 模值很小相位噪声会被放大。稳妥的做法是把n_phase先做平滑或者只在 |S21| 大于 0.05 的频段应用分支修正。2.3 反演参数的物理约束待求量物理约束反演中如何处理z无源材料要求 Re(z) 0取两个候选符号中实部为正的那个n无源材料要求 Im(n) 0若虚部为负取共轭m整数用连续性约束和相位基线求 roundε, μ极高频段要趋近于 1远离工作频段时候选分支更难选建议只反演局域频带另外需要注意当目标的等效折射率为负也就是左手材料常见情况时上面的 Im(n) 0 判据仍然成立但 Re(z) 0、Re(ε) 0、Re(μ) 0 同时存在会让 z 的正号判据失效。此时如果建模目标确实是负折射率我会在分支选择后额外检查 ε 和 μ 实部是否同号若同号且都为负再把 n 的实部取反。大多数超材料反演源码不会处理这一步所以会出现“仿真能用、反演得到负介电常数但折射率却为正”的怪异结果。3. get_S_Parameter.m 首个难关把 CST 导出的 S 参数读进 MATLAB3.1 CST 导出前的端口设置与参考面在 CST 里得到 S 参数的路径很多但如果你最终要交给 MATLAB 做参数反演必须先把参考面定准。波端口默认安装在端口面上而端口面距离超材料表面通常有一段空气层用于避免高次模这段空气层会引入额外相位延迟。正确做法是在波端口属性里勾选 Deembed把参考面平移到介质表面或者在导出后在 MATLAB 里用传输线理论做 phase de-embedding。get_S_Parameter.m源码里一般不包含 de-embedding 计算它默认你已经处理好参考面这一步最容易忽略。另一个设置是频率范围。CST 默认扫描频点可能不够密导致 S21 相位在相邻频点之间跳变超过一个周期MATLAB 的unwrap很可能无法正确还原。我一般把频率步长设置为工作频带的百分之一以下例如 2-18 GHz 范围取 0.02 GHz 步长。零点不需要参与反演因为 0 Hz 处 S 参数对超材料等效参数没有意义反而可能让后面acos产生除零问题。3.2 导出 Touchstone 格式CST 支持导出 Touchstone.s2p文件也支持导出 ASCII 文本。建议导出.s2pMATLAB 处理这个格式相对成熟。导出时最关键的是数据格式选择S-Parameter Format 有 Magnitude/Phase、Real/Imaginary 等选项。我通常选 Real/Imaginary也就是 RI 格式原因有两点一是能直接拿到复平面上的实部和虚部构建复数一步到位二是幅度和相位格式在相位跨越 ±180° 时容易产生人为不连续虽然数据本身没问题但后续unwrap处理多一道风险。一个完整的.s2p文件头部通常长这样! Created by CST Studio Suite # GHz S RI R 50 10.000000 -0.345620 0.012345 0.876512 -0.123456 ... 10.020000 -0.345512 0.012333 0.876498 -0.123455 ...第一行#之后依次是频率单位、S 参数格式、电阻/阻抗类型、参考阻抗。R 50表示参考电阻为 50 Ω。读取时需要跳过!和#开头的注释行然后按每行 9 个实数的格式解析。CST 有时会导出包含[Two-Port Data]的分段格式但标准单频扫描的.s2p通常是这种行式布局。CST 导出选项推荐值原因文件格式Touchstone.s2pMATLAB 解析简单信息完整数据格式Real/Imaginary避免相位折叠引入额外处理参考阻抗50 Ω和 NRW 公式默认归一化一致频率步长工作频带 / 100 以下保证unwrap没有跨周期跳变3.3 一个不依赖 RF 工具箱的读取函数实验室里并不是每台机器都有 MATLAB RF Toolbox所以get_S_Parameter.m这类源码更常见的做法是自己写解析器。下面这段代码可以直接放到脚本里读取 CST 导出的.s2p或普通文本function data read_s2p_ri(filename) fid fopen(filename, r); if fid -1 error(无法打开文件: %s, filename); end vals []; while ~feof(fid) line strtrim(fgetl(fid)); if isempty(line) || line(1) ! || line(1) # continue; end if line(1) [ continue; % 分段格式这里先跳过 end nums sscanf(line, %f); if numel(nums) 9 vals(end1, :) nums(1:9); %#ok elseif numel(nums) 5 vals(end1, :) [nums(1:5), NaN, NaN, NaN, NaN]; %#ok end end fclose(fid); if isempty(vals) error(文件中没有有效S参数数据); end data.f vals(:,1); data.S11 vals(:,2) 1i*vals(:,3); data.S21 vals(:,4) 1i*vals(:,5); if size(vals,2) 9 data.S12 vals(:,6) 1i*vals(:,7); data.S22 vals(:,8) 1i*vals(:,9); end end这段代码做了三件重要的事第一用strtrim去掉行首尾空格避免sscanf把空行读成空数组第二遇到[开头的分段格式行直接跳过防止解析器崩溃第三对于只有 S11 和 S21 的非完整数据用 NaN 补齐 S12、S22保证主程序里不会出现维度不匹配。之后在主脚本里直接调用d 0.0016; % 结构厚度单位米 sp read_s2p_ri(unit_cell.s2p); f sp.f; S11 sp.S11; S21 sp.S21; [eps_r, mu_r, n, z] invert_nrw(S11, S21, f, d);这里d是唯一需要手动确认的量它必须和 CST 模型在电磁波入射方向上的实际厚度一致。很多初学者把晶格周期当作d导致算出来的 ε 整段向上平移。对于有介质基板的超表面d就是基板厚度加金属贴片厚度如果金属可视为零厚度则直接取基板厚度。建议在 CST 里量一下模型包围盒的 Z 方向尺寸不要凭设计参数猜测。4. 参数反演的两步法从复 S 参数到 ε、μ、n4.1 完整反演函数读取到 S11、S21 后核心反演逻辑可以封装成一个独立函数。下面这段代码就是我整理的 NRW 实现主体也是get_S_Parameter.m在反演阶段最常见的写法function [eps_r, mu_r, n, z] invert_nrw(S11, S21, f, d) c0 299792458; k0 2 * pi * f / c0; k0d k0 * d; % 第一步阻抗 z z1 sqrt(((1 S11).^2 - S21.^2) ./ ((1 - S11).^2 - S21.^2)); z2 -z1; z z1; z(real(z) 0) z2(real(z) 0); % 第二步反余弦折射率 cosnk (1 - S11.^2 S21.^2) ./ (2 * S21); cosnk min(max(cosnk, -1), 1); % 数值截断防止 acos 出现 NaN n0 acos(cosnk) ./ k0d; % 第三步用 S21 相位展开作为分支参考 phase_S21 unwrap(angle(S21)); n_phase -phase_S21 ./ k0d; m round((real(n_phase) - real(n0)) ./ (2 * pi ./ k0d)); n n0 m .* (2 * pi ./ k0d); % 第四步无源约束 n(imag(n) 0) conj(n(imag(n) 0)); eps_r n ./ z; mu_r n .* z; end这个函数把上一章的分支修正直接放到了代码里。代码里z(real(z) 0) z2(real(z) 0)用的是逻辑索引比 for 循环简洁但要求 S11、S21 都是列向量。如果读取函数返回的是行向量运行时会报错所以建议在整个流程开头统一f f(:); S11 S11(:); S21 S21(:);。cosnk的截断范围是 [-1, 1]因为反余弦定义域就在这个区间CST 算出的 S 参数在高损耗频点可能让表达式模值略大于 1不截断的话acos会返回复数后面的unwrap和round全乱套。4.2 分支修正与相位展开的关系unwrap(angle(S21))解决的是“角度折叠”n0中的反余弦还带来“空间周期折叠”两者不能混为一谈。角度折叠是单位圆上的 2π 跳变MATLAB 的unwrap只能在相邻采样间隔不超过 π 时可靠工作空间周期折叠则是传播路径在长度 d 上多累积了 2π 相位它随频率变化。所以即便 S21 相位被unwrap完全展开n0仍可能整体偏移若干个2π/(k0d)。现象来源处理方法S21 相位在 ±π 间跳变反正切函数的周期性MATLABunwrapn 的整体偏移反余弦多值性用n_phase作为基线round求整数倍n 在个别频点突变|S21| 太小噪声放大平滑n_phase或限制频段ε、μ 的实部为负且 n 为正负折射率材料分支错误检查 ε、μ 实部同号必要时将 n 实部取反表格给出四类典型反演异常。反演完成后我建议先画 n 的实部和虚部曲线。理想情况下实部应当平滑且随频率变化符合物理直觉虚部应当非负。如果虚部出现锯齿多半是unwrap在某个频点失去相位连续性处理办法是减小 CST 频率步长后重新导出.s2p而不是在 MATLAB 里做暴力平滑那只会掩盖问题。4.3 负折射率材料的分支处理上面代码里的z(real(z) 0) z2(real(z) 0)假设了无源材料 Re(z) 0这在普通介质和大多数超材料中都成立。但对于负折射率设计Re(n) 0 和 Re(ε) 0、Re(μ) 0 同时存在z 本身仍可保持正值。此时分支修正的 m 计算可能出错n_phase来自传输相位它的实部天然为负而n0来自反余弦主值在除以 k0d 后通常是正数两者相差的不再是整数倍周期。这种情况我一般先按常规流程反演一次得到 n 后检查real(eps).*real(mu)的符号。如果乘积为负说明 n 的符号和解不一致直接把 m 整体减 1重新计算 n、ε、μ。更严谨的做法是回到麦克斯韦方程组的因果性约束但工程上这个“先反演、再判号、再平移一个周期”的流程足够处理大多数左手材料单胞也更容易在脚本里自动实现。5. 最后一公里批量扫描、结果验证与常见返工5.1 用 CST 宏和 MATLAB 循环做批量扫描超材料设计几乎都要扫参数方环宽度、开口间隙、介质厚度。每次手动在 CST 里改参数、导出.s2p、再跑 MATLAB会让人崩溃。CST 的 VBA 宏可以把整个流程录制下来我通常会录制一次仿真和导出操作然后把参数改成一个循环变量。下面是一个典型骨架不同 CST 版本方法名略有差异最好先宏录制一次再替换循环体Sub BatchScan() Dim app As Object Dim i As Integer Set app CSTStudio.Application For i 1 To 8 app.StoreParameter gap, 0.1 * i app.Rebuild app.Solver.Start app.ExportSParameterFile D:\sim\gap_ i .s2p Next i End Sub在 CST 中运行宏之前先确保模型参数名是gap并且已经在参数列表里定义。ExportSParameterFile的路径必须提前存在否则会静默失败。宏跑完后MATLAB 侧用dir批量读取并反演files dir(D:\sim\*.s2p); d 0.0016; for k 1:numel(files) sp read_s2p_ri(fullfile(files(k).folder, files(k).name)); [eps_r, mu_r, n, z] invert_nrw(sp.S11, sp.S21, sp.f, d); [~, tag] fileparts(files(k).name); save(sprintf(%s_invert.mat, tag), eps_r, mu_r, n, z, sp); end这里每轮循环会把原始 S 参数和反演结果一起存进.mat文件后续画损耗角正切、色散曲线都不用重新读.s2p。做参数扫描时我习惯把gap值写进文件名这样从文件名就能重建扫描条件和频点省掉一份额外记录表。5.2 反演结果自检拿 ε、μ 回填 CST最有效的验证不是看曲线是否平滑而是把反演得到的 ε_r、μ_r 作为均匀介质块的参数填回 CST再算一次 S 参数然后与原始超材料单胞的 S 参数对比。具体做法是在 CST 中新建一个尺寸相同的矩形块材料类型设为 Dispersive Material把某个频点反演出的复 ε、复 μ 填进去仿真同一频点如果 S21 幅度差超过 0.1、相位差超过 5°说明这个超材料单元在目标频带内不能等效为均匀介质需要缩小单元尺寸或重新检查亚波长条件。这是我在所有反演脚本里都会保留的一步。很多源码包不提供验证脚本导致你反演出的 ε 看起来再漂亮也可能只是分支选择“恰好连续”而已。如果回填后 S21 对不上优先检查的不是 MATLAB 代码而是 CST 里的 reference plane 有没有校准以及超材料单元是否真的满足周期边界下电尺寸小于 λ/4 的要求。最后提醒一个非常实际的操作批量扫描时不要让 MATLAB 和 CST 同时打开同一批.s2p文件。CST 导出文件时可能还没写完缓冲区MATLAB 读到的数据会缺行我吃过这个亏后会在 CST 宏里对每个导出文件做一次文件大小检查或者在 MATLAB 侧用pause(2)延迟读取确认文件大小不再变化后再进反演流程。本文还有配套的精品资源点击获取