简介本资源是一份面向光学工程、光纤传感及通信专业初学者与科研人员的MATLAB仿真工具包聚焦光纤布拉格光栅FBG核心特性建模解决反射谱与透射谱可视化理解与参数影响分析的实际需求。压缩包为RAR格式仅含1个关键文件——MATLAB脚本.m体积仅614B轻量简洁适用于快速运行、参数调试与原理验证。该脚本基于布拉格条件λ_Bragg 2nEff·Λ构建理论模型通过矩阵传输法或傅里叶变换实现反射/透射系数计算并直接绘制波长-强度响应曲线支持用户灵活修改光栅周期、有效折射率、调制深度及光栅长度等参数直观观察其对谱宽、峰值、边模抑制比等性能指标的影响。目前已有307人学习下载是理解FBG工作机理、开展课程设计、辅助实验教学及初步器件仿真的实用入门级代码资源。 我拿到过不少FBG的仿真程序其中FBG.rar这种命名十有八九是一份打包好的Matlab工程打开之后要么是.m脚本要么是一堆数据和几张反射谱/透射谱图。很多人解压后第一件事是点运行看能不能直接画出反射峰然后就被一堆矩阵运算搞懵。这篇文章就把FBG的反射谱和透射谱是怎么回事、Matlab程序怎么从零写出来、调参时哪些坑最容易踩一次性讲透。不管你是从网盘下载了FBG.rar想复现还是实验室里需要自己搭一个光栅仿真脚本下面的内容都能直接照着做。1. 拆开FBG.rar之前先想清楚反射谱/透射谱要算的是什么1.1 FBG像一个波长选择器光纤布拉格光栅Fiber Bragg Grating, FBG的本质是在纤芯里写入一段周期性折射率调制。这段调制相当于一个很窄的反射式滤波器波长满足布拉格条件的光会被反射回去其他波长正常通过。布拉格波长由光栅周期和有效折射率决定λ_B 2 n_eff Λ对于1550 nm波段的通信光纤n_eff通常取1.45对应的周期Λ534.48 nm。如果是紫外写入的FBG这个周期就是几百纳米量级。仿真要算的反射谱就是R(λ)这个函数——它在λ_B附近出现峰值透射谱则是T(λ)它和反射谱互补所以带内出现凹陷。这里必须强调一下反射谱和透射谱并不是互不相关的两条曲线。不考虑吸收和散射时R(λ)T(λ)1。我在后面还会拿这个等式当调试工具如果算完两条谱相加不等于1代码里必然有问题。1.2 反射谱与透射谱一张图表的两个视角在通信和传感领域通常先看反射谱因为传感器读取解调仪时接的就是反射光但在器件设计时透射谱也很常用因为它能直观反映插入损耗和阻带深度。大多数Matlab程序会把两个谱都画在同一张图里左边轴反射率右边轴透射率这样能一眼看出带宽和旁瓣。用传输矩阵法仿真时R和T实际上来自同一个2×2矩阵并不是两次独立计算。这也是为什么很多FBG.rar里的脚本结构看起来一模一样输入参数、切段、循环矩阵、画图。理解了这一点你拿到任何人的脚本都能快速看懂主逻辑。1.3 什么时候抄现成脚本什么时候自己从头写如果你只是需要一张均匀FBG的反射谱示意图直接跑现成脚本就够了。但如果你要做啁啾光栅、相移光栅、切趾光栅或者要把温度和应变影响同时加进模型那么现成代码往往需要大改。此时建议自己从传输矩阵法写一个通用的核心函数我下面给出的脚本就是按这个思路设计的。这个函数很短几十行但足以覆盖均匀、切趾和啁啾三类常见FBG。2. 传输矩阵法里的“为什么”反射谱和透射谱背后的耦合模理论2.1 从Maxwell到耦合模中间砍掉了哪些细节严格求解光在周期波导里的传输需要解矢量亥姆霍兹方程数值上可以用FDTD、有限元等但速度慢。耦合模理论利用“正向模和反向模之间能量交换”这个物理图像把问题化简为一组一阶常微分方程。光栅的折射率扰动一般只影响模式的有效折射率不改变模场分布所以单模光纤里的FBG用一个近似标量模型就够。对于均匀周期光栅耦合模方程有解析解对于非均匀光栅如切趾或啁啾需要把整段光栅切成很多小段每段当均匀光栅处理再用传输矩阵级联。这就是传输矩阵法Transfer Matrix Method, TMM。它兼顾精度和速度是当前FBG仿真最主流的工具。2.2 σ、κ、δ到底都是什么含义在Matlab代码里你会看到几个变量delta失谐量、kappa交流耦合系数、sigma_hat直流自耦合系数。这三个参量搞不明白程序就只能黑盒跑。用式子表达失谐量 δ 2π n_eff (1/λ - 1/λ_B)交流耦合系数 κ π/λ · v · Δnv是条纹可见度均匀光栅取1直流自耦合系数 σ_hat 2π/λ · Δn_effΔn_eff是有效折射率平均变化若折射率调制均值为0通常忽略这里最容易混的是Δn到底是峰值折射率调制幅度还是平均变化在光栅写入过程中通常用 Δn 表示折射率调制的幅度即 n(z) n_eff Δn·cos(2πz/Λ) 的振幅。此时平均变化可以认为是0所以 σ_hat ≈ 0除非有额外的平均折射率漂移。但如果光栅不是理想正弦调制比如折射率变化是 n_eff Δn/2 Δn/2 cos(...)那么直流项就是 Δn/2。很多现成脚本里把σ_hat设成0就是这个原因。还有一个关键系数γ sqrt(kappa^2 - sigma_hat^2)它决定矩阵里三角函数的自变量是真三角函数还是双曲函数。当kappa sigma_hat时γ是实数反射谱在带内振荡否则是一个单调衰减。别乱开根号要写成复数的sqrt。2.3 传输矩阵为什么能一段一段级联将长度为L的光栅分成N段每段长度ΔzL/N。每一小段内的光栅看成均匀周期光栅用2×2矩阵M_i描述输入输出场关系。整段光栅的传输矩阵就是所有小段矩阵的乘积M M_N · ... · M_1。这个思想类似电路里的二端口网络级联前一段的输出就是后一段的输入方向必须固定。每段矩阵元在Matlab里的常见写法是gamma sqrt(kappa^2 - sigma_hat^2); A cosh(gamma * dz) - 1i * sigma_hat/gamma * sinh(gamma * dz); B -1i * kappa/gamma * sinh(gamma * dz); C 1i * kappa/gamma * sinh(gamma * dz); D cosh(gamma * dz) 1i * sigma_hat/gamma * sinh(gamma * dz);如果gamma是虚数cosh(gammadz)实际会变成cos(|gamma|dz)sinh(gammadz)变成isin(...)/? 所以直接用复数sqrt不会出问题。2.4 边界条件反射率和透射率怎么从矩阵元里取出来光从光栅左侧入射假设反射波只向后传最右端没有反向波输入。设左侧前向场为E_f(0)反向场为E_b(0)右侧前向场为E_f(L)反向场为E_b(L)0。传输矩阵写成[E_f(0); E_b(0)] M * [E_f(L); 0]这里M用小段矩阵倒序相乘得到。于是反射系数 r E_b(0)/E_f(0) M[2,1] / M[1,1]透射系数 t E_f(L)/E_f(0) 1 / M[1,1]。反射率 R|r|²透射率 T|t|²。事实上由能量守恒RT1。注意有的脚本把方向写反矩阵乘序不同最后取的矩阵元也不同。你只要记住这个规矩最终用到的只有M(1,1)、M(2,1)两个元素。3. 手写一个可运行的matlab脚本反射谱/透射谱的完整实现3.1 参数区波长、光栅结构、折射率调制我现在给出一份可以直接运行的Matlab脚本用于计算均匀FBG的反射谱和透射谱。参数先设为1550nm中心波长n_eff1.45光栅长度L5mm折射率调制幅度Δn2e-4条纹可见度v1。扫描波长范围取1549~1551nm点数1000。注意单位统一用米。lambda_B 1550e-9; n_eff 1.45; Lambda lambda_B / (2 * n_eff); L 5e-3; dn 2e-4; v 1.0; N 500; dz L / N; lambda linspace(1549e-9, 1551e-9, 1000); R zeros(size(lambda)); T zeros(size(lambda)); for idx 1:length(lambda) lam lambda(idx); delta 2*pi*n_eff*(1/lam - 1/lambda_B); kappa pi/lam * v * dn; sigma_hat 0; gamma sqrt(kappa^2 - sigma_hat^2); Mtot eye(2); for m 1:N A cosh(gamma*dz) - 1i*sigma_hat/gamma*sinh(gamma*dz); B -1i*kappa/gamma*sinh(gamma*dz); C 1i*kappa/gamma*sinh(gamma*dz); D cosh(gamma*dz) 1i*sigma_hat/gamma*sinh(gamma*dz); M [A B; C D]; Mtot Mtot * M; end r Mtot(2,1) / Mtot(1,1); t 1 / Mtot(1,1); R(idx) abs(r)^2; T(idx) abs(t)^2; end figure; plot(lambda*1e9, R, b-, LineWidth, 1.5); hold on; plot(lambda*1e9, T, r--, LineWidth, 1.5); xlabel(波长 (nm)); ylabel(反射率 / 透射率); legend(R, T); grid on;这段代码最外层的循环是波长扫描内层是N个小段矩阵级联。N取500时对1000个波长点计算单次大概几十毫秒到一二百毫秒完全够用。3.2 切分N段每段的矩阵怎么算N的选择很讲究。N太小反射谱会出现振荡和台阶N太大计算变慢。经验值是每段长度要远小于光栅周期对应的相位变化量。对于均匀光栅N200~1000足够。对于啁啾光栅需要保证每段长度内啁啾造成的相位变化远小于π通常N取2000以上。上面代码中每个小段矩阵假设没有啁啾项和切趾项所以每段的kappa和delta都一样。如果你要计算切趾光栅只需在循环内把dn替换成dnexp(-((m-0.5)dz - L/2)^2/(2(L0.5)^2))之类的包络。参数sigma_hat的取值要根据光栅平均折射率变化调整在后续调参部分我会专门举例。3.3 级联相乘从起点到终点传递光场矩阵相乘的顺序很容易出错。正确的方向是从z0到zL也就是从入射端到出射端。如果你的小段数组是按z从0到L排列那么Mtot应该左乘每一段矩阵就像上面代码里的Mtot Mtot * M。如果搞反了反射率会算成透射率谱形完全不对。另外当gamma为纯虚数时cosh和sinh实际上转换为cos和sin双曲函数不会报错但可能产生数值极大的中间值。如果N过大且L很长Mtot的每个元素可能接近Inf导致r计算失败。此时建议改用归一化矩阵传播法或者用传递矩阵的散射矩阵形式。不过对于常见的毫米级光栅N500不会爆。3.4 提取反射率R和透射率T并画图上一节给出了r和t的公式。务必用abs(r)^2而不是直接取Mtot(2,1)^2因为r是复数平方后相位与1/M(1,1)有关。很多新手的错误是直接R abs(Mtot(2,1))^2这在某些边界条件下恰好数值正确但原理不对因为归一化分母被忽略了。画图时建议在同一个figure中双y轴绘制或者把反射谱用线性坐标、透射谱用dB坐标分开看图。如果你是设计FBG传感器关心的是反射峰移动量用线性坐标就够如果是做滤波器件透射深度用dB更直观。3.5 如果只想用现成的FBG.rar怎么快速验证它靠谱拿到别人的FBG.rar先别急着改参数用默认参数运行一遍然后用下面几个方法验证检查RT是否等于1在数值精度内。将中心波长输出与根据公式λ_B2n_effΛ计算的值比较差异应在扫描步长内。用已知解析解对比均匀弱光栅的峰值反射率约等于 tanh²(κL)。例如Δn2e-4λ_B1550nmκπ/λ*Δn≈405.4 m⁻¹L5mmκL≈2.03tanh²(2.03)≈0.93。如果程序结果接近这个值说明靠谱。如果三者都没问题那这份程序基本可信。如果RT不是1先看是不是在画图前做了归一化再看矩阵乘逆顺序。4. 调参实验为什么长度和调制深度会把谱形变成完全不同的样子4.1 固定调制深度拉长光栅峰值反射率迅速逼近1我固定Δn1e-4分别取L1mm、2mm、5mm、10mm计算峰值反射率。κπ/λ*Δn ≈202.7 m⁻¹。κL分别是0.203、0.405、1.014、2.027。根据R_maxtanh²(κL)对应结果如下光栅长度 LκL峰值反射率 R_max谱线特征1 mm0.200.04弱反射近透明2 mm0.410.15反射峰出现带宽较宽5 mm1.010.63峰明显带宽变窄10 mm2.030.93强反射旁瓣明显可以看出L从1mm到10mm反射率从几乎透明变成强反射但带宽也会变窄。原因是光栅长度变长后参与相长干涉的周期数增多波长选择性更强。所以做传感器时长度和灵敏度要平衡。L太长反射峰虽高但窄解调仪不一定跟得上L太短反射信号弱。通常传感用FBG长度取3~10mm。4.2 固定长度增大Δn带宽明显展宽旁瓣抬头另一个调参方向是折射率调制深度Δn。固定L5mmΔn从1e-4增到1e-3。κL从1.01增到10.1。峰值反射率接近1但带宽从几十pm变成几百pm。这是因为强光栅的禁带宽度由κ决定κ越大反射带宽越宽。同时均匀光栅的旁瓣会非常明显。旁瓣峰值出现在带外约十几个pm处甚至第一旁瓣反射率可达几个百分点在传感应用中可能被误认为第二个反射峰。这就是为什么实用FBG大多要做切趾。4.3 切趾光栅用高斯包络换掉矩形包络旁瓣能被压下去切趾的思想很简单把均匀分布的折射率调制幅度改成沿长度方向平滑变化一般是中间大、两端小。高斯切趾的包络函数是exp(-(z-L/2)²/(2σ_z²))其中σ_z是高斯宽度。代码里可以这样改动for m 1:N z_m (m - 0.5) * dz; apod exp(-(z_m - L/2)^2 / (2 * (L/4)^2)); kappa_m pi/lam * v * dn * apod; sigma_hat 0; gamma sqrt(kappa_m^2 - sigma_hat^2); A cosh(gamma*dz) - 1i*sigma_hat/gamma*sinh(gamma*dz); B -1i*kappa_m/gamma*sinh(gamma*dz); C 1i*kappa_m/gamma*sinh(gamma*dz); D cosh(gamma*dz) 1i*sigma_hat/gamma*sinh(gamma*dz); M [A B; C D]; Mtot Mtot * M; end高斯宽度选L/4时两端调制幅度降到约exp(-2)0.135足够把旁瓣压到-40dB以下。代价是主峰带宽略有展宽峰值反射率略降。对于传感应用这种旁瓣抑制非常重要。4.4 用同一段代码对比均匀FBG vs 切趾FBG我习惯把均匀和切趾放在同一个脚本里跑输出两张谱图做对比。均匀光栅在带外有一串旁瓣第一旁瓣约-20dB切趾后旁瓣几乎贴着零轴。如果你拿到的FBG.rar里只有均匀光栅代码可以按上面的apod变量自己扩展基本不用改其他部分。5. 跑仿真时躲不开的几个坑边界、网格和程序包问题5.1 透射谱带外不等于1说明矩阵算错还是归一化有问题带外没有反射T1才对。如果跑出来T带外是0.9或1.1先检查扫描波长范围是否远离中心。均匀光栅中远离布拉格波长时δ很大γ虚部很大矩阵元素相位快速振荡。如果N不足数值上会产生虚反射带外T就不是1。把N从500提高到2000这个问题一般会消失。另一个可能是矩阵乘序错。如果Mtot乘反相当于把透射波和反射波方向弄反R/T谱和真实光栅的左右关系会乱。5.2 N取太小会出来“台阶状”反射谱如果反射谱出现不规则毛刺或台阶大概率是N太小。传输矩阵法用分段常数近似连续光栅段数太小会让谱线呈现周期性伪峰。我实测L5mmN50时反射谱会像梳子一样N500时平滑。所以建议起步N取500追求精度取2000。注意N太大会让内层循环占用大量内存但对1000个波长点、N2000也没压力。5.3 RAR包里的程序最常见的三个运行报错拿到FBG.rar解压后运行程序时常见几个报错未定义变量脚本里用了中文变量名或缺少参数初始化。很多分享版程序写得不规范需要你补上行clear; clc; close all;并检查路径。矩阵维度不匹配通常是因为lambda是行向量而计算中把它当标量用或.*写成了*。图形窗口闪退老程序可能用plot(lambda, R)而lambda是nm单位画出来的横轴数值是1.55e-6看起来像0。改成plot(lambda*1e9, R)即可。5.4 验证脚本的黄金标准RT1以及和解析解对比最后无论你用的是FBG.rar里的代码还是自己写的脚本发布或分享前一定要做一次自检。最有效的自检是检查能量守恒RT是否恒等于1在非吸收条件下。在均匀光栅情况下峰值反射率和解析解tanh²(κL)对比带外反射应趋于0。如果这两项通过你的程序基本可靠。我个人的习惯是在代码末尾加一个断言assert(max(abs(R T - 1)) 1e-10, RT不等于1);如果矩阵乘用的是标准传输矩阵物理上RT应严格等于1但由于数值误差会略有偏差1e-10的阈值足够。如果断言失败优先检查矩阵方向。FBG.rar这种打包分享的Matlab程序最大的价值不是那几十行代码而是让你从运行中接受一次物理建模的训练。我自己第一次跑出反射谱的时候第一反应是把N改成5000看它会不会更漂亮结果发现谱形几乎不变反而意识到N500已经足够。这种直觉不是看公式能得到的必须亲手调几次参数、踩几个坑才能真正建立起来。希望这篇文章能帮你少走这点弯路。本文还有配套的精品资源点击获取
FBG反射谱与透射谱的Matlab仿真：从传输矩阵到调参避坑
简介本资源是一份面向光学工程、光纤传感及通信专业初学者与科研人员的MATLAB仿真工具包聚焦光纤布拉格光栅FBG核心特性建模解决反射谱与透射谱可视化理解与参数影响分析的实际需求。压缩包为RAR格式仅含1个关键文件——MATLAB脚本.m体积仅614B轻量简洁适用于快速运行、参数调试与原理验证。该脚本基于布拉格条件λ_Bragg 2nEff·Λ构建理论模型通过矩阵传输法或傅里叶变换实现反射/透射系数计算并直接绘制波长-强度响应曲线支持用户灵活修改光栅周期、有效折射率、调制深度及光栅长度等参数直观观察其对谱宽、峰值、边模抑制比等性能指标的影响。目前已有307人学习下载是理解FBG工作机理、开展课程设计、辅助实验教学及初步器件仿真的实用入门级代码资源。 我拿到过不少FBG的仿真程序其中FBG.rar这种命名十有八九是一份打包好的Matlab工程打开之后要么是.m脚本要么是一堆数据和几张反射谱/透射谱图。很多人解压后第一件事是点运行看能不能直接画出反射峰然后就被一堆矩阵运算搞懵。这篇文章就把FBG的反射谱和透射谱是怎么回事、Matlab程序怎么从零写出来、调参时哪些坑最容易踩一次性讲透。不管你是从网盘下载了FBG.rar想复现还是实验室里需要自己搭一个光栅仿真脚本下面的内容都能直接照着做。1. 拆开FBG.rar之前先想清楚反射谱/透射谱要算的是什么1.1 FBG像一个波长选择器光纤布拉格光栅Fiber Bragg Grating, FBG的本质是在纤芯里写入一段周期性折射率调制。这段调制相当于一个很窄的反射式滤波器波长满足布拉格条件的光会被反射回去其他波长正常通过。布拉格波长由光栅周期和有效折射率决定λ_B 2 n_eff Λ对于1550 nm波段的通信光纤n_eff通常取1.45对应的周期Λ534.48 nm。如果是紫外写入的FBG这个周期就是几百纳米量级。仿真要算的反射谱就是R(λ)这个函数——它在λ_B附近出现峰值透射谱则是T(λ)它和反射谱互补所以带内出现凹陷。这里必须强调一下反射谱和透射谱并不是互不相关的两条曲线。不考虑吸收和散射时R(λ)T(λ)1。我在后面还会拿这个等式当调试工具如果算完两条谱相加不等于1代码里必然有问题。1.2 反射谱与透射谱一张图表的两个视角在通信和传感领域通常先看反射谱因为传感器读取解调仪时接的就是反射光但在器件设计时透射谱也很常用因为它能直观反映插入损耗和阻带深度。大多数Matlab程序会把两个谱都画在同一张图里左边轴反射率右边轴透射率这样能一眼看出带宽和旁瓣。用传输矩阵法仿真时R和T实际上来自同一个2×2矩阵并不是两次独立计算。这也是为什么很多FBG.rar里的脚本结构看起来一模一样输入参数、切段、循环矩阵、画图。理解了这一点你拿到任何人的脚本都能快速看懂主逻辑。1.3 什么时候抄现成脚本什么时候自己从头写如果你只是需要一张均匀FBG的反射谱示意图直接跑现成脚本就够了。但如果你要做啁啾光栅、相移光栅、切趾光栅或者要把温度和应变影响同时加进模型那么现成代码往往需要大改。此时建议自己从传输矩阵法写一个通用的核心函数我下面给出的脚本就是按这个思路设计的。这个函数很短几十行但足以覆盖均匀、切趾和啁啾三类常见FBG。2. 传输矩阵法里的“为什么”反射谱和透射谱背后的耦合模理论2.1 从Maxwell到耦合模中间砍掉了哪些细节严格求解光在周期波导里的传输需要解矢量亥姆霍兹方程数值上可以用FDTD、有限元等但速度慢。耦合模理论利用“正向模和反向模之间能量交换”这个物理图像把问题化简为一组一阶常微分方程。光栅的折射率扰动一般只影响模式的有效折射率不改变模场分布所以单模光纤里的FBG用一个近似标量模型就够。对于均匀周期光栅耦合模方程有解析解对于非均匀光栅如切趾或啁啾需要把整段光栅切成很多小段每段当均匀光栅处理再用传输矩阵级联。这就是传输矩阵法Transfer Matrix Method, TMM。它兼顾精度和速度是当前FBG仿真最主流的工具。2.2 σ、κ、δ到底都是什么含义在Matlab代码里你会看到几个变量delta失谐量、kappa交流耦合系数、sigma_hat直流自耦合系数。这三个参量搞不明白程序就只能黑盒跑。用式子表达失谐量 δ 2π n_eff (1/λ - 1/λ_B)交流耦合系数 κ π/λ · v · Δnv是条纹可见度均匀光栅取1直流自耦合系数 σ_hat 2π/λ · Δn_effΔn_eff是有效折射率平均变化若折射率调制均值为0通常忽略这里最容易混的是Δn到底是峰值折射率调制幅度还是平均变化在光栅写入过程中通常用 Δn 表示折射率调制的幅度即 n(z) n_eff Δn·cos(2πz/Λ) 的振幅。此时平均变化可以认为是0所以 σ_hat ≈ 0除非有额外的平均折射率漂移。但如果光栅不是理想正弦调制比如折射率变化是 n_eff Δn/2 Δn/2 cos(...)那么直流项就是 Δn/2。很多现成脚本里把σ_hat设成0就是这个原因。还有一个关键系数γ sqrt(kappa^2 - sigma_hat^2)它决定矩阵里三角函数的自变量是真三角函数还是双曲函数。当kappa sigma_hat时γ是实数反射谱在带内振荡否则是一个单调衰减。别乱开根号要写成复数的sqrt。2.3 传输矩阵为什么能一段一段级联将长度为L的光栅分成N段每段长度ΔzL/N。每一小段内的光栅看成均匀周期光栅用2×2矩阵M_i描述输入输出场关系。整段光栅的传输矩阵就是所有小段矩阵的乘积M M_N · ... · M_1。这个思想类似电路里的二端口网络级联前一段的输出就是后一段的输入方向必须固定。每段矩阵元在Matlab里的常见写法是gamma sqrt(kappa^2 - sigma_hat^2); A cosh(gamma * dz) - 1i * sigma_hat/gamma * sinh(gamma * dz); B -1i * kappa/gamma * sinh(gamma * dz); C 1i * kappa/gamma * sinh(gamma * dz); D cosh(gamma * dz) 1i * sigma_hat/gamma * sinh(gamma * dz);如果gamma是虚数cosh(gammadz)实际会变成cos(|gamma|dz)sinh(gammadz)变成isin(...)/? 所以直接用复数sqrt不会出问题。2.4 边界条件反射率和透射率怎么从矩阵元里取出来光从光栅左侧入射假设反射波只向后传最右端没有反向波输入。设左侧前向场为E_f(0)反向场为E_b(0)右侧前向场为E_f(L)反向场为E_b(L)0。传输矩阵写成[E_f(0); E_b(0)] M * [E_f(L); 0]这里M用小段矩阵倒序相乘得到。于是反射系数 r E_b(0)/E_f(0) M[2,1] / M[1,1]透射系数 t E_f(L)/E_f(0) 1 / M[1,1]。反射率 R|r|²透射率 T|t|²。事实上由能量守恒RT1。注意有的脚本把方向写反矩阵乘序不同最后取的矩阵元也不同。你只要记住这个规矩最终用到的只有M(1,1)、M(2,1)两个元素。3. 手写一个可运行的matlab脚本反射谱/透射谱的完整实现3.1 参数区波长、光栅结构、折射率调制我现在给出一份可以直接运行的Matlab脚本用于计算均匀FBG的反射谱和透射谱。参数先设为1550nm中心波长n_eff1.45光栅长度L5mm折射率调制幅度Δn2e-4条纹可见度v1。扫描波长范围取1549~1551nm点数1000。注意单位统一用米。lambda_B 1550e-9; n_eff 1.45; Lambda lambda_B / (2 * n_eff); L 5e-3; dn 2e-4; v 1.0; N 500; dz L / N; lambda linspace(1549e-9, 1551e-9, 1000); R zeros(size(lambda)); T zeros(size(lambda)); for idx 1:length(lambda) lam lambda(idx); delta 2*pi*n_eff*(1/lam - 1/lambda_B); kappa pi/lam * v * dn; sigma_hat 0; gamma sqrt(kappa^2 - sigma_hat^2); Mtot eye(2); for m 1:N A cosh(gamma*dz) - 1i*sigma_hat/gamma*sinh(gamma*dz); B -1i*kappa/gamma*sinh(gamma*dz); C 1i*kappa/gamma*sinh(gamma*dz); D cosh(gamma*dz) 1i*sigma_hat/gamma*sinh(gamma*dz); M [A B; C D]; Mtot Mtot * M; end r Mtot(2,1) / Mtot(1,1); t 1 / Mtot(1,1); R(idx) abs(r)^2; T(idx) abs(t)^2; end figure; plot(lambda*1e9, R, b-, LineWidth, 1.5); hold on; plot(lambda*1e9, T, r--, LineWidth, 1.5); xlabel(波长 (nm)); ylabel(反射率 / 透射率); legend(R, T); grid on;这段代码最外层的循环是波长扫描内层是N个小段矩阵级联。N取500时对1000个波长点计算单次大概几十毫秒到一二百毫秒完全够用。3.2 切分N段每段的矩阵怎么算N的选择很讲究。N太小反射谱会出现振荡和台阶N太大计算变慢。经验值是每段长度要远小于光栅周期对应的相位变化量。对于均匀光栅N200~1000足够。对于啁啾光栅需要保证每段长度内啁啾造成的相位变化远小于π通常N取2000以上。上面代码中每个小段矩阵假设没有啁啾项和切趾项所以每段的kappa和delta都一样。如果你要计算切趾光栅只需在循环内把dn替换成dnexp(-((m-0.5)dz - L/2)^2/(2(L0.5)^2))之类的包络。参数sigma_hat的取值要根据光栅平均折射率变化调整在后续调参部分我会专门举例。3.3 级联相乘从起点到终点传递光场矩阵相乘的顺序很容易出错。正确的方向是从z0到zL也就是从入射端到出射端。如果你的小段数组是按z从0到L排列那么Mtot应该左乘每一段矩阵就像上面代码里的Mtot Mtot * M。如果搞反了反射率会算成透射率谱形完全不对。另外当gamma为纯虚数时cosh和sinh实际上转换为cos和sin双曲函数不会报错但可能产生数值极大的中间值。如果N过大且L很长Mtot的每个元素可能接近Inf导致r计算失败。此时建议改用归一化矩阵传播法或者用传递矩阵的散射矩阵形式。不过对于常见的毫米级光栅N500不会爆。3.4 提取反射率R和透射率T并画图上一节给出了r和t的公式。务必用abs(r)^2而不是直接取Mtot(2,1)^2因为r是复数平方后相位与1/M(1,1)有关。很多新手的错误是直接R abs(Mtot(2,1))^2这在某些边界条件下恰好数值正确但原理不对因为归一化分母被忽略了。画图时建议在同一个figure中双y轴绘制或者把反射谱用线性坐标、透射谱用dB坐标分开看图。如果你是设计FBG传感器关心的是反射峰移动量用线性坐标就够如果是做滤波器件透射深度用dB更直观。3.5 如果只想用现成的FBG.rar怎么快速验证它靠谱拿到别人的FBG.rar先别急着改参数用默认参数运行一遍然后用下面几个方法验证检查RT是否等于1在数值精度内。将中心波长输出与根据公式λ_B2n_effΛ计算的值比较差异应在扫描步长内。用已知解析解对比均匀弱光栅的峰值反射率约等于 tanh²(κL)。例如Δn2e-4λ_B1550nmκπ/λ*Δn≈405.4 m⁻¹L5mmκL≈2.03tanh²(2.03)≈0.93。如果程序结果接近这个值说明靠谱。如果三者都没问题那这份程序基本可信。如果RT不是1先看是不是在画图前做了归一化再看矩阵乘逆顺序。4. 调参实验为什么长度和调制深度会把谱形变成完全不同的样子4.1 固定调制深度拉长光栅峰值反射率迅速逼近1我固定Δn1e-4分别取L1mm、2mm、5mm、10mm计算峰值反射率。κπ/λ*Δn ≈202.7 m⁻¹。κL分别是0.203、0.405、1.014、2.027。根据R_maxtanh²(κL)对应结果如下光栅长度 LκL峰值反射率 R_max谱线特征1 mm0.200.04弱反射近透明2 mm0.410.15反射峰出现带宽较宽5 mm1.010.63峰明显带宽变窄10 mm2.030.93强反射旁瓣明显可以看出L从1mm到10mm反射率从几乎透明变成强反射但带宽也会变窄。原因是光栅长度变长后参与相长干涉的周期数增多波长选择性更强。所以做传感器时长度和灵敏度要平衡。L太长反射峰虽高但窄解调仪不一定跟得上L太短反射信号弱。通常传感用FBG长度取3~10mm。4.2 固定长度增大Δn带宽明显展宽旁瓣抬头另一个调参方向是折射率调制深度Δn。固定L5mmΔn从1e-4增到1e-3。κL从1.01增到10.1。峰值反射率接近1但带宽从几十pm变成几百pm。这是因为强光栅的禁带宽度由κ决定κ越大反射带宽越宽。同时均匀光栅的旁瓣会非常明显。旁瓣峰值出现在带外约十几个pm处甚至第一旁瓣反射率可达几个百分点在传感应用中可能被误认为第二个反射峰。这就是为什么实用FBG大多要做切趾。4.3 切趾光栅用高斯包络换掉矩形包络旁瓣能被压下去切趾的思想很简单把均匀分布的折射率调制幅度改成沿长度方向平滑变化一般是中间大、两端小。高斯切趾的包络函数是exp(-(z-L/2)²/(2σ_z²))其中σ_z是高斯宽度。代码里可以这样改动for m 1:N z_m (m - 0.5) * dz; apod exp(-(z_m - L/2)^2 / (2 * (L/4)^2)); kappa_m pi/lam * v * dn * apod; sigma_hat 0; gamma sqrt(kappa_m^2 - sigma_hat^2); A cosh(gamma*dz) - 1i*sigma_hat/gamma*sinh(gamma*dz); B -1i*kappa_m/gamma*sinh(gamma*dz); C 1i*kappa_m/gamma*sinh(gamma*dz); D cosh(gamma*dz) 1i*sigma_hat/gamma*sinh(gamma*dz); M [A B; C D]; Mtot Mtot * M; end高斯宽度选L/4时两端调制幅度降到约exp(-2)0.135足够把旁瓣压到-40dB以下。代价是主峰带宽略有展宽峰值反射率略降。对于传感应用这种旁瓣抑制非常重要。4.4 用同一段代码对比均匀FBG vs 切趾FBG我习惯把均匀和切趾放在同一个脚本里跑输出两张谱图做对比。均匀光栅在带外有一串旁瓣第一旁瓣约-20dB切趾后旁瓣几乎贴着零轴。如果你拿到的FBG.rar里只有均匀光栅代码可以按上面的apod变量自己扩展基本不用改其他部分。5. 跑仿真时躲不开的几个坑边界、网格和程序包问题5.1 透射谱带外不等于1说明矩阵算错还是归一化有问题带外没有反射T1才对。如果跑出来T带外是0.9或1.1先检查扫描波长范围是否远离中心。均匀光栅中远离布拉格波长时δ很大γ虚部很大矩阵元素相位快速振荡。如果N不足数值上会产生虚反射带外T就不是1。把N从500提高到2000这个问题一般会消失。另一个可能是矩阵乘序错。如果Mtot乘反相当于把透射波和反射波方向弄反R/T谱和真实光栅的左右关系会乱。5.2 N取太小会出来“台阶状”反射谱如果反射谱出现不规则毛刺或台阶大概率是N太小。传输矩阵法用分段常数近似连续光栅段数太小会让谱线呈现周期性伪峰。我实测L5mmN50时反射谱会像梳子一样N500时平滑。所以建议起步N取500追求精度取2000。注意N太大会让内层循环占用大量内存但对1000个波长点、N2000也没压力。5.3 RAR包里的程序最常见的三个运行报错拿到FBG.rar解压后运行程序时常见几个报错未定义变量脚本里用了中文变量名或缺少参数初始化。很多分享版程序写得不规范需要你补上行clear; clc; close all;并检查路径。矩阵维度不匹配通常是因为lambda是行向量而计算中把它当标量用或.*写成了*。图形窗口闪退老程序可能用plot(lambda, R)而lambda是nm单位画出来的横轴数值是1.55e-6看起来像0。改成plot(lambda*1e9, R)即可。5.4 验证脚本的黄金标准RT1以及和解析解对比最后无论你用的是FBG.rar里的代码还是自己写的脚本发布或分享前一定要做一次自检。最有效的自检是检查能量守恒RT是否恒等于1在非吸收条件下。在均匀光栅情况下峰值反射率和解析解tanh²(κL)对比带外反射应趋于0。如果这两项通过你的程序基本可靠。我个人的习惯是在代码末尾加一个断言assert(max(abs(R T - 1)) 1e-10, RT不等于1);如果矩阵乘用的是标准传输矩阵物理上RT应严格等于1但由于数值误差会略有偏差1e-10的阈值足够。如果断言失败优先检查矩阵方向。FBG.rar这种打包分享的Matlab程序最大的价值不是那几十行代码而是让你从运行中接受一次物理建模的训练。我自己第一次跑出反射谱的时候第一反应是把N改成5000看它会不会更漂亮结果发现谱形几乎不变反而意识到N500已经足够。这种直觉不是看公式能得到的必须亲手调几次参数、踩几个坑才能真正建立起来。希望这篇文章能帮你少走这点弯路。本文还有配套的精品资源点击获取
