简介基于Matlab的MRR微环谐振器仿真图源码资料面向光通信、电子信息、数学等专业需要完成课程设计或毕业设计的学生。资源聚焦微环半径、透射系数与耦合系数对传输特性的影响覆盖普通微环半径5~30μm、t0.96/0.963、k0.01~0.5等典型参数组合可帮助读者理解微环滤波原理并快速搭建类似仿真。压缩包内共7个文件以fig格式的仿真结果图为主附1个TypicalMicroRing.m主程序及txt说明文档便于对照代码与图形输出整体仅315KB轻量实用。解压后目录结构简单便于快速定位文件。通过阅读主程序可学习参数扫描循环、绘图样式控制与图像保存等实用Matlab技巧方便迁移到其他光学器件仿真中。已有617人学习查看适合有一定Matlab基础、希望参考完整绘图思路的读者作为课设或毕设的参考资料。1. 基于Matlab实现MRR-微环仿真图先搞清楚这条谱线到底在说什么第一次用Matlab画MRR微环仿真图我猜你不是卡在物理上而是卡在“为什么我算出来的透射谱长这样”。MRR即微环谐振器是硅光、集成光电子里最常用的滤波与调制结构一根直波导和一个环形波导靠在一起满足共振条件的波长会陷进环里直通端口出现一连串陷波谷。所谓仿真图核心就是这条透射谱。我用传输矩阵法从头写不依赖任何光子学工具箱只靠Matlab自带语法跑通适合课程设计、毕业设计、流片前快速验证的工程师和学生。全套代码几十行改半径、耦合系数、损耗就能立刻看到谱线怎么变这也是我推荐先拿它练手的原因。2. 先把物理模型立住微环谐振器的传输矩阵与三个核心参数2.1 微环为什么会出现共振谷相位条件与自由光谱范围光在环形波导里绕一圈相位积累量是 βL其中 β 是传播常数L 是环周长。当这个相位恰好是 2π 的整数倍时环内形成相长干涉能量不断在环里增强直通端的输出就会掉下去形成一个共振谷。用波长写就是2π × neff × L / λ 2πm化简后得到 neff × L m × λ。这里的 neff 是模式等效折射率m 是整数。换句话说只有特定的几个波长能满足这个条件所以透射谱上不是只掉一个谷而是一串等间隔的谷。相邻两个谷之间的距离叫自由光谱范围即 FSR。一阶近似下FSR ≈ λ² / (ng × L)注意这里必须用群折射率 ng而不是相折射率 neff。ng neff − λ × dneff/dλ它跟模式色散有关。拿一个常见尺寸来算半径 R 10 μm环周长 L ≈ 62.83 μmng 4.2工作波长 1.55 μm那么 FSR ≈ 9.1 nm。这个数字在设计里非常重要它决定了这个环在一个波段内能放下几个共振峰也决定了热调谐或者电调谐能把峰推多远。很多人在这一步就会翻车用 neff 去算 FSR结果算出来偏得离谱而且画出来的谱线周期也对不上。原因后面避坑章会展开这里先记住一个结论仿真谱线的横轴周期由 ng 决定不是由 neff 决定。2.2 用传输矩阵把直波导和环形波导写进 Matlab微环仿真的经典做法是传输矩阵法也叫散射矩阵法。我不在这里堆公式只讲怎么落地成 Matlab 代码。直波导和环形波导之间有一个耦合区。这个耦合区可以用两个参数描述自耦合系数 t 和交叉耦合系数 κ。物理上|t|² 是直通端直接过去的功率比例|κ|² 是耦合进环的功率比例。无损耗耦合时两者满足|t|² |κ|² 1光在环内走一圈还会遇到波导损耗。设单圈往返的场传输幅度为 a它与功率损耗系数 α 的关系是a exp(−α × L / 2)之所以是 L/2因为场幅度衰减是功率衰减的一半a 描述的是电场幅度。把耦合区和环内传播串起来可以得到全通微环直通端的功率传输函数T (t² a² − 2at·cosφ) / (1 a²t² − 2at·cosφ)其中 φ 2π × neff × L / λ 是单圈相位。这个公式非常关键整个全通微环的仿真图都是从它来的。当 t a 时分子在共振点φ 2πm处恰好等于 0也就是消光比趋向无穷大物理学上叫临界耦合。t a 是欠耦合消光比变小t a 是过耦合谷底会被抬起来。2.3 耦合系数、损耗系数与群折射率一张参数表和一个选参套路我先给出一组在硅光 SOI 平台、1550 nm 波段常见的参数量级后面所有代码都用这套量级。它不是某个流片厂家的精确值但足够让仿真趋势和实物一致。参数符号典型量级怎么确定环半径R520 μm由 FSR 和工艺最小弯曲半径决定等效折射率neff2.32.6用模式求解器或经验公式取近似值群折射率ng3.84.5用色散数据拟合或者由 FSR 实测反推单圈场衰减a0.980.999由波导损耗 α 和周长决定自耦合系数t0.80.999由耦合区 gap、长度决定可调我一般会先按量级设一组参数跑通再往实测靠。比如第一版经常用 R 10 μm、neff 2.4、ng 4.2、a 0.995、t 0.995这样 t 和 a 刚好相等能看到最漂亮的临界耦合谷。a 0.995 对应的波导损耗大约是 7 dB/cm对 SOI 平台来说是一个偏保守但完全合理的量级。把这组参数写进 Matlab 之前先统一单位。波长用米半径用米损耗系数用 1/m。下面这个脚本就是参数初始化% mrr_params.m % MRR 微环仿真图参数初始化SOI 平台近似参数 R 10e-6; % 环半径单位米 L 2 * pi * R; % 环周长单位米 neff 2.4; % 1550 nm 附近等效折射率 ng 4.2; % 群折射率决定 FSR 周期 dneff_dT 1.8e-4; % 热光系数单位1/K后面对热调谐有用 lambda linspace(1.50e-6, 1.58e-6, 200001); % 扫描波长范围单位米 a 0.995; % 单圈场传输幅度对应约 7 dB/cm 波导损耗 t 0.995; % 耦合区自耦合系数ta 为临界耦合这段代码里最关键的是lambda的密度。200001 个点看起来吓人但 Matlab 对向量化运算非常快几毫秒就能算完。点数太少的话共振谷的尖峰和半高宽会被采样磨平后面提取 Q 值时会发现问题。a和t我故意取成相等是为了第一版就能看到临界耦合谷如果你想要更接近实物通常 t 会和 a 有一个小的失配消光比是有限的。参数说明就一句话所有单位统一到米和 1/m别用纳米和 dB/cm 混着算这是微环仿真里最常见的低级错误。3. 用 Matlab 跑通第一张 MRR 透射谱最小可复现脚本3.1 主脚本从波长扫描到 plot 出图这一节我把最小可复现脚本拆成两部分主脚本负责扫描波长、调用核心函数、出图核心函数负责传输矩阵的计算。好处是后面做参数扫描时不需要改动主逻辑只需要在函数外面循环。先看主脚本% mrr_demo.m % 全通 MRR 微环仿真图临界/欠/过耦合三条曲线对比 lambda linspace(1.50e-6, 1.58e-6, 200001); R 10e-6; neff 2.4; a 0.995; T_crit mrr_thru(lambda, R, neff, a, 0.995); % 临界耦合 ta T_under mrr_thru(lambda, R, neff, a, 0.999); % 欠耦合 ta T_over mrr_thru(lambda, R, neff, a, 0.950); % 过耦合 ta lambda_nm lambda * 1e9; % 转成纳米画图易读 plot(lambda_nm, 10*log10(max(T_crit, 1e-12)), LineWidth, 1.2); hold on; plot(lambda_nm, 10*log10(max(T_under, 1e-12)), LineWidth, 1.2); plot(lambda_nm, 10*log10(max(T_over, 1e-12)), LineWidth, 1.2); xlabel(Wavelength (nm)); ylabel(Transmission (dB)); legend(Critical ta, Under ta, Over ta, Location, southwest); grid on; hold off;这段代码的逻辑很直白定义波长扫描向量调用自写函数生成透射率再用10*log10转成 dB。注意max(T_crit, 1e-12)这个细节透射率在临界耦合的谷底理论上是 0直接取log10会得到-Inf曲线画出来会掉到图底干扰视角。加一个 1e-12 的下限相当于把谷底限制在 −120 dB既不会破坏曲线形状又避免无穷大问题。参数说明T_under用了 t 0.999比 a 大一点点T_over用了 t 0.950比 a 小不少。这样三条曲线在共振谷附近的表现差异会很明显。lambda_nm只是单位换算不改变物理结果。3.2 传输矩阵计算函数手写而不是调工具箱核心函数mrr_thru就是上一章传输函数的直接翻译function T mrr_thru(lambda, R, neff, a, t) % MRR 全通微环直通端透射率 % lambda : 波长向量单位米 % R : 环半径单位米 % neff : 等效折射率近似值 % a : 单圈往返场传输幅度 % t : 直波导-环耦合区自耦合系数 L 2 * pi * R; phi 2 * pi * neff * L ./ lambda; % 单圈相位向量 T (t.^2 a.^2 - 2*a*t.*cos(phi)) ... ./ (1 a.^2 * t.^2 - 2*a*t.*cos(phi)); % 透射功率 end这里用到的是./和.*这种点运算因为lambda是向量phi是向量必须逐元素计算。很多新手第一次写微环仿真就是栽在没加点号上Matlab 要么报矩阵维度错误要么算出一个莫名其妙的标量。参数说明neff我在这里取常数这是一个刻意简化。严格来说 neff 随波长变化但在 1550 nm 附近一个 80 nm 的扫描窗口里neff 的变化量只有百分之几对共振谷的位置和形状影响不大。如果你想更精确可以把 neff 做成波长查表或者用一个线性色散模型代替常数这个我在 4.2 节讲热调谐时会顺带说。还有一个容易被忽略的点phi的分母用的是lambda本身所以共振谷在短波方向会比长波方向稍微密集一点这是正常的色散效应不是代码 bug。3.3 先看结果临界耦合、欠耦合、过耦合三条曲线跑完上面的脚本你应该看到一组周期性的透射谱每个周期内有一个向下的尖峰相邻尖峰间距大约 9.1 nm和前面用 FSR 公式算的一致。三条曲线的差异集中在共振谷附近临界耦合曲线 t a 0.995谷底能掉到 −80 dB 以下看起来非常锋利。欠耦合曲线 t 0.999谷底变浅消光比可能只有 20 dB 左右但谷的宽度比临界耦合时更窄。过耦合曲线 t 0.95谷底被明显抬高消光比进一步变差而且谷的形状变得宽而平。这三个状态在实际器件里都有意义。临界耦合是做滤波器、调制器最想要的消光比高欠耦合 Q 值高适合做传感和慢光过耦合在接收机前端有时用来展宽带宽。所以仿真图的目的不只是画一条漂亮的线而是帮你看出 t 和 a 的匹配关系。拿到一组实测透射谱如果谷底不够深你第一反应应该是t 和 a 不匹配要么是耦合间隙不对要么是波导损耗比设计值大。4. 把仿真图做得更接近实物损耗、调制与上下路微环4.1 加波导损耗和弯曲损耗消光比不再无限大上一章把 a 直接设成了常数算起来方便但如果你要跟实测对照最好把 a 跟波导损耗关联起来。常见做法是先用 dB/cm 描述损耗再换算成单圈场衰减 a% mrr_loss_scan.m % 固定自耦合系数 t扫描波导损耗观察消光比变化 lambda linspace(1.50e-6, 1.58e-6, 200001); R 10e-6; L 2*pi*R; neff 2.4; t 0.995; alpha_dB_list [2, 6, 15]; % 损耗依次增大单位 dB/cm for i 1:numel(alpha_dB_list) alpha_dB alpha_dB_list(i); alpha alpha_dB * log(10) / 10 / 0.01; % dB/cm - 1/m a_i exp(-alpha * L / 2); % 单圈场衰减 T_i mrr_thru(lambda, R, neff, a_i, t); plot(lambda*1e9, 10*log10(max(T_i,1e-12)), LineWidth, 1.2); hold on; end xlabel(Wavelength (nm)); ylabel(Transmission (dB)); legend(2 dB/cm, 6 dB/cm, 15 dB/cm); grid on; hold off;注意里面的单位换算。alpha_dB * log(10) / 10是把 dB/cm 转成每厘米的功率衰减系数再除以 0.01 才是每米的功率衰减系数最后在exp(-alpha * L / 2)里除以 2是因为我们关心的是电场幅度衰减。这一套换算我每次写都要在心里过一遍因为它太容易错有人忘了除 0.01有人忘了除 2结果算出来的 a 要么等于 1要么小到离谱。参数说明当 t 0.995 时损耗 2 dB/cm 对应 a ≈ 0.9985属于欠耦合消光比有限损耗 6 dB/cm 对应 a ≈ 0.9955接近临界谷底非常深损耗 15 dB/cm 对应 a ≈ 0.9889进入过耦合谷底抬升明显。这组脚本很适合用来做设计预算给定耦合间隙决定 t你能一眼看出波导损耗允许多大才能达到目标消光比。4.2 热调谐与电调谐共振波长漂移的 Matlab 实现微环最常用的功能是调谐。热调谐的原理是硅的热光效应温度变化改变材料折射率进而改变 neff共振波长随之移动。仿真里不需要真的建热模型直接把 neff 加上一个温度相关的变化量就行% mrr_thermal_tuning.m % 热调谐仿真不同温度增量下透射谱偏移 lambda linspace(1.50e-6, 1.58e-6, 200001); R 10e-6; L 2*pi*R; neff0 2.4; a 0.995; t 0.995; dneff_dT 1.8e-4; % 热光系数1/K dT_list [0, 10, 20, 40]; % 温度增量单位 K for i 1:numel(dT_list) neff neff0 dneff_dT * dT_list(i); % 只改折射率 phi 2 * pi * neff * L ./ lambda; T (t.^2 a.^2 - 2*a*t.*cos(phi)) ... ./ (1 a^2*t^2 - 2*a*t.*cos(phi)); plot(lambda*1e9, 10*log10(max(T,1e-12)), LineWidth, 1.2); hold on; end xlabel(Wavelength (nm)); ylabel(Transmission (dB)); legend(dT0K, dT10K, dT20K, dT40K); grid on; hold off;这段代码里热调谐只体现在neff neff0 dneff_dT * dT_list(i)这一行。因为共振条件 neff × L m × λneff 变大λ 也变大所以所有谷整体向长波方向平移。近似算一下每 10 K 温度增量大约平移 0.66 nm这个量级跟硅光流片实测很接近。如果你做的是电调谐也就是通过载流子注入改变折射率逻辑完全一样只需要把dneff_dT * dT换成电致折射率变化量dneff_electrical。在仿真图层面这两种调谐只是改一个数而已真正复杂的反而在后面的调制带宽和损耗分析。参数说明lambda扫描范围和点密不用改因为调谐只是平移不会改变曲线的包络周期。如果你发现温度高了以后谱线周期也变了那不是代码错而是neff变化导致 FSR 略微变化这个二阶效应在 40 K 范围内可以忽略。4.3 从全通到 add-drop透射谱和 Drop 端谱线一起画全通微环只有一个直波导实际系统里更常见的是 add-drop 微环上下两根直波导谐振波长从输入端进去直接从 Drop 端出来直通端则空出这个波长。Add-drop 结构可以同时看两条谱线适合做滤波器和波分复用解复用器。function [T_thru, T_drop] mrr_add_drop(lambda, R, neff, a, t1, t2) % Add-drop 微环返回直通端与 Drop 端透射率 % t1, t2 : 上下两个耦合区的自耦合系数 L 2 * pi * R; k1 sqrt(1 - t1.^2); % 输入耦合区交叉系数 k2 sqrt(1 - t2.^2); % 输出耦合区交叉系数 phi 2 * pi * neff * L ./ lambda; % 单圈相位 z exp(1j * phi); z_half exp(1j * phi / 2); % 半圈相位 % 直通端传输 thru (t1 - t2 .* a .* z) ./ (1 - t1 .* t2 .* a .* z); % Drop 端传输经过输入耦合、半圈、输出耦合 drop (k1 .* k2 .* sqrt(a) .* z_half) ./ (1 - t1 .* t2 .* a .* z); T_thru abs(thru).^2; T_drop abs(drop).^2; end这里用到了复数传输。z exp(1j * phi)表示一圈的相位积累sqrt(a)是半圈的场衰减z_half是半圈相位。分母1 - t1*t2*a*z描述的是环内多圈叠加是传输矩阵法处理环腔的标准形式。注意a不是平方关系因为在分母里 a 表示整圈的场传输幅度这个约定和全通函数里保持一致。调用这个函数并画图% mrr_add_drop_demo.m lambda linspace(1.50e-6, 1.58e-6, 200001); [Thru, Drop] mrr_add_drop(lambda, 10e-6, 2.4, 0.995, 0.985, 0.985); subplot(2,1,1); plot(lambda*1e9, 10*log10(max(Thru,1e-12)), LineWidth, 1.2); ylabel(Thru (dB)); grid on; title(Add-drop MRR); subplot(2,1,2); plot(lambda*1e9, 10*log10(max(Drop,1e-12)), LineWidth, 1.2); xlabel(Wavelength (nm)); ylabel(Drop (dB)); grid on;这段代码里我取 t1 t2 0.985表示上下耦合区对称。在共振波长处Drop 端接近 1直通端接近 0两个图形成互补。如果 t1 和 t2 不相等Drop 端的峰值会下降直通端的谷底也会抬起来。实际设计 add-drop 滤波器时一个经典的调参套路是先定 t2 决定 Drop 端带宽再定 t1 平衡端口反射和通道串扰。5. MRR 仿真避坑手册5 个让新手翻车的细节5.1 横轴用波长还是频率FSR 周期性为什么看起来不对现象用波长当横轴画出来的透射谱相邻谷的间距看起来忽大忽小怎么跟公式对不上。原因FSR 在频率轴上是严格均匀的在波长轴上只是近似均匀。微环共振条件是频率的等间隔序列映射到波长轴后因为 λ c/f同样的频率间隔在不同波长处对应的波长间隔不一样。再加上 neff 本身有色散波长轴上的不均匀性会更明显。解决如果要证明仿真结果和理论一致先用频率轴画一次确认峰间隔均匀再切回波长轴画图用于论文展示。Matlab 里切换很简单把lambda换成c ./ lambda即可。检查的时候不要肉眼盯间距直接在数据里取相邻两个谷的坐标用 diff 看一下差值的相对波动通常 80 nm 扫描范围内波动不超过 1%就说明模型没问题。5.2 群折射率写错相邻谐振谷间距系统性偏差现象透射谱画出来了形状很好看但谷与谷之间的间距比理论值大了一倍或者整体缩水。原因FSR 公式用的是群折射率 ng不是等效折射率 neff。很多教程一开始只提 neff新手就顺手用 neff 去算。SOI 波导在 1550 nm 附近 neff 大约 2.4ng 大约 4.2差得非常多直接导致 FSR 计算和仿真对不上。解决仿真函数里可以用 neff 做共振位置计算但做 FSR 验证、提取实验数据时必须用 ng。没有 ng 数据时一个靠谱的替代方案是用两个相邻谷的实测位置反推 ngng λ² / (FSR × L)。先测后算比拍脑袋准得多。5.3 临界耦合时消光比算不出“无穷大”数值精度和损耗清零的坑现象设定 t a 后理论上谷底透射率应该是 0消光比无穷大但仿真出来的谷底只有 −60 dB 或者 −80 dB怎么压都压不下去。原因Matlab 用的是双精度浮点两个几乎相等的大数相减时有效位数被截断结果达不到真正的 0。另一个常见操作是直接把 a 设为 1 来“模拟零损耗”这会让分母和分子在共振处都趋于 0数值上反而更不稳定。解决不要追求谷底绝对深度把注意力放在谷底相对变化上。画图时用max(T, 1e-12)截断能让视觉上保持整洁。如果你确实需要验证临界耦合的物理极限可以把 t 和 a 都改成精确值比如 t 0.995、a 0.995然后在共振波长处单独用符号运算或者高精度定点计算验证不要指望普通 plot 曲线的数值能体现无穷大。5.4 导出图片时线宽和字体全变了仿真图的最后一公里现象在 Matlab 里看曲线很清晰存成 PNG 插到论文或者课程报告里线变细、字变小、背景发灰甚至出现锯齿。原因很多人习惯用print -dpng或者直接截图这会使用屏幕渲染设置而不是出版级渲染设置。另外 Matlab 默认的白色背景在老的绘图脚本里可能没有显式设置存出来的图带着灰色背景。解决用exportgraphics。手动改一下效果会好非常多% 出版级导出示例 set(gcf, Color, w); % 白底 exportgraphics(gca, mrr_critical.png, Resolution, 300);这一步其实就是大家常说的 Matlab 图像处理流程的一部分从仿真数据到可投稿的矢量图先保证曲线属性再控制导出分辨率。300 dpi 是期刊最低要求600 dpi 更保险。如果还想进一步加工把数据存成 CSV 拖进 Origin 绘制也是常见做法但曲线形状的一致性要靠仿真数据本身足够密来保证。5.5 中文注释乱码与脚本编码Matlab 版本差异的坑现象代码里写了中文注释换个电脑打开工程注释变成乱码或者直接报错无法运行。原因Matlab 旧版本对 UTF-8 支持不完整Windows 系统下默认使用系统本地编码比如 GBK。代码文件用 UTF-8 保存在旧版 Matlab 里打开就会被错误解析。这个坑特别折腾尤其是把工程打包发给导师或者同事时。解决两个方案。方案一是所有注释和字符串统一用英文这是最保守的做法任何版本都不会出错。方案二是把脚本保存为 UTF-8并在文件开头加一行feature(DefaultCharacterSet, UTF-8)但旧版本未必认这个函数。我现在的习惯是自己用的临时脚本随便写中文注释凡是准备打包给别人或者上传的源码全部改成英文注释。反正代码逻辑不长英文注释成本不高。6. 让仿真图对得上实测数据Q 值提取与损耗反演技巧拿到一条实测透射谱以后第一件事不是看消光比而是提取 Q 值和 FSR。Q 值定义是共振波长除以 3 dB 带宽也就是Q λ0 / Δλ提取方法很简单先用findpeaks找谷的位置再量谷底的半高全宽。对临界耦合曲线这个操作在仿真里同样适用而且可以用来反推损耗。反推的思路是保持 t 不变扫描损耗参数 a让仿真 Q 值与实测 Q 值一致这时候的损耗就是波导实际损耗。这个操作我每次都用屡试不爽。% 从透射谱提取 Q 值示例 T_dB 10*log10(max(T_crit, 1e-12)); [~, idx] min(T_dB); % 共振谷位置 lambda0 lambda_nm(idx); half_level T_dB(idx) 3; % 谷底向上 3 dB % 在谷两侧找到与 half_level 最近的点取波长差 left find(T_dB(1:idx) half_level, 1, last); right find(T_dB(idx:end) half_level, 1, first) idx - 1; Delta_lambda lambda_nm(right) - lambda_nm(left); Q lambda0 / Delta_lambda;这段代码的局限是只适用于单谷提取。实际透射谱里可能多个谷重叠或者有噪声建议先用smoothdata做一次轻度平滑再用findpeaks(-T_dB, MinPeakProminence, 10)自动找谷。提取到 Q 和 FSR 之后核验两个量是否自洽FSR 反推 ngQ 反推损耗。如果这两个参数能闭环你的 MRR 微环仿真图才真正具备设计价值。我现在的习惯是拿到任何微环仿真结果都先跑一遍这个验证流程而不是直接看图形好不好看。这个顺序帮我挡掉了至少三次参数错乱希望帮到你。本文还有配套的精品资源点击获取