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 · Δn（v是条纹可见度，均匀光栅取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段，每段长度Δz=L/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|²。事实上由能量守恒R+T=1。

注意：有的脚本把方向写反，矩阵乘序不同，最后取的矩阵元也不同。你只要记住这个规矩：最终用到的只有M(1,1)、M(2,1)两个元素。

3. 手写一个可运行的matlab脚本：反射谱/透射谱的完整实现

3.1 参数区：波长、光栅结构、折射率调制

我现在给出一份可以直接运行的Matlab脚本，用于计算均匀FBG的反射谱和透射谱。参数先设为1550nm中心波长，n_eff=1.45，光栅长度L=5mm，折射率调制幅度Δn=2e-4，条纹可见度v=1。扫描波长范围取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太大，计算变慢。经验值是每段长度要远小于光栅周期对应的相位变化量。对于均匀光栅，N=200~1000足够。对于啁啾光栅，需要保证每段长度内啁啾造成的相位变化远小于π，通常N取2000以上。

上面代码中每个小段矩阵假设没有啁啾项和切趾项，所以每段的kappa和delta都一样。如果你要计算切趾光栅，只需在循环内把dn替换成dnexp(-((m-0.5)dz - L/2)^2/(2(L0.5)^2))之类的包络。参数sigma_hat的取值要根据光栅平均折射率变化调整，在后续调参部分我会专门举例。

3.3 级联相乘：从起点到终点传递光场

矩阵相乘的顺序很容易出错。正确的方向是从z=0到z=L，也就是从入射端到出射端。如果你的小段数组是按z从0到L排列，那么Mtot应该左乘每一段矩阵，就像上面代码里的Mtot = Mtot * M。如果搞反了，反射率会算成透射率，谱形完全不对。

另外，当gamma为纯虚数时，cosh和sinh实际上转换为cos和sin，双曲函数不会报错，但可能产生数值极大的中间值。如果N过大且L很长，Mtot的每个元素可能接近Inf，导致r计算失败。此时建议改用归一化矩阵传播法，或者用传递矩阵的散射矩阵形式。不过对于常见的毫米级光栅，N=500不会爆。

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，先别急着改参数，用默认参数运行一遍，然后用下面几个方法验证：

• 检查R+T是否等于1（在数值精度内）。
• 将中心波长输出与根据公式λ_B=2n_effΛ计算的值比较，差异应在扫描步长内。
• 用已知解析解对比：均匀弱光栅的峰值反射率约等于 tanh²(κL)。例如Δn=2e-4，λ_B=1550nm，κ=π/λ*Δn≈405.4 m⁻¹，L=5mm，κL≈2.03，tanh²(2.03)≈0.93。如果程序结果接近这个值，说明靠谱。

如果三者都没问题，那这份程序基本可信。如果R+T不是1，先看是不是在画图前做了归一化，再看矩阵乘逆顺序。

4. 调参实验：为什么长度和调制深度会把谱形变成完全不同的样子

4.1 固定调制深度，拉长光栅：峰值反射率迅速逼近1

我固定Δn=1e-4，分别取L=1mm、2mm、5mm、10mm，计算峰值反射率。κ=π/λ*Δn ≈202.7 m⁻¹。κL分别是0.203、0.405、1.014、2.027。根据R_max=tanh²(κL)，对应结果如下：

可以看出：L从1mm到10mm，反射率从几乎透明变成强反射，但带宽也会变窄。原因是光栅长度变长后，参与相长干涉的周期数增多，波长选择性更强。

所以做传感器时，长度和灵敏度要平衡。L太长，反射峰虽高但窄，解调仪不一定跟得上；L太短，反射信号弱。通常传感用FBG长度取3~10mm。

4.2 固定长度，增大Δn：带宽明显展宽，旁瓣抬头

另一个调参方向是折射率调制深度Δn。固定L=5mm，Δ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，说明矩阵算错还是归一化有问题

带外没有反射，T=1才对。如果跑出来T带外是0.9或1.1，先检查扫描波长范围是否远离中心。均匀光栅中，远离布拉格波长时δ很大，γ虚部很大，矩阵元素相位快速振荡。如果N不足，数值上会产生虚反射，带外T就不是1。把N从500提高到2000，这个问题一般会消失。

另一个可能是矩阵乘序错。如果Mtot乘反，相当于把透射波和反射波方向弄反，R/T谱和真实光栅的左右关系会乱。

5.2 N取太小会出来“台阶状”反射谱

如果反射谱出现不规则毛刺或台阶，大概率是N太小。传输矩阵法用分段常数近似连续光栅，段数太小会让谱线呈现周期性伪峰。我实测L=5mm，N=50时，反射谱会像梳子一样；N=500时平滑。所以建议起步N取500，追求精度取2000。注意N太大会让内层循环占用大量内存，但对1000个波长点、N=2000也没压力。

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 验证脚本的黄金标准：R+T=1，以及和解析解对比

最后，无论你用的是FBG.rar里的代码，还是自己写的脚本，发布或分享前一定要做一次自检。最有效的自检是检查能量守恒R+T是否恒等于1（在非吸收条件下）。在均匀光栅情况下，峰值反射率和解析解tanh²(κL)对比；带外反射应趋于0。如果这两项通过，你的程序基本可靠。

我个人的习惯是在代码末尾加一个断言：

assert(max(abs(R + T - 1)) < 1e-10, 'R+T不等于1');

如果矩阵乘用的是标准传输矩阵，物理上R+T应严格等于1，但由于数值误差会略有偏差，1e-10的阈值足够。如果断言失败，优先检查矩阵方向。

FBG.rar这种打包分享的Matlab程序，最大的价值不是那几十行代码，而是让你从运行中接受一次物理建模的训练。我自己第一次跑出反射谱的时候，第一反应是把N改成5000看它会不会更漂亮，结果发现谱形几乎不变，反而意识到N=500已经足够。这种直觉不是看公式能得到的，必须亲手调几次参数、踩几个坑才能真正建立起来。希望这篇文章能帮你少走这点弯路。

本文还有配套的精品资源，点击获取