阵列波导光栅AWG的Matlab仿真:从相位模型到FFT加速
简介面向光通信与集成光子方向的课程设计、毕业设计人群提供阵列波导光栅AWG模型及其在Matlab中的仿真完整项目。代码围绕材料折射率Si、SiO2、Si3N4等、平板波导模场分析、有效折射率计算与光谱响应等关键步骤展开包含AWGSimulator交互应用及相关示例脚本从基础工具函数到整体仿真流程均有实现。压缩包共22个文件以17个m源文件为主辅以mlx实时脚本、mlapp应用与说明文档打包体积1.07MB目录结构清晰便于按模块查阅和二次开发。已有173人学习下载。项目代码均测试运行成功适合作为课设答辩、毕设演示或项目初期原型参考通信、电子信息等专业学生可在此基础上修改拓展。1. 阵列波导光栅AWG模型及其在Matlab中的仿真到底要做什么拿到这个题目我第一反应不是去找一个现成的AWG工具箱而是先把“AWG仿真”到底算到哪一步搞清楚。AWG属于平面光波导器件做工艺级仿真要动到FDTD或BPM写在毕设题里通常没那么深。真正会被要求输出的是给定中心波长和信道数设计出一组阵列波导参数在Matlab里扫波长得到不同输出波导的透射谱再读FSR、插损、串扰和带宽这组器件指标。也就是说模型在物理上只做“相位叠加”不碰麦克斯韦方程组的数值解。这个定位想清楚后面代码就写得很快。2. AWG模型的第一性原理相位延迟、衍射级次与FSR计算2.1 为什么AWG靠相位差就能分波AWG结构上可以拆成三段输入自由传播区、阵列波导区、输出自由传播区。光从输入波导出来后在自由传播区自然发散照到一组长度递增的阵列波导上。每条阵列波导的长度相差ΔL光在里面走的时间就不同于是到达输出自由传播区时相邻波导之间带上一个k0·n_eff·ΔL的相位延迟。输出自由传播区再次把各通道的光重叠在一起发生干涉。某个波长如果正好让阵列波导间的相位差是2π的整数倍就会在对应的输出波导位置形成相干增强其他位置相互抵消这就是AWG分波的本质。如果把输出位置的角度记为θ干涉极大条件可以写成n_eff·ΔL n_s·d·sinθ m·λ这个式子很关键。它说明三件事第一θ一旦确定中心波长就由m和ΔL唯一决定第二同一级次m下不同θ对应不同λ所以不同输出口天然对应不同波长第三m是衍射级次不是波导编号。在Matlab建模时我会先把这条式子写成独立的function后面所有参数设计都从它反推。2.2 用中心波长和衍射级次反推ΔL设计AWG的第一步是给定中心波长λ0选一个衍射级次m然后把相邻阵列波导的长度差ΔL算出来。原因很简单ΔL决定FSR而FSR要大于整个波分复用带宽不然相邻级次的光会串进来。反过来ΔL又不能太小否则m太小器件看起来更像普通衍射光栅相位延迟优势体现不出来。在Matlab里这个计算只需要几行lambda0 1.55e-6; % 中心波长单位 m m 85; % 衍射级次无量纲 n_eff 2.2; % 阵列波导的有效折射率 DeltaL m * lambda0 / n_eff; fprintf(DeltaL %.3f um\n, DeltaL * 1e6);这里的m85是随手先给的值需要验证。把输入数字代入DeltaL约为59.9μm恰好落在当前硅光或二氧化硅AWG常见的几十微米量级。注意fprintf里乘1e6是把米换成微米输出单位要单独标清楚不然之后在图上读到参数时容易把1e-6和1e-9弄混。2.3 FSR与信道间隔的取舍FSR只和中心波长、群折射率、ΔL有关FSR λ0² / (n_g·ΔL)其中n_g是群折射率它把材料色散也囊括进去数值上比n_eff略高。用上一小节的ΔL59.9μm、n_g2.5估算FSR约等于16nm。这个数值意味着一个器件最多承载约20个100GHz信道因为100GHz信道间隔对应0.8nm。我把这个取舍做成一张小表方便在不同毕设要求间切换衍射级次mΔL (μm)FSR (nm)可容纳100GHz信道数5035.227.3348559.916.12012084.511.414写这段表的时候要提醒一句m越大阵列波导越长工艺容差越差但FSR越小m越小FSR越大但相邻信道的相位差变小对自由传播区的尺寸要求更高。大部分毕设里m取60到100之间是稳妥的。若题目里指定了信道数比如32信道先按FSR≥1.2倍总带宽去倒推ΔL再回算m这比拍脑袋选m要可靠。3. 在Matlab中把AWG拆成三段式模型坐标、相位与透射谱3.1 三段式输入自由传播区、阵列波导、输出自由传播区我在Matlab里建AWG不做严格的光场传播而是用“采样点叠加”的方式做数值近似。第一步给阵列波导一个横向坐标xx通常是对称分布例如从-N-1/2到N-1/2再乘上阵列间距d。输入波导的能量耦合进每条阵列波导的幅度用一条高斯包络来近似tail exp(-xx.^2 / w_a^2)w_a是光斑在阵列端展开后的半宽。这个近似在毕设精度下完全够用因为它抓住了AWG最重要的两点中心波导耦合最强、两侧逐渐变弱。如果题目要求均匀变迹把tail全部设为1即可但旁瓣会变高。输出口的场叠加则写为E_out sum(tail .* exp(1j * phi))其中phi是每条阵列波导在输出端贡献的总相位。把这条式子展开就是阵列波导延迟相位加上自由传播区由位置差带来的线性相位。摆清楚这三段在代码里的变量对应关系比背公式有用。3.2 相位构建最容易翻车的地方少算了线性相位不少学生搭AWG模型时只在phi里放n_eff·ΔL这一项输出口全用同一个位置去算结果谱线只有一条完全看不出波长和输出口的对应关系。问题在于“输出位置”没有进到相位里。实际从第k条阵列波导到某个输出口O多走的路径近似为d·x_O/f换成相位是k0·d·x_O/f这会直接决定峰值位置。在3.3的代码我会把这个线性项加进去。调试的时候如果发现频谱图左右不对称优先检查这一项的正负号。3.3 最小可运行代码无工具箱扫出透射谱下面这版代码是给毕业设计用的最小版本不依赖任何Matlab工具箱只需要一个.m脚本。我故意把循环写得更直观先保证能跑通性能优化放到后面一章。clear; clc; lambda0 1.55e-6; % 中心波长 1550 nm n_eff 2.2; % 有效折射率 n_g 2.5; % 群折射率 N 32; % 阵列波导数 d 15e-6; % 阵列波导间距 m 85; % 衍射级次 DeltaL m * lambda0 / n_eff; f 200e-6; % 自由传播区长度 w_a 20e-6; % 阵列端光斑半宽 xx ((1:N) - (N1)/2) * d; % 对称排布 tail exp(-xx.^2 / w_a^2); % 高斯变迹包络 lam linspace(1540e-9, 1560e-9, 2001); x_out linspace(-30e-6, 30e-6, 7); Spec zeros(length(x_out), length(lam)); for oi 1:length(x_out) for wi 1:length(lam) k0 2 * pi / lam(wi); phi k0 * n_eff * (0:N-1) * DeltaL ... k0 * xx * x_out(oi) / f; E sum(tail .* exp(1j * phi)); Spec(oi, wi) abs(E).^2; end end Spec Spec / max(Spec(:)); norm Spec(:, lam 1.55e-6); % 中心波长处的空间分布这段代码的循环逻辑值得解释外循环固定输出口位置内循环扫描波长每一轮都重新计算k0等于把n_eff随波长的缓慢变化留给了k0去吸收适合只分析谱线形状的场景。如果后面要精确算材料色散再引入Sellmeier方程。参数方面f200μm决定线性相位的曲率f太小会让输出口的位置分辨率变差d决定阵列波导对空间频率的采样周期d选得太大会让离散栅瓣进入输出区tail的w_a控制旁瓣。这四个参数加上N、m就是整个模型的可调面。3.4 参数表变量、单位、在模型里干什么整理成一个参数表作为仿真报告里的表格模板参数符号典型值作用调大后的效果阵列波导数N32参与叠加的波导数量串扰下降计算量上升阵列间距d15 μm空间采样周期栅瓣位置外移但器件尺寸变大衍射级次m85决定ΔL和FSRFSR变小传播区长度f200 μm线性相位缩放输出位置分辨率变化光斑半宽w_a20 μm高斯包络宽度越大旁瓣越高长度差ΔL59.9 μm波导延迟由m和λ0确定这个表在写毕业设计文档时可以直接复制到需求分析或者参数设计小节比放一屏代码更直观。4. AWG仿真核心参数调优串扰、带宽与“发散”排查4.1 从透射谱里读四个指标中心波长、FSR、串扰、3dB带宽仿真做出来不是跑完就结束要能回答“这个AWG参数好不好”。我会在Matlab里写一个评估函数输入是Spec矩阵和lam数组输出四个指标。中心波长指的是每个输出口峰值所在的波长值用find(Spec(oi,:) max(...))就能找到。FSR是相邻两个峰值间的波长差直接在谱线上读两个峰之间的距离。串扰是目标口主峰功率与它旁边通道在同一波长的功率之比在线性坐标下是比值在图上一般用dB显示为负值。3dB带宽则是主峰两侧下降到一半处对应的波长宽度反映通道对波长偏移的容忍度。这四个指标从直觉上就能对应研究目标FSR决定能容纳多少信道串扰决定相邻信道隔不隔离带宽决定对激光器波长漂移的敏感度。毕业设计回答“我设计的AWG性能如何”时至少要有这四个数字。4.2 三个必调参数N、m、高斯基底宽度先看N。阵列波导数N从16加到80透射谱的旁瓣会逐渐降低。原因是从“有限长度光栅”角度看N就是参与叠加的周期数周期数越多干涉越锐利。但N也不是无限加大就好仿真时间线性增长而且到了某个值以后旁瓣降幅趋缓收益变薄。我在毕设里一般会扫N16、32、64三档把串扰变化列成表选一个“再往上加也不明显”的点。再看m。m增大导致FSR变小信道间隔不变的情况下可容纳信道数变少。反过来m太小会让ΔL变到几十微米以下相位延迟不够输出口难以分辨相邻波长。所以m的调法通常是给定信道间隔后反推而不是盲目试。最后是高斯基底宽度w_a。w_a越小高斯包络越尖旁瓣压低但主峰变宽w_a越大越接近均匀阵列主峰更窄第一旁瓣逼近-13dB。这个参数是“带宽与串扰之间的权衡旋钮”。调整w_a时有人会误把w_a设成和d一个量级导致除了中心波导其余全被衰减掉谱线变成一个大包络。我习惯让w_a落在d的1到2倍。4.3 仿真发散最常见的三个原因及其排查Matlab仿真发散常常不是物理发散而是数值处理出了岔子。第一个是相位缠绕累积。当扫描波长范围跨过大半个FSRphi里的k0·ΔL项反复超过2π从-π跳到π。在求和之前应该对每个输出口的相位做unwrap不然叠加结果里出现毛刺。第二个是波长扫描步长太粗。100GHz信道间隔是0.8nm如果步长取1nm两个峰之间采不到点频谱图直接“漏峰”看起来像发散。我会把步长压到0.02nm以下再跑一遍对比峰值是否稳定。第三个是log10(0)带来的-inf把整张图拉扁。在画对数坐标时Spec里有零值时先做Spec(Spec0)eps再取log10。这三条排查顺序是我踩了几次坑后总结的先看是否存在-m再看是否存在0值最后才怀疑物理模型本身。因为数值问题出现的概率远大于模型假设错误。4.4 用一段脚本直接读结果并标注峰值调完参数后我用下面这段脚本把峰值位置和串扰直接打印出来省得每次都在图上肉眼读数for oi 1:size(Spec, 1) [~, pi_] max(Spec(oi, :)); fprintf(输出口%02d 峰值波长 %.2f nm\n, ... oi, lam(pi_) * 1e9); end % 中心波长处的空间串扰 [~, ci] max(Spec); % 每波长下最强输出口 crosstalk_dB 10 * log10(Spec(:, ci(1)) ./ Spec(:, ci(2)));打印结果时如果发现相邻输出口的峰值波长差不等于0.8nm先怀疑x_out的采样范围是否覆盖了整个成像区再检查f是否偏小。相比之下直接打印数值比分次看图要可靠这也是把仿真结论写进报告最省事的方式。5. 验证与提速用FFT把AWG仿真加速上百倍5.1 把输出口的叠加写成傅里叶变换回到3.3里的叠加式E sum(tail·exp(j·k0·n_eff·(0:N-1)·ΔL)·exp(j·k0·xx·x_out/f))。等一下第一项exp(j·k0·n_eff·(0:N-1)·ΔL)不依赖x_out对整个输出分布是常数因子第二项里的k0·xx·x_out/f正是xx随x_out的线性相位。傅里叶变换恰好就是把一个序列乘上线性相位再求和。因此一次FFT就能得到所有输出口的场强而不必对x_out循环。这个改写带来的收益在N大到64、输出口有几十个时非常明显。原本是两层循环嵌套现在变成一次fft配合矩阵乘法处理波长扫描速度能提升一到两个量级。5.2 用FFT版本验证原来的逐点叠加下面的脚本以3.3节的模型为基础用FFT计算单个波长下的输出分布再和循环版对比k0 2 * pi / lambda0; phi k0 * n_eff * (0:N-1) * DeltaL; seq tail .* exp(1j * phi); % 直接对序列做FFT得到输出分布 Spec_fft fftshift(ifft(fftshift(seq))); % 对应的输出位置坐标 x_fft (-N/2:N/2-1) * (lambda0 * f / (N * d)); % 和逐点叠加对比 Spec_loop zeros(1, N); for oi 1:N xo (oi - N/2) * d * 0.1; ph k0 * xx * xo / f; Spec_loop(oi) abs(sum(tail .* exp(1j * (phi ph)))).^2; end err max(abs(Spec_fft).^2 - Spec_loop); fprintf(max|diff| %e\n, err);注意这里用ifft还是fft取决于x_out的符号约定。我习惯把输出位置坐标和空间频率对应好x_fft的单位由λ·f/(N·d)给出可以借此直接换算输出口的物理间距。err量级在1e-14附近说明实现正确。5.3 三个在毕业设计里能落地的验证技巧一旦FFT版本跑通我通常会做这三件事一是把FFT结果和逐点叠加画在同一张图里重叠曲线直接作为“代码正确性验证”的截图二是把波长扫描也向量化用lam作为其中一个维度一次性得到完整的Spec矩阵画成彩色surf图能直接看到波长与输出口的映射关系三是把FSR实测值和理论值并列打印差值不超过2%就可以在答辩时说模型验证完毕。需要留意的是FFT版在输出口数量等于N的整数倍时会丢失某些中间位置的细节所以严格验证还是保留3.3的循环版作为对照。用一个开关变量控制跑哪一版例如用use_ffttrue切到FFT分支这样既能提速又能随时回到原始实现自查。命令行里跑一下err如果出现1e-10以下的差值那代码和物理模型之间的一致性已经有了闭环证据接下来画图、写报告基本不用返工。本文还有配套的精品资源点击获取