啁啾光纤光栅MATLAB仿真：反射谱与群时延计算全解析

简介：本资源是一份面向光学工程、光通信及信号处理方向初学者与实践者的MATLAB仿真工具包，聚焦啁啾光纤光栅（CFBG）核心特性建模，解决反射谱展宽机制与时延响应分析等关键学习难点。压缩包仅含1个MATLAB脚本文件（.m），大小仅1KB，代码精炼，完整实现啁啾率、光栅长度、折射率等参数设置，基于布拉格反射原理与傅里叶变换计算反射谱，并同步输出群时延曲线，支持直观理解脉冲整形与色散补偿原理。已有364人学习下载，适用于高校课程设计、光纤传感实验预研及光器件仿真入门。用户可直接运行脚本观察不同啁啾参数对反射带宽与时延斜率的影响，获取可复用的数值建模框架与可视化结果，无需额外依赖库，是掌握CFBG物理本质与MATLAB光学仿真实践的高效起点。 光栅时延怎么用MATLAB算？啁啾光纤光栅的反射谱和时延又是怎么对应起来的？这是很多刚接触光纤布拉格光栅仿真的朋友都会卡住的问题。手上刚好有一套在MATLAB里做啁啾光栅（Chirped Fiber Bragg Grating, CFBG）反射谱和群时延计算的脚本，跑过不少参数，也踩过不少坑，今天把这套东西完整拆开聊聊。不管你是课程设计要做“matlab啁啾光栅的反射谱、时延”，还是项目里需要评估色散补偿能力，这篇文章都可以直接拿来当参考。

先说结论：这个代码项目解决的是“已知光栅参数，如何得到反射谱和时延曲线”这件事。核心算法是传输矩阵法（Transfer Matrix Method, TMM），相比直接解耦合模方程，它的实现简单、参数直观，适合处理啁啾、切趾、非均匀光栅这一类实际问题。我会把原理、代码、结果解读、排错经验一条条过一遍，尽量让你看完就能自己跑出结果。

1. 项目整体设计与思路拆解

1.1 啁啾光栅的核心物理：为什么反射谱和时延相关

普通光纤布拉格光栅（FBG）的周期是均匀的，只有满足布拉格条件的波长才会被反射，反射谱就是一条窄峰。啁啾光纤光栅的光栅周期沿长度方向线性变化，不同位置对应的布拉格波长不同，于是不同波长的光被反射的位置也不同。

这一点非常关键，因为它同时决定了两件事：

• 反射谱的带宽：周期变化覆盖的波长范围，决定了反射谱宽度。周期变化范围越宽，反射谱越宽。
• 群时延特性：不同波长的反射点在光栅内部深度不同，往返传播距离不同，于是不同波长之间存在时间延迟差。这就是色散补偿的原理。

一台典型的啁啾光栅，长波段和短波段反射点的位置可能相差几厘米到十几厘米，光在光纤中往返这几厘米，就产生了纳秒量级的群时延差。这个时延差除以波长带宽，就是光栅能提供的色散量，单位通常是 ps/nm。

所以，仿真啁啾光栅不能只看反射谱，时延曲线才是它的“灵魂”。设计色散补偿器时，我们关心的是时延曲线是否线性、纹波大不大。时延纹波大会导致脉冲信号畸变，这个问题后面会详细讲。

1.2 传输矩阵法为什么适合这个场景

计算啁啾光栅响应，常见的方法有几种：耦合模方程直接数值积分、Rouard方法（逐层反射叠加）、传输矩阵法（TMM），还有商业软件如OptiGrating。实际写代码验证，我最推荐传输矩阵法。

它把长度为L的光栅切成M段，每一小段视为一个均匀FBG，用一个2x2矩阵描述该段的传输特性，然后把所有矩阵按顺序相乘，最终得到整个光栅的响应。这个过程特别像一个“管道串联模型”：每一小段光栅是一个管道元件，整根光栅就是一堆元件串起来。

为什么选它？

• 实现思路直白，每段只算一个矩阵乘法，没有复杂的微分方程数值求解。
• 天然支持啁啾：每段的布拉格波长不一样，直接按位置设置局部布拉格波长即可。
• 天然支持切趾（apodization）：只需要在每段的耦合系数上乘一个窗函数值。
• 计算量可控：M取几百到几千，波长扫描几千个点，MATLAB运行时间在几秒到几十秒，完全可接受。

相比之下，直接积分耦合模方程实现麻烦，Rouard方法在段数多时效率低。TMM是教学和工程验证之间最平衡的方案。

1.3 仿真的输入参数怎么定

仿真之前，先想清楚要这几个参数：

这里有个容易踩坑的点：啁啾量单位。常见写法是 nm/cm，表示沿光栅每前进1厘米，布拉格波长变化多少纳米。计算带宽的公式很简单：

总带宽 ≈ 啁啾量 × 光栅长度

比如啁啾量 0.5 nm/cm、长度 10 cm，总带宽约 5 nm。如果用10 cm光栅补偿5 nm带宽的信号，做出来的延迟差大约是：

Δτ = 2 × n_eff × L / c = 2 × 1.45 × 0.1 / 3e8 ≈ 0.97 ns

这对应约 193 ps/nm 的色散量，大约能补偿 11 km 标准单模光纤（17 ps/nm/km）。这些估算在写方案时非常有用，可以快速判断参数是否合理，再做代码仿真验证。

2. 从耦合模公式到 MATLAB 实现

2.1 光栅离散化与啁啾建模

写代码的第一步，是把连续的光栅离散成M段。这里的关键不是简单地把长度等分，而是每一段都要有自己独立的布拉格波长。

我把光栅从 z=0 到 z=L 均匀切M段，第 i 段中心位置记作 z_i，则该段的布拉格波长为：

λ_B(z_i) = λ₀ + C × (z_i - L/2)

C 的单位要先用 SI 制，即从 nm/cm 换算成 m/m。换算注意：如果 C_input = 0.5 nm/cm，那么 C_SI = 0.5 × 1e-9 / 0.01 = 5e-8，这个数是布拉格波长沿光栅长度方向的斜率。

然后在每一段用局部布拉格波长计算失谐量：

δ_i(λ) = 2π × n_eff / λ - π / Λ_i = 2π × n_eff × (1/λ - 1/λ_B(z_i))

其中 λ_B(z_i) = 2 × n_eff × Λ_i。

这样啁啾的物理效果就自然进来了：每个小段的布拉格波长不同，同一个仿真波长 λ 在不同段里的失谐量不同，于是反射强度分布沿光栅错开，长波在前端反射、短波在后端反射，时延差由此产生。

2.2 均匀光栅传输矩阵的代码实现

每一小段均匀光栅的传输矩阵长这样：

function T = uniform_grating_T(delta, kappa, dz) % 均匀FBG段的传输矩阵 % delta: 失谐量 (rad/m) % kappa: 耦合系数 (rad/m) % dz: 段长度 (m) gamma = sqrt(kappa.^2 - delta.^2); sg = sinh(gamma * dz); cg = cosh(gamma * dz); T = [ cg - 1i*delta/gamma*sg, -1i*kappa/gamma*sg ; 1i*kappa/gamma*sg, cg + 1i*delta/gamma*sg ]; end

注意几个细节：

• gamma 可能是复数。当 |delta| 大于 kappa 时，gamma 是纯虚数，sinh 和 cosh 的复数参数会自动处理，不用特殊判断。MATLAB 对复数 sinh/cosh 支持很好。
• kappa 的定义必须是 弧度/米。常用的近似是 κ = π × Δn / λ_B，这里 Δn 是折射率调制幅度。实际一些资料还会加入条纹可见度系数 v，通常取1即可。
• 这个矩阵是在“局部坐标”下把输入波映射到输出波。多个段连乘时顺序不能反。

2.3 整段光栅的响应计算

从入射端到出射端，逐段乘矩阵：

T_total = eye(2); for i = 1:M % 获取该段局部参数 lambdaB = lambdaB_z(i); delta = 2*pi*n_eff * (1/lambda - 1/lambdaB); kappa = pi * dn / lambda; % 也可以固定用 lambda0 T_local = uniform_grating_T(delta, kappa, dz); T_total = T_local * T_total; % 注意顺序，不能写反 end

边界条件是：入射端只有前向波，出射端没有反向波。设入射波为1，出射端反向波为0，推一下反射系数：

r = -T_total(2,1) / T_total(2,2)

这个表达式很多教材直接用，但建议自己推一遍。推导思路是列方程组，利用 [A(L); 0] = T × [1; B(0)]，由第二行解出 B(0)，再除以入射波1。

反射率就是 |r|²。得到了复反射系数 r，相位和群时延也就有了。

2.4 切趾函数是怎么加进去的

实际啁啾光栅如果不做切趾，光栅两端折射率的突变会产生类似法布里-珀罗腔的效应，反射谱带内会出现严重的振荡，时延纹波也会非常大。这个问题在实测中很常见，所以仿真阶段就要把切趾考虑进去。

最简单的做法是给每段的 kappa 乘上一个窗函数值。比如高斯切趾：

% 高斯切趾，sigma 控制切趾宽度 profile = exp(-((z - L/2).^2) / (2*sigma_z^2)); kappa_i = kappa0 * profile(i);

余弦切趾也常用：

profile = cos(pi * (z - L/2) / L);

切趾之后，光栅两端的耦合强度逐渐降到0，反射谱旁瓣会被大幅压低，时延纹波也会明显改善。代价是反射谱带宽边缘会变缓、有效反射带宽略有缩小，实际设计时要权衡。

3. 完整仿真流程与结果解读

3.1 一套可以直接跑的参考代码

下面这整套代码结构，我建议你直接抄到一个 .m 文件里跑一下，然后根据实际需求改参数。

clear; clc; close all; % ========== 光栅参数 ========== lambda0 = 1550e-9; % 中心波长，m n_eff = 1.45; % 有效折射率 dn = 1.5e-4; % 折射率调制深度 L = 0.10; % 光栅长度，m（10cm） C_input = 0.5; % 啁啾量，nm/cm M = 1000; % 分段数 % 换算啁啾量为 SI 制，d(lambdaB)/dz，无量纲 C_SI = C_input * 1e-9 / 0.01; % ========== 波长扫描范围 ========== BW = C_input * (L*100); % 总带宽 nm lambda_scan = (lambda0 - BW/2*0.9e-9) : 1e-12 : (lambda0 + BW/2*0.9e-9); lambda_scan = lambda_scan(:); % 列向量 % ========== 预计算位置与每段布拉格波长 ========== z = linspace(0, L, M+1); z_mid = (z(1:end-1) + z(2:end)) / 2; % 每段中点 lambdaB_z = lambda0 + C_SI * (z_mid - L/2); % 切趾：高斯型 sigma_z = L / 3.5; apod = exp(-((z_mid - L/2).^2) / (2*sigma_z^2)); % ========== 波长循环计算 ========== r = zeros(size(lambda_scan)); % 复数反射系数 for k = 1:length(lambda_scan) lambda = lambda_scan(k); T_total = eye(2); for i = 1:M % 局部失谐 delta = 2*pi*n_eff * (1/lambda - 1/lambdaB_z(i)); % 耦合系数，并乘切趾 kappa = pi * dn / lambda * apod(i); gamma = sqrt(kappa^2 - delta^2); sg = sinh(gamma * L/M); cg = cosh(gamma * L/M); T_i = [ cg - 1i*delta/gamma*sg, -1i*kappa/gamma*sg ; 1i*kappa/gamma*sg, cg + 1i*delta/gamma*sg ]; T_total = T_i * T_total; end r(k) = -T_total(2,1) / T_total(2,2); end % ========== 反射谱与群时延 ========== R = abs(r).^2; phase_r = unwrap(angle(r)); % 群时延，单位 ps c0 = 2.99792458e8; tau = -lambda_scan.^2 / (2*pi*c0) .* gradient(phase_r, lambda_scan); tau_ps = tau * 1e12; % ========== 画图 ========== figure('Position', [100 100 800 600]); subplot(2,1,1); plot(lambda_scan*1e9, 10*log10(R), 'LineWidth', 1.2); xlabel('波长 (nm)'); ylabel('反射率 (dB)'); title('啁啾光栅反射谱'); grid on; subplot(2,1,2); plot(lambda_scan*1e9, tau_ps, 'LineWidth', 1.2); xlabel('波长 (nm)'); ylabel('群时延 (ps)'); title('啁啾光栅群时延'); grid on;

这里波长扫描步长取了 1 pm（1e-12 m），对10nm带宽来说需要约9000个点，每个点循环1000段，MATLAB跑起来可能要几十秒。如果觉得慢，可以先把步长放宽到5pm甚至10pm看趋势，确定没毛病再用细步长出图。

3.2 反射谱结果怎么看

跑出来的反射谱，正常情况下是一个带内相对平坦的矩形谱，带外快速衰减。

如果一切正常，你会看到：

• 反射谱带宽接近设计值：0.5 nm/cm × 10 cm = 5 nm。
• 中心波长在1550 nm附近，但由于折射率调制会略微偏移，这个是正常的。
• 加切趾后，谱边缘会比较圆滑，带内振荡很小。如果没加切趾，带内会出现周期性的波动。

有一个细节要注意：反射率最大值不一定等于1。反射率峰值取决于折射率调制深度 Δn 和光栅长度 L 的乘积。如果 Δn 太小，光栅太短，反射率可能只有70%、80%。算一下耦合强度 κL 就能估算：

κL = π × Δn / λ_B × L

当 κL 在2到3时，反射率已经接近饱和。如果发现反射率太低，优先加 Δn 而不是加 L。

3.3 时延曲线和色散量怎么算

群时延曲线是这个仿真最重要的输出。对理想线性啁啾光栅，时延-波长曲线应该接近一条直线。

从图形上可以读出：

• 长波长端时延大，短波长端时延小（取决于啁啾方向，代码里默认长波对应光栅前端，短波对应后端）。
• 时延跨度约 2·n_eff·L/c。比如 L=10cm 时大约是 0.97 ns，和理论估算一致。
• 斜率就是色散量。对时延曲线做线性拟合，斜率 × 1e3（把nm换成pm）得到 ps/nm。

在MATLAB里可以用 polyfit 来拟合并计算色散：

% 只取反射率较高区域拟合 idx = R > 0.5; p = polyfit(lambda_scan(idx)*1e9, tau_ps(idx), 1); disp(['色散量: ', num2str(p(1)), ' ps/nm']);

拟合出的斜率就是色散量。假如是 -200 ps/nm，负号表示反常色散，可以用来补偿标准单模光纤的正色散。这个数值和设计估算要对得上，对不上就是参数有误。

4. 踩坑总结与常见问题排查

4.1 反射谱出现严重振荡，这是什么原因

如果跑出来的反射谱带内像锯齿一样上下跳动，最可能的原因是没加切趾或者切趾太弱。光栅两端折射率突变产生的边界反射和带内反射形成干涉，就会在谱上叠加振荡。

解决方法按优先级排列：

1. 先加高斯或余弦切趾，把两端的耦合强度压下去。
2. 检查段数M是否太少。啁啾光栅总带宽越大，每段的布拉格波长变化越大，需要更多段来逼近连续变化。1000段通常够，2000段更稳。
3. 检查波长扫描步长是否太粗，至少在带宽内扫500个点以上。

还有一种情况是 kappa 或 delta 单位搞错了，导致某些段处于“过耦合”状态，谱形也会异常。建议先拿均匀光栅（C=0）验证代码，和解析解对比反射谱峰值位置和带宽，确认基础公式没问题之后再引入啁啾。

4.2 群时延曲线出现很多毛刺

时延毛刺绝大多数来自相位处理的问题。反射系数 r 的相位在波长变化时会绕过 2π 边界，如果不做 unwrap，差分出来的时延就是一大堆尖峰。

代码里已经写了 unwrap(angle(r))，但还有另一个坑：当反射率很低时（比如谱边缘），相位本身非常敏感，数值噪声会被放大，时延曲线在谱边缘疯狂抖动。

处理办法是：只看反射率高于某一阈值的波长区域，比如保留 R > 0.3 或 R > 0.5 的部分，边缘区域直接扔掉不分析。群时延纹波统计也应在有效反射带宽内做，否则会把谱边缘的随机噪声算进去，得到毫无意义的结论。

另外，如果布拉格波长梯度很大，相位在相邻波段变化剧烈，gradient 差分可能不够精确。这时可以加密波长步长，或者用中心差分代替一阶差分。

4.3 带宽不对，和设计值差很多

带宽由啁啾量×长度决定，但很多人换算单位时出错。

你的啁啾量如果是 0.5 nm/cm，长度是 10 cm，带宽就是 5 nm。但实际如果你用 0.5 nm/cm 当作 SI 单位直接参与计算，算出带宽就会差10万倍。这个坑我踩过不止一次。

建议统一这样做：先算理论带宽，再设置扫描范围。代码里先用：

BW = C_input * (L*100); % 带宽 nm

然后扫描范围设为 lambda0 ± BW/2×0.9，留出10%的余量。这样可以快速发现参数设置是否出现数量级错误。

另一个容易出错的是局部布拉格波长的位置偏移。如果代码里没用 z_mid（每段中点），而是用了段起点，啁啾中心位置会整体偏移，反射谱中心也会偏移，带宽可能看起来不对称。用段中点计算局部布拉格波长更标准。

4.4 计算效率太低，怎么优化

初期调试可以降低段数和波长点数，等参数都对了再跑精细版本。如果确实需要高精度扫描，有几个优化手段：

• 预计算可以用矩阵形式避免最内层重复计算三角函数。把每个波长对应的 delta 算成 M×N 矩阵，然后按列循环，MATLAB向量化后能快不少。
• 使用 parfor 把波长循环并行化。MATLAB的并行计算工具箱可以轻松提速，但要注意循环内不能有依赖。
• 切趾和布拉格波长向量先算好，不要放在波长循环里重复计算。

我自己常用的做法是：先用 M=500、波长点 2000 快速调参，最后用 M=2000、波长点 10000 出正式结果。这样既能快调参数，又能保证最终图件精度。

4.5 实测光栅和仿真对不上

如果你做的是实际刻写光栅的仿真验证，还有一个常见差异来源：实际光栅的折射率调制深度不是均匀的，切趾分布也不是理想函数。此外，光纤材料色散在窄带例子中影响小，但在宽带（几十nm）情况下，n_eff 随波长变化就会带来明显偏差。

建议在仿真中加入材料色散近似，至少要考虑纤芯和包层折射率随波长变化导致的有效折射率变化。更简单的做法是把 n_eff 设为与波长相关的多项式拟合，哪怕只在1550nm附近取一阶斜率，也能显著改善宽带时延曲线的吻合度。

最后说一个实际设计中的经验：判断一组啁啾光栅参数是否合理，先看反射谱带宽和时延斜率这两个宏观指标，它们准了，再抠纹波和边沿形状。我见过太多人一上来就纠结反射谱边缘不够陡，其实在色散补偿应用里，带内线性度和纹波才是真正决定系统性能的指标。这个顺序搞反了，后面全是无用功。

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