简介：这份资源面向光纤通信与光学传感方向的学习者和科研人员，聚焦光纤布拉格光栅（FBG）的反射谱与透射谱模拟，帮助理解布拉格条件、有效折射率与光栅周期对光谱响应的影响。压缩包内共1个文件，为MATLAB脚本（.m），整体仅614B，体量轻巧，便于直接运行与二次修改。脚本围绕初始化光栅参数、计算各波长反射与透射系数、借助傅里叶变换处理周期性结构并绘制谱线等环节展开，读者可通过调整周期、折射率调制深度和光栅长度，直观观察反射峰位置与谱宽变化。目前已有308人学习，适合作为FBG原理验证与参数扫描的入门工具，也可为后续传感与滤波设计提供参考。
1. 从一份 FBG.rar 说起：反射谱和透射谱到底怎么在 MATLAB 里跑出来
拿到一个叫FBG.rar的压缩包，里面躺着几份 MATLAB 脚本，目标很明确：把光纤布拉格光栅的反射谱和透射谱画出来。这件事听起来像课程作业，但真正做过的人都知道，它是后续做波长解调、传感标定、滤波器设计的第一块地基。FBG 的本质是一段周期性折射率调制，满足布拉格条件的波长会被反射回去，其余透射过去，所以反射谱和透射谱是一体两面。MATLAB 在这里的角色不是“画图工具”，而是把耦合模方程或者传输矩阵法算成一条可用的曲线。适合谁看？做光纤传感、光通信器件仿真、或者刚接手一个 FBG 解调项目，需要先自己把谱形跑通的人。下面不聊虚的，直接从参数、代码、踩坑一路推到能复现。
2. 先把 FBG 的谱算对：耦合模方程与传输矩阵的选型理由
2.1 为什么大多数 MATLAB 脚本用传输矩阵法而不是直接解微分方程
FBG 的反射谱计算，常见做法有两类：一类是直接数值求解耦合模方程，另一类是传输矩阵法。直接解微分方程在均匀光栅里没问题，但一旦要做切趾、相移、啁啾，边界条件会变得很烦。传输矩阵法的好处是把光栅切成 N 段，每段用一个 2×2 矩阵表示，最后把所有矩阵乘起来，得到整个光栅的反射系数和透射系数。这个思路在 MATLAB 里特别顺手，因为矩阵乘法是它的强项，而且分段数 N 可以随参数调整，精度和速度之间好平衡。
我一般会先确认三个物理量：中心波长、折射率调制深度、光栅长度。中心波长决定了谱峰位置，调制深度决定了反射率峰值和带宽，光栅长度决定了整体反射强度。很多网上流传的脚本只给一个lambda_B就开算，结果谱形对不上，问题往往出在没把有效折射率和调制深度分开设置。
2.2 传输矩阵法的核心参数表
在动手写代码前，先把参数含义和典型值列清楚。下面这张表是我在多个 FBG 项目里反复用到的起点，具体数值按你的实际光纤和写入条件改。
| 参数名 | 含义 | 典型值 | 对谱形的影响 |
|---|---|---|---|
lambda_B | 布拉格中心波长 | 1550 nm | 决定反射峰位置 |
n_eff | 有效折射率 | 1.447 | 与周期共同决定中心波长 |
delta_n | 折射率调制深度 | 1e-4 ~ 5e-4 | 越大反射率越高，带宽越宽 |
L | 光栅长度 | 5 mm ~ 20 mm | 越长反射率越高，峰越窄 |
N | 分段数 | 200 ~ 1000 | 越大越精确，但计算变慢 |
apodization | 切趾函数 | 高斯/升余弦 | 抑制旁瓣，改变谱形 |
这张表里最容易翻车的是delta_n和L的配合。有人把delta_n设成 1e-3，L又给 20 mm，算出来反射率直接饱和到 1，透射谱掉到零，看起来“完美”，但实际光栅根本做不到那么强的调制。血泪经验是：先固定L，再调delta_n，让峰值反射率落在 0.5 到 0.9 之间，这样谱形才有观察价值。
2.3 用 MATLAB 跑通均匀 FBG 的最小代码
下面这段代码是均匀 FBG 的传输矩阵法实现，直接复制到 MATLAB 里就能跑。它不依赖任何工具箱，纯基础语法。
% 均匀FBG传输矩阵法计算反射谱与透射谱 clear; clc; % 物理参数 lambda_B = 1550e-9; % 中心波长，单位米 n_eff = 1.447; % 有效折射率 delta_n = 2e-4; % 折射率调制深度 L = 10e-3; % 光栅长度，10 mm N = 500; % 分段数 % 波长扫描范围 lambda = linspace(1548e-9, 1552e-9, 2000); % 计算光栅周期 Lambda = lambda_B / (2 * n_eff); % 耦合系数（简化形式，适用于弱调制） kappa = pi * delta_n / lambda_B; % 初始化反射系数 r = zeros(size(lambda)); for idx = 1:length(lambda) lam = lambda(idx); % 相位失配量 delta = 2 * pi * n_eff * (1/lam - 1/lambda_B); % 分段长度 dz = L / N; % 传输矩阵初始化 M = eye(2); for k = 1:N % 每段的传输矩阵 gamma = sqrt(kappa^2 - delta^2 + 0i); if abs(gamma) < 1e-12 % 避免除零，用极限形式 T11 = 1 - 1i * delta * dz; T12 = -1i * kappa * dz; T21 = 1i * kappa * dz; T22 = 1 + 1i * delta * dz; else T11 = cosh(gamma * dz) - 1i * (delta / gamma) * sinh(gamma * dz); T12 = -1i * (kappa / gamma) * sinh(gamma * dz); T21 = 1i * (kappa / gamma) * sinh(gamma * dz); T22 = cosh(gamma * dz) + 1i * (delta / gamma) * sinh(gamma * dz); end M = M * [T11, T12; T21, T22]; end % 边界条件：入射端反射为0，出射端无反射 r(idx) = -M(2,1) / M(2,2); end % 反射率和透射率 R = abs(r).^2; T = 1 - R; % 绘图 figure; plot(lambda*1e9, R, 'b', 'LineWidth', 1.5); hold on; plot(lambda*1e9, T, 'r--', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('反射率 / 透射率'); legend('反射谱', '透射谱'); title('均匀FBG反射谱与透射谱'); grid on;这段代码的逻辑说明：外层循环扫描波长，内层循环把光栅分成 N 段，每段用传输矩阵表示，最后累乘得到总矩阵。gamma是传播常数，当它接近零时用泰勒展开避免除零错误。r是反射系数，取模平方得到反射率，透射率用 1 减反射率。参数方面，N取 500 在普通笔记本上大约几秒跑完，如果谱形有毛刺，先加到 1000 看看是不是分段不够。lambda的范围要覆盖中心波长前后至少 2 nm，否则看不到完整的旁瓣。
2.4 切趾和啁啾怎么加进去
均匀光栅的反射谱旁瓣很高，实际器件往往要做切趾。切趾的本质是让delta_n沿光栅长度变化，比如高斯切趾：delta_n(z) = delta_n * exp(-((z - L/2)/(L/4))^2)。在传输矩阵法里，只需要在每段的kappa上乘这个系数。啁啾则是让周期Lambda沿长度线性变化，对应lambda_B随位置偏移。这两步加进去后，谱形会明显不同：切趾压旁瓣，啁啾展宽带宽。我一般会先跑均匀的，确认代码没问题，再逐步加切趾和啁啾，这样出问题容易定位。
3. 从脚本到可复用函数：把 FBG 反射谱计算封装成 MATLAB 函数
3.1 为什么要封装成函数而不是每次改脚本
直接改脚本的后果是，每换一组参数就要复制一份新文件，最后文件夹里全是fbg_test1.m、fbg_test2.m，自己都分不清哪个是最终版。封装成函数的第一个好处是参数通过输入传递，调用时一目了然；第二个好处是可以批量扫描参数，比如扫delta_n从 1e-4 到 5e-4，看反射率峰值怎么变。MATLAB 的函数文件写起来不复杂，关键是把物理参数和计算逻辑分开。
3.2 函数接口设计与参数校验
下面是一个我常用的函数签名，输入是结构体params，输出是波长、反射率、透射率三个数组。结构体的好处是参数名自带说明，调用时不容易传错位置。
function [lambda, R, T] = fbg_spectrum(params) % FBG_SPECTRUM 计算光纤布拉格光栅的反射谱和透射谱 % 输入：params 结构体，包含以下字段： % lambda_B - 中心波长 (m) % n_eff - 有效折射率 % delta_n - 折射率调制深度 % L - 光栅长度 (m) % N - 分段数 % lambda_min - 扫描起始波长 (m) % lambda_max - 扫描终止波长 (m) % num_points - 波长点数 % apod - 切趾类型：'none' | 'gaussian' | 'raisedcosine' % 输出： % lambda - 波长数组 (m) % R - 反射率数组 % T - 透射率数组 % 参数校验 if ~isfield(params, 'lambda_B') || ~isfield(params, 'n_eff') error('缺少必要参数：lambda_B 或 n_eff'); end if params.delta_n <= 0 || params.L <= 0 error('delta_n 和 L 必须为正数'); end if params.N < 50 warning('分段数 N 小于 50，结果可能不准确'); end % 解包参数 lambda_B = params.lambda_B; n_eff = params.n_eff; delta_n = params.delta_n; L = params.L; N = params.N; lambda = linspace(params.lambda_min, params.lambda_max, params.num_points); % 光栅周期 Lambda = lambda_B / (2 * n_eff); kappa0 = pi * delta_n / lambda_B; % 切趾系数 z = linspace(0, L, N); switch lower(params.apod) case 'gaussian' apod = exp(-((z - L/2) / (L/4)).^2); case 'raisedcosine' apod = 0.5 * (1 + cos(2*pi*(z - L/2)/L)); otherwise apod = ones(1, N); end r = zeros(size(lambda)); for idx = 1:length(lambda) lam = lambda(idx); delta = 2 * pi * n_eff * (1/lam - 1/lambda_B); dz = L / N; M = eye(2); for k = 1:N kappa = kappa0 * apod(k); gamma = sqrt(kappa^2 - delta^2 + 0i); if abs(gamma) < 1e-12 T11 = 1 - 1i * delta * dz; T12 = -1i * kappa * dz; T21 = 1i * kappa * dz; T22 = 1 + 1i * delta * dz; else T11 = cosh(gamma * dz) - 1i * (delta / gamma) * sinh(gamma * dz); T12 = -1i * (kappa / gamma) * sinh(gamma * dz); T21 = 1i * (kappa / gamma) * sinh(gamma * dz); T22 = cosh(gamma * dz) + 1i * (delta / gamma) * sinh(gamma * dz); end M = M * [T11, T12; T21, T22]; end r(idx) = -M(2,1) / M(2,2); end R = abs(r).^2; T = 1 - R; end逻辑说明：函数开头做参数校验，避免传错结构体字段导致莫名其妙的结果。切趾系数在循环外算好，循环内直接乘到kappa上，这样比每段重新算切趾函数快。apod字段用字符串控制，默认'none'就是均匀光栅。参数方面，num_points建议不低于 1000，否则谱峰可能被采样点跳过。N和num_points是两回事：N是光栅分段数，num_points是波长扫描点数，两个都要够大。
3.3 调用函数并对比不同切趾效果
封装好之后，调用就变成几行。下面这段代码同时算均匀、高斯切趾、升余弦切趾三种情况，画在一张图里对比。
% 对比不同切趾类型的反射谱 params.lambda_B = 1550e-9; params.n_eff = 1.447; params.delta_n = 2e-4; params.L = 10e-3; params.N = 500; params.lambda_min = 1548e-9; params.lambda_max = 1552e-9; params.num_points = 2000; apod_types = {'none', 'gaussian', 'raisedcosine'}; colors = {'b', 'r', 'g'}; figure; hold on; for i = 1:length(apod_types) params.apod = apod_types{i}; [lambda, R, ~] = fbg_spectrum(params); plot(lambda*1e9, R, colors{i}, 'LineWidth', 1.5); end xlabel('波长 (nm)'); ylabel('反射率'); legend('均匀', '高斯切趾', '升余弦切趾'); title('不同切趾类型的FBG反射谱对比'); grid on;运行后你会看到，均匀光栅的旁瓣明显，高斯切趾把旁瓣压下去但主峰略微展宽，升余弦切趾的旁瓣抑制更强。这里的关键参数是切趾函数的宽度，我用的L/4是经验值，如果旁瓣压不下去，可以把分母改小，让切趾更剧烈。但注意，切趾太强会导致有效光栅长度变短，反射率峰值下降，需要回头调delta_n或L补偿。
4. 避坑与排查：FBG 谱形计算里最容易翻车的 5 个地方
4.1 反射率大于 1 或者出现负值
现象：画出来的反射谱某些点超过 1，或者透射谱出现负数。原因通常是传输矩阵累乘时数值溢出，或者gamma在零附近处理不当。解决：检查delta_n是否过大导致kappa远超实际，把delta_n降到 5e-4 以下；同时在abs(gamma) < 1e-12的分支里，确认泰勒展开的符号没写反。另一个可能是N太小，分段粗糙导致矩阵乘积失真，把N加到 1000 试试。
4.2 谱峰位置对不上中心波长
现象：反射峰不在lambda_B上，偏了几纳米。原因多半是Lambda算错了。Lambda = lambda_B / (2 * n_eff)这个关系里，n_eff用的是有效折射率，不是纤芯折射率。如果你填的是 1.45 而实际有效折射率是 1.447，偏差会累积到谱峰上。解决：用lambda_B和Lambda反推n_eff，确保三者自洽。另外，波长扫描范围要覆盖lambda_B附近至少 ±2 nm，否则峰被截断。
4.3 切趾后反射率掉得太厉害
现象：加了高斯切趾，反射率从 0.9 掉到 0.3。原因：切趾函数把光栅两端的kappa压到接近零，等效于缩短了光栅长度。解决：要么增大delta_n或L补偿，要么把切趾函数的宽度放宽，比如把L/4改成L/3。我一般会先算均匀光栅的反射率，再算切趾后的，确保峰值下降不超过 20%。
4.4 循环太慢，跑一次要几分钟
现象：N取 1000、num_points取 5000 时，MATLAB 跑一次要好几分钟。原因：双层循环里每次都在算gamma和矩阵乘法，没有向量化。解决：把波长循环向量化，或者把内层分段循环用矩阵连乘代替。更实际的做法是先用N=200、num_points=1000快速看谱形，确认参数合理后再加细。如果经常要扫参数，考虑用parfor并行，但要注意 MATLAB 并行工具箱的配置。
4.5 透射谱和反射谱加起来不等于 1
现象：R + T在某些波长上明显偏离 1。原因：传输矩阵法假设无损耗，如果代码里不小心引入了复数的错误分支，或者r的模平方算错，就会破坏能量守恒。解决：检查r = -M(2,1)/M(2,2)这个边界条件是否用对，确认M是总传输矩阵而不是散射矩阵。另外，如果加了啁啾，lambda_B随位置变化，delta的表达式要相应修改，否则相位失配算错，能量也会不守恒。
5. 进阶技巧：用 MATLAB 做 FBG 参数扫描与谱形验证
5.1 扫折射率调制深度看反射率峰值变化
参数扫描是验证 FBG 设计是否合理的最快方式。下面这段代码扫delta_n从 1e-4 到 5e-4，每次算反射率峰值，画成曲线。
% 扫描折射率调制深度，观察反射率峰值变化 delta_n_list = linspace(1e-4, 5e-4, 20); R_peak = zeros(size(delta_n_list)); params.lambda_B = 1550e-9; params.n_eff = 1.447; params.L = 10e-3; params.N = 300; params.lambda_min = 1548e-9; params.lambda_max = 1552e-9; params.num_points = 1000; params.apod = 'none'; for i = 1:length(delta_n_list) params.delta_n = delta_n_list(i); [~, R, ~] = fbg_spectrum(params); R_peak(i) = max(R); end figure; plot(delta_n_list*1e4, R_peak, 'bo-', 'LineWidth', 1.5); xlabel('折射率调制深度 \times 10^{-4}'); ylabel('峰值反射率'); title('FBG峰值反射率随调制深度的变化'); grid on;逻辑说明：每次改变delta_n，调用fbg_spectrum算反射谱，取最大值存下来。参数方面，N可以降到 300 加快扫描，因为这里只关心峰值趋势，不关心旁瓣细节。跑完后你会看到反射率随delta_n单调上升，但上升速度逐渐变缓，这就是饱和效应。实际设计时，选在曲线拐点附近，既保证反射率够高，又不会对调制深度过于敏感。
5.2 用透射谱零点验证切趾效果
切趾效果好不好，看反射谱旁瓣是一种方法，看透射谱零点更直接。均匀光栅的透射谱在反射峰位置有一个深零点，切趾后零点变宽但深度可能变浅。下面这段代码提取透射谱最小值，对比不同切趾类型。
% 对比不同切趾类型的透射谱零点深度 apod_types = {'none', 'gaussian', 'raisedcosine'}; T_min = zeros(size(apod_types)); params.lambda_B = 1550e-9; params.n_eff = 1.447; params.delta_n = 2e-4; params.L = 10e-3; params.N = 500; params.lambda_min = 1549e-9; params.lambda_max = 1551e-9; params.num_points = 2000; for i = 1:length(apod_types) params.apod = apod_types{i}; [~, ~, T] = fbg_spectrum(params); T_min(i) = min(T); end disp('透射谱最小透射率：'); for i = 1:length(apod_types) fprintf('%s: %.4f\n', apod_types{i}, T_min(i)); end运行后，均匀光栅的T_min接近 0，高斯切趾的T_min可能升到 0.05 左右，升余弦切趾更高。这说明切趾在压旁瓣的同时，牺牲了主峰的反射深度。如果你的应用对透射零点深度敏感，比如做滤波，就要在切趾和反射率之间折中。我一般会先定一个可接受的T_min下限，再选切趾类型和宽度。
5.3 一个我常犯的错误：忘了改波长范围
最后说一个我踩过多次的坑。每次换lambda_B，比如从 1550 nm 换到 1530 nm，只改了lambda_B，忘了改lambda_min和lambda_max，结果画出来的谱是空的，或者只有一半。MATLAB 不会报错，因为linspace照样生成数组，只是范围不对。后来我养成了一个习惯：在函数里加一句检查，如果lambda_B不在[lambda_min, lambda_max]范围内，直接warning提示。这个习惯帮我省了很多来回改参数的时间。希望帮到你。
本文还有配套的精品资源，点击获取
