简介这套MATLAB光栅反射仿真代码包面向光学工程、光通信与光谱分析领域的研究者与学生涵盖均匀布拉格光栅、切趾光栅、相移布拉格光栅、取样光栅及啁啾光栅等多种类型可快速掌握不同光栅结构的反射谱仿真方法。压缩包共12个文件包含6个MATLAB脚本和6张仿真结果图整体仅88KB结果图直观展示反射谱/透射谱曲线便于对照验证。已有297人学习下载适合作为课程设计、论文仿真及项目预研的参考素材。通过调整代码中的光栅周期、折射率调制深度、占空比等参数读者能深入理解各光栅特性对反射谱的影响并基于MATLAB的FFT及光栅方程快速获得衍射强度分布为后续光学系统设计与优化奠定基础。1. 取样光栅反射谱仿真MATLAB里先要跨过的三道坎取样光栅不是把普通光栅的周期拉长而是对连续光栅做周期性的“开窗”采样每个采样窗口内有一段均匀光栅窗口之间是空白。这种结构在梳状滤波器、DFB激光器波长选择里很常见反射谱是一簇等间距的峰而不是单峰。在MATLAB里仿真这个结构看着像是改一改折射率数组的事但很多人会卡在三个地方取样周期和占空比如何映射到折射率分布传输矩阵法层数太多时计算慢到怀疑人生以及把多峰反射谱误判成数值误差。这篇文章直接从取样光栅的傅里叶展开切入给出一套可运行、可调参的MATLAB代码并说明参数扫描时怎么读谱、怎么验证传输矩阵实现没有写错。2. 取样光栅的折射率调制先用傅里叶展开算出梳状谱的间隔2.1 取样光栅的折射率分布模型取样光栅沿传播方向z的折射率分布最常见的写法是n(z) n_eff Δn · s(z) · cos(2πz/Λ)其中Λ是光栅小周期也就是每个窗口内条纹的周期Δn是窗口内的折射率调制深度s(z)是取样函数它是一个周期为P的方波窗口内为1窗口外为0。占空比duty定义为窗口内光栅长度与取样周期的比值即duty L_s / P。窗口外的区域也有n_eff只是没有Δn调制。在MATLAB里生成这个分布关键是判断z坐标是否落在取样窗口内。判断条件可以写成mod(z, P) dutyP。由于dutyP是窗口长度这个条件会把每个取样周期的前段变成窗口后段变成间隙。窗口起点在z0便于后续传输矩阵相位对齐。2.2 傅里叶展开如何决定反射峰的位置与幅度因为s(z)是周期为P的方波它可以展开成傅里叶级数。当s(z)乘以cos(2πz/Λ)时频率域产生卷积。s(z)的第m次谐波对应空间角频率2πm/P所以等效的周期调制波矢变成2π/Λ 2πm/P。这个谐波光栅的周期近似为1/Λ_m 1/Λ m/P第m个峰的布拉格波长近似为λ_m ≈ λ_B m · λ_B · Λ / P其中λ_B 2·n_eff·Λ。因此相邻谐波峰之间的波长间隔是Δλ ≈ λ_B · Λ / P这个间隔是梳状谱最基本的单位。方波取样函数的第m级傅里叶系数幅度为c_m duty · sinc(m·duty)其中sinc(x)sin(πx)/(πx)。不同占空比下前几级傅里叶系数的相对幅度如下表mduty0.25duty0.5duty0.7500.250.50.75±10.2250.3180.225±20.15900.159±30.0750.1060.075当duty0.5时除m0外的偶次谐波全部为0所以反射谱中每两个峰之间会缺一个。这一点会在后续仿真中直接看到。2.3 在MATLAB里生成取样光栅折射率曲线的代码用向量化方式生成折射率分布function nz sampled_grating_profile(n_eff, dn, Lambda, P, duty, L) % 生成取样光栅的折射率分布 % n_eff: 有效折射率 % dn: 折射率调制幅度 % Lambda: 光栅小周期 (nm) % P: 取样周期 (nm) % duty: 占空比 0~1 % L: 光栅总长度 (nm) % 输出 nz: 与坐标 z 等长的折射率数组 dz Lambda / 20; % 空间步长 z 0:dz:L; % 离散坐标 mask mod(z, P) duty * P; % 取样窗口掩码 nz n_eff * ones(size(z)); nz(mask) nz(mask) dn * cos(2 * pi * z(mask) / Lambda); end说明dz取Lambda/20是空间离散的常用选择它保证每个光栅周期内有20个采样点足够精确描述余弦调制。mask是一个逻辑数组为1的位置就是取样窗口。向量化避免了循环在L达到10mm、数组长度达到几十万时依然很快。调用示例n_eff 1.452; dn 1e-4; Lambda 534; % nm P 50e3; % 50 μm duty 0.5; L 10e6; % 10 mm nz sampled_grating_profile(n_eff, dn, Lambda, P, duty, L); plot((0:length(nz)-1)*dz/1e3, nz); xlabel(z (μm)); ylabel(折射率);运行后可以看到窗口内是余弦条纹窗口外是一条直线。如果窗口内包含的周期数不是整数边带谱会出现轻微不对称这是正常现象。3. 用传输矩阵法计算取样光栅反射谱的最小MATLAB实现3.1 传输矩阵法的离散模型传输矩阵法TMM把光栅沿轴向切成N段每段厚度为dz折射率近似为常数。每一层对应一个2×2矩阵M_i [ cos(k_i·dz), j·sin(k_i·dz)/n_i ; j·n_i·sin(k_i·dz), cos(k_i·dz) ]其中k_i 2π·n_i/λλ是真空波长n_i是第i层的平均折射率。整个光栅的总矩阵是各层矩阵按光传播方向依次相乘。要得到反射系数不能直接用M(2,1)/M(1,1)——那是只在特定边界条件下成立的简化。入射介质和出射介质都是n_eff时正确做法是计算vec M * [1; n_eff]; B vec(1); C vec(2); r (n_eff·B - C) / (n_eff·B C)反射率R |r|²。3.2 完整MATLAB函数采样周期、占空比、长度都可以调下面是计算取样光栅反射谱的完整函数它直接复用上一节的sampled_grating_profilefunction [wl, R] compute_sg_reflectance(n_eff, dn, Lambda, P, duty, L, wl_range) % 计算取样光栅反射谱 % 输入: % n_eff : 有效折射率 % dn : 折射率调制深度 % Lambda : 光栅小周期 (nm) % P : 取样周期 (nm) % duty : 占空比 0~1 % L : 光栅总长度 (nm) % wl_range : 波长范围 [wl_min, wl_max] (nm) % 输出: % wl : 波长序列 (nm) % R : 反射率 (0-1) dz Lambda / 20; nz sampled_grating_profile(n_eff, dn, Lambda, P, duty, L); N length(nz) - 1; nlambda 400; % 波长采样点数 wl linspace(wl_range(1), wl_range(2), nlambda); R zeros(size(wl)); for idx 1:nlambda lam wl(idx); M eye(2); for i 1:N ni 0.5 * (nz(i) nz(i1)); % 层内平均折射率 ki 2 * pi * ni / lam; Mi [cos(ki*dz), 1j*sin(ki*dz)/ni; 1j*ni*sin(ki*dz), cos(ki*dz)]; M M * Mi; end vec M * [1; n_eff]; B vec(1); C vec(2); r (n_eff*B - C) / (n_eff*B C); R(idx) abs(r)^2; end end代码里层间折射率取了相邻网格点的平均值这让折射率突变发生在层与层之间更接近物理模型。矩阵相乘顺序是从左到右即第1层在最靠近入射侧这个顺序一旦颠倒反射谱会出现镜像畸变。外层波长循环和内层分层循环是计算量的主要来源。当L10mm、dz26.7nm时N约37万400个波长点就是1.5亿次矩阵乘在普通PC上需要几十秒。可以先跑200个点看轮廓再在感兴趣的范围局部加密。调用示例[wl, R] compute_sg_reflectance(1.452, 1e-4, 534, 50e3, 0.5, 10e6, [1500,1600]); plot(wl, R, LineWidth, 1.2); xlabel(波长 (nm)); ylabel(反射率); grid on;3.3 运行结果怎么看反射谱中的多个峰值是否正常运行后会看到在1551.8nm附近有一个主峰两侧有多个次级峰。首先检查主峰位置2×1.452×534≈1551.8nm如果偏差超过零点几个纳米多半是单位换算错误。然后检查峰间距理论Δλ λ_B·Λ/P代入约1552×534/50000≈16.6nm。对于duty0.5偶次谐波消失所以实际可见的相邻峰间隔是2Δλ≈33.2nm而不是16.6nm。如果你看到16.6nm间隔说明占空比实际生效的不是0.5或者mask判断出现了错误。另一个需要注意的现象是反射谱顶部出现细密振荡。这种振荡通常来源于取样窗口边缘的菲涅耳反射而不是真正的布拉格反射。如果振荡幅度过大可以减小dz到Lambda/30或者增大吸收项。4. 取样光栅参数扫描周期、占空比、长度对反射谱的影响4.1 扫描取样周期P峰值间距如何变化取样周期P直接控制梳状谱的间隔。固定duty0.5分别取P20μm、50μm、100μm计算同一波长范围的反射谱。为了对比可以把三条曲线纵向偏移开P_scan [20e3, 50e3, 100e3]; % nm figure; hold on; colors lines(3); for idx 1:length(P_scan) [wl_i, R_i] compute_sg_reflectance(1.452, 1e-4, 534, P_scan(idx), 0.5, 10e6, [1500,1600]); plot(wl_i, R_i idx*0.3, Color, colors(idx,:)); end legend(P20um, P50um, P100um); xlabel(波长 (nm)); ylabel(反射率纵向偏移); grid on;P50μm时Δλ16.6nmduty0.5的可见峰间隔为33.2nm所以在1500-1600nm范围内能看到主峰附近两个一级边带和两个三级边带。P20μm时Δλ41.5nm可见峰间隔83nm一个100nm的窗口里往往只看到主峰和一个边带。P100μm时可见峰间隔16.6nm峰很多但高阶峰幅度衰减很快。注意主峰反射率理论上与P无关。因为取样窗口总长度是duty×L与P无关。如果你的仿真结果显示主峰反射率随P明显变化说明窗口内的相位对齐出了问题通常是取样窗口起点没有统一对齐或者空间步长dz没有整除取样周期。4.2 扫描占空比反射强度与包络的变化占空比决定傅里叶系数的分布。用下面的脚本扫描duty0.25、0.5、0.75duty_scan [0.25, 0.5, 0.75]; figure; hold on; for idx 1:3 [wl_i, R_i] compute_sg_reflectance(1.452, 1e-4, 534, 50e3, duty_scan(idx), 10e6, [1540,1580]); plot(wl_i, R_i, LineWidth, 1.2); end legend(duty0.25, duty0.5, duty0.75); xlabel(波长 (nm)); ylabel(反射率); grid on;运行后可以看到三个规律。第一主峰反射率随duty增大而增大因为窗口内总调制长度变大。第二duty0.5时偶次傅里叶分量为0所以只有m±1、±3、±5等奇数峰主峰与最近边带的距离是Δλ但相邻可见峰之间的间隔是2Δλ。第三duty0.25时偶次分量都存在峰密度更高但高阶峰幅度低。duty0.75时偶次分量不为0但幅度较小谱线形态介于前两者之间。4.3 用表格汇总扫描结果来定位设计参数把上面分析整理成表设计阶段可以直接查参数配置Δλ (nm)出现的m相邻可见峰间隔 (nm)主峰反射率(约)P20μm, duty0.541.51,3,5...830.58P50μm, duty0.516.61,3,5...33.20.58P100μm, duty0.58.31,3,5...16.60.58duty0.25, P50μm16.61,2,3...16.60.11duty0.75, P50μm16.61,2,3...16.60.68这里主峰反射率是近似值实际会因窗口内周期数是否整数而出现几个百分点的波动。表格给出一条清晰的设计路径要增加梳状谱间隔减小P要抑制偶数级峰用duty0.5要抬高主峰反射率增大duty但这会同时增强其他级峰边带抑制比会下降。实际项目中通常先确定Δλ反推P再调整duty平衡峰间均匀性。5. 验证传输矩阵实现的一个技巧把取样光栅退化到均匀光栅传输矩阵法代码写完最怕矩阵乘法顺序或反射系数公式出错。一个有效的验证手段是把取样光栅退化到均匀光栅设置duty1P10×L这样整个光栅段内都有调制等价于均匀布拉格光栅。均匀光栅的峰值反射率有解析解R_theory tanh²(π·dn·L/λ_B)交叉验证代码n_eff 1.452; dn 1e-4; Lambda 534; L 1e6; % 1 mm P 10*L; % 远大于L duty 1; [wl, R] compute_sg_reflectance(n_eff, dn, Lambda, P, duty, L, [1545,1560]); lambda_B 2*n_eff*Lambda; kappa pi*dn/lambda_B; R_theory tanh(kappa*L)^2; R_sim max(R); fprintf(仿真峰值 %f, 解析峰值 %f, 相对误差 %e\n, ... R_sim, R_theory, abs(R_sim-R_theory)/R_theory);如果相对误差小于1e-3基本上可以认为TMM实现正确。如果误差较大优先检查反射系数计算公式以及层间折射率取的是端点值还是平均值。用0.5*(nz(i)nz(i1))取平均误差通常最小。验证通过后再改回取样光栅参数。日常调试时建议把计算好的反射谱保存下来save(sg_reflectance_result.mat, wl, R);下次直接load避免重复计算。如果频繁固定结构参数只扫波长可以在compute_sg_reflectance里用persistent变量缓存折射率数组。这样每次调用函数时如果Lambda、P、duty、L没有变化就跳过重新生成nz的那一步只更新波长循环能让参数优化耗时减少近一半。本文还有配套的精品资源点击获取