企业官网建设流程全解析

简介：利用分步傅立叶积分法计算光纤中色散与自相位调制（SPM）的MATLAB源码及说明文档，面向光纤通信方向的研究生、工程师或对非线性光纤光学感兴趣的开发者。资源包共2个文件，含1个.m脚本与1个txt说明，压缩包仅2KB，内容精简，可快速下载部署。代码采用经典的分步思想，将非线性薛定谔方程拆成线性与非线性两步处理，通过逐步迭代在频域中计算色散、在时域中计算非线性效应，能模拟不同波长分量速度差导致的脉冲展宽，以及强光脉冲引起的折射率变化与相位畸变；说明文档中还梳理了模式色散、材料色散、波导色散以及孤子效应的背景推导，便于对照MATLAB代码理解每一步频域和时域交替计算逻辑。已有554人学习下载，适合作为理解色散管理与SPM影响的入门代码参考，帮助读者掌握从光纤传输方程到数值实现的基本流程。

1. 光纤中色散与SPM仿真：从分步傅立叶积分法讲起

拿到一个“光脉冲在光纤里跑”的问题，色散和自相位调制（SPM）同时起作用，手推解析公式基本走不通。分步傅立叶积分法是目前工程上最顺手的解法：把非线性薛定谔方程沿传播方向切成许多小段，每段只让线性色散和非线性相位分别作用，再用 FFT 把偏微分变成频域乘法。一个皮秒级高斯脉冲跑 20 公里，也就是几十到几百次循环的事。适合刚接触光通信仿真的人，也适合要做参数扫描的熟手。这篇笔记从原理、MATLAB 实现、验证基准到五个常见翻车点按顺序过一遍，照着搭一遍就能改成自己的光纤参数和脉冲形状。

2. 分步傅里叶的算符拆分：为什么线性项在频域、非线性项在时域

在动手敲代码之前，先把要拆的东西说清楚。光脉冲在单模光纤里的传输，普遍用慢变包络近似下的非线性薛定谔方程描述。我习惯写成下面这种形式：

∂u/∂z = i(β2/2)∂²u/∂t² + iγ|u|²u

其中 u 是脉冲复包络，z 是传播距离，t 是随脉冲移动的本地时间，β2 是群速度色散系数，γ 是非线性系数。损耗项 αu/2 可以先不放进去，后面需要对比实验再加也不迟。

2.1 色散算符本质上是“全局的”，SPM 算符本质上是“局部的”

右边两项看着像两个作用力，但它们的数学性质完全不同。色散项是对时间求二阶导，在时域里它会把脉冲每一点的相位和附近点耦合起来；可一旦变换到频域，它变成一个纯粹的乘法：每个频率分量各自按-iβ2ω²/2传播，互不干扰。这是傅里叶变换最漂亮的地方——把一个偏微分算符对角化。

SPM 项则是瞬时的、局部的：t 时刻的相位变化只取决于该点自己的瞬时功率|u(t)|²。在频域里它反而会混叠成卷积，非常难处理；但在时域里它就是一个复指数乘法exp(iγ|u|²Δz)。

如果同时处理这两项，既没有闭式解，有限差分又受数值色散拖累。分步傅里叶的思路是：既然色散在频域好算、非线在时域好算，那就把光纤切成小段，每段内部把它们的耦合暂时忽略，分开算再拼起来。

2.2 对称分步格式：半个色散步 + 一个非线性步 + 半个色散步

工程上几乎不用最朴素的“先色散后非线性”的一阶格式，因为那样每个步长只精确到一阶。比较稳妥的是对称分步格式（Strang 式分裂），在一个步长 h 内这么安排：

1. 先让色散单独作用 h/2；
2. 再让 SPM 单独作用 h；
3. 最后让色散单独作用 h/2。

写成算符就是L(h/2) · N(h) · L(h/2)，其中色散算符 L 在频域实现，非线性算符 N 在时域实现。这样整个传播算法的局部误差是 O(h³)，累积到总长度上就是 O(h²)，比简单分裂高一个量级，而实现代价只多了一次 FFT。

这里有一个容易误解的点：为什么两个半色散步合起来还是一个整色散步？因为线性算符在频域是精确指数，L(h/2)·L(h/2) = L(h)。但每轮循环结束时不急着把它合并，保持这种对称结构，程序的可读性和扩展性都会好一些。

2.3 步长先按物理尺度定：色散长度、非线性长度

步长 h 不能拍脑袋。工程上有两个参考长度：

LD = T0² / |β2|，色散长度，表示脉冲因色散明显展宽的特征距离；
LNL = 1 / (γ P0)，非线性长度，表示 SPM 相位积累到 1 rad 量级的特征距离。

SSFM 的步长要远小于二者中的较小者。我一般取h < min(LD, LNL)/20。举个例子：T0 = 10 ps，β2 = -2.17e-26 s²/m，那么 LD ≈ 4.6 km；P0 = 100 mW，γ = 1.3e-3 W⁻¹m⁻¹，那么 LNL ≈ 7.7 km。两者取小是 4.6 km，除以 20 是 230 m，所以取 50～100 m 都很安全。

如果只盯着总长度 20 km，固定步长取 100 m，循环 200 次，每次 2~3 次 FFT，毫秒级就出结果。但如果功率很高，非线性相位在单个步长内累积超过 π，分裂误差会立刻浮现，这个在第五章会专门讲。先记住：步长不要用总长度除以一个好看的整数，而是用物理尺度去约束。

3. MATLAB 实现与参数表：一套可跑的色散 + SPM 仿真骨架

第 2 章把原理说透了，这一章直接上能跑的代码。为了不让单位换算变成坑，我把所有量纲统一到米、秒、瓦特。很多仿真翻车都发生在“ps²/km”和“s²/m”混用上，这里先花两分钟把单位钉死。

3.1 参数定义：所有单位统一到 SI 制

以下是参数表，对应 1550 nm 窗口的标准单模光纤：

MATLAB 代码第一段：

%% 物理常数与光纤参数 c = 3e8; lambda = 1550e-9; % 中心波长，m D = 17e-6; % 17 ps/(nm*km) 换算为 s/m^2 beta2 = -lambda^2/(2*pi*c)*D; % 群速度色散，s^2/m gamma = 1.3e-3; % 非线性系数，W^-1*m^-1 alpha = 0; % 线性损耗，Np/m，先置零 %% 时间网格与频率网格 N = 2^15; Twin = 1280e-12; % 时间窗总宽度，s dt = Twin/N; t = (-N/2:N/2-1)*dt; % 时间坐标，s freq = (-N/2:N/2-1)/(N*dt); % 频率坐标，Hz omega = 2*pi*freq; % 角频率，rad/s %% 输入脉冲 T0 = 10e-12; P0 = 0.1; u0 = sqrt(P0)*exp(-t.^2/(2*T0^2)); % 高斯脉冲，峰值功率为 P0

这段代码里的 D 从 17 ps/(nm·km) 换算成 SI 制是 17e-6 s/m²，别怀疑这个数值。β2 那个负号来自反常色散，1550 nm 标准光纤的 β2 本来就是负的；如果符号写反，色散作用就从展宽变成压缩，输出形态会非常奇怪。

时间窗 1280 ps 是 10 ps 脉宽的 128 倍，对固定步长仿真来说足够宽。后面第五章会讲窗口边缘的“折返”问题，这里先留个印象：窗口宁可大一点，FFT 点数多花不了多少时间，结果的可信度提升却是实打实的。

另外注意t的构造用了(-N/2:N/2-1)*dt，不是linspace(-Twin/2, Twin/2, N)。两者相差一个采样点，这个细节直接决定频域数组和 FFT 输出能不能对齐。

3.2 传播主循环：对称分步格式的实现

有了u0和频率轴，下一步构造频域线性算子。这里有一个非常容易踩的索引坑：我用fftshift(fft(u))把零频搬到数组中心，那么omega数组必须也按“负频率到正频率”排列，代码里omega = 2*pi*freq已经满足这个顺序。下面这段是主循环：

%% 色散线性算子：半步步长 lin_op = exp(-1i * beta2 * omega.^2 * h/2); if alpha ~= 0 loss_op = exp(-0.5 * alpha * h); else loss_op = 1; % 无损耗时保持单位为 1 end %% 循环传播 u = u0; Nz = round(20e3 / h); save_step = round(Nz/10); % 大概每 2 km 记录一次 for k = 1:Nz % 前半步线性色散 uh = fftshift(fft(u)); uh = uh .* lin_op; u = ifft(ifftshift(uh)); % 损耗与非线性的时域作用 u = u .* loss_op; u = u .* exp(1i * gamma * abs(u).^2 * h); % 后半步线性色散 uh = fftshift(fft(u)); uh = uh .* lin_op; u = ifft(ifftshift(uh)); end

逻辑并不复杂：每轮先做 h/2 的色散，在频域乘一个指数；回到时域乘 SPM 相位；再回到频域做剩下的 h/2 色散。abs(u).^2是当前时刻每个采样点的瞬时功率，gamma*abs(u).^2*h就是这一小段累积的非线性相位，整个过程没有对脉冲波形做任何“截断”或“近似归一化”，这就是 SSFM 最让人放心的地方。

lin_op只依赖于频率轴和 h，和距离无关，所以放在循环外先算好。如果后面改成自适应步长，lin_op就要在步长变化时更新，这一点在第六章还会提。loss_op放在时域乘，对应的是功率损耗，所以指数上是-0.5*alpha*h，不是-alpha*h。

3.3 每段距离的记录与能量守恒检查

很多仿真跑完只看最后波形，中间过程一概不管。我一般会在循环里顺带记录功率和能量，这样即使结果错了，也能知道从哪一段开始偏的。

%% 预分配记录变量 z_axis = linspace(0, 20e3, 11); record_idx = round(z_axis / h) + 1; P_record = zeros(11, N); E_record = zeros(1, 11); %% 循环内，在 k 满足记录条件时执行 % 假设上面主循环中加一句： % if any(k == record_idx) % P_record(num, :) = abs(u).^2; % E_record(num) = sum(abs(u).^2) * dt; % end

记录点一般取 10～20 个即可，没必要每步都存，否则脉冲窗口很大的时候内存会吃紧。能量守恒是检验数值格式是否稳定的第一道闸门：没有损耗和增益时，sum(abs(u).^2)*dt应该恒定不变。如果能量漂移超过千分之一，先怀疑步长太大，再检查 FFT 的归一化有没有被自己改坏。

4. 验证从线性基准开始：解析解、能量守恒与 SPM 频谱展宽

写完代码不能直接拿非线性结果去交差，先做两个对照实验，分别验证色散步和非线性步的正确性。这是过去我被导师要求养成的习惯，现在我觉得这是最快定位 bug 的方式。

4.1 线性情形：高斯脉冲的解析展宽公式

先把 γ 置 0，让程序退化成纯色散问题。高斯脉冲在色散介质中的传播是教科书上唯一有闭式解的例子里最适合拿来当基准的一个。传播到距离 L 后，理论上：

u(L,t) = u0 * T0 / sqrt(T0^2 + i*beta2*L) * exp(-t^2 / (2*(T0^2 + i*beta2*L)))

这个公式在代码里可以这么对着验证：

%% 无非线性对照 gamma0 = 0; % 用上面主循环跑完 20 km，得到 u_lin %% 解析公式 u_theory = sqrt(P0)*T0/sqrt(T0^2 + 1i*beta2*20e3) ... .* exp(-t.^2./(2*(T0^2 + 1i*beta2*20e3))); % 对比输出端时域波形 max(abs(u_lin - u_theory).^2)

这个对比误差应该在 1e-6 量级以下，因为对称式 SSFM 对线性步的频域乘法是精确的，唯一误差来自中间非线性算符被强制设为 1。如果这里有明显偏差，检查频率轴排序和lin_op的符号就够了。

在 20 km 末端，10 ps 高斯脉冲的峰值功率会掉到原来的1/sqrt(1 + (beta2*L/T0^2)^2)左右。代入 β2 = -2.17e-26，L = 2e4，T0² = 1e-22，比值为 4.34，算下来峰值功率大约掉到 0.1/√(1+4.34²) ≈ 0.022 W，也就是约 -6.5 dB。这个数可以作为快速心算的参考，跑完代码如果差太多，大概率是单位换算出了偏差。

4.2 打开 SPM：频谱展宽与边瓣结构

线性基准通过后，把 γ 恢复到 1.3e-3，同样跑 20 km。由于 β2 为负且 |β2| 和 γP0 都处于同一量级，SPM 与色散之间会出现明显的相互作用，而不是简单的“先展宽再变谱”。

观察频谱时，把输出端的光场做 FFT：

%% 输出频谱 U_out = fftshift(fft(u)); P_spectrum = abs(U_out).^2 * dt; plot(fftshift(freq), P_spectrum);

峰值非线性相位可以提前估算：phi_max = gamma * P0 * L，代入 γP0 = 1.3e-4（1/m），L = 2e4 m，得到 2.6 rad。这个相位量级足够让光谱明显开口展宽，并在频谱两侧出现振荡边瓣。边瓣是 SPM 导致的“瞬时频率随时间变化”在频域里的投影，不是算法噪声。判断方法是：把 γ 慢慢增大，边瓣间距按连续规律变化，而如果时间窗不够大，看到的伪结构往往是突然出现的周期性尖峰。

这里我习惯把线性结果和 SPM 结果的频谱画在同一张图上。正常情况是：中心频率附近的主瓣保持类似前 1/e² 半径的形状，两侧对称长出新的峰，整个谱宽比线性情形展宽到 3～5 倍。如果不对称，先怀疑 β2 符号；如果展宽几十倍，多半是步长不够小。

4.3 用功率扫描检查非线性强度与理论量级是否一致

还有一个快速自检方法：把峰值功率从 10 mW 扫到 200 mW，记录输出光谱的 99% 能量带宽。理论上功率每增加一倍，非线性长度缩一半，SPM 效应显著增强，谱宽会明显变宽。如果程序输出对功率完全不敏感，通常是在某个地方不小心对u做了归一化，导致abs(u).^2不再是真实物理功率。

这个扫描不需要重新编译，只需要在参数段循环改 P0，每次重新构造u0即可。扫描结果可以用一张峰值功率 - 输出谱宽的曲线展示，比单看一眼波形更能说明代码对非线性项的响应是正确的。

5. 避坑手册：分步傅里叶仿真常见的五个翻车现场

这段是过去实际踩过的坑汇总。每条都按现象、原因、解决三个环节写，方便你直接对照自己的输出定位。

5.1 脉冲尾巴跑到时间窗口边缘，被“折返”回窗口开头

现象：时域波形两侧出现原脉冲的影子，整体看起来像是两个脉冲叠在一起，频谱上出现等间距细密条纹。

原因：FFT 隐式假定时域信号是周期延拓的。脉冲色散展宽后超出Twin，从窗口右端消失的部分又从左端进入，形成虚假干涉。口算时觉得 1280 ps 已经很大，但 20 km 后高斯脉冲 RMS 展宽了好几倍，仍然可能撞边。

解决：把时间窗扩大。最稳妥的做法是先跑一次线性情形看输出脉冲的实际包络范围，让Twin至少是输出端脉冲有效宽度的 4 倍。如果为了省内存不扩窗口，就改用局部色散吸收边界，但工程上还是扩窗口最简单，N 从 2^15 涨到 2^16，耗时几乎没变化。

5.2 结果对步长极其敏感，缩小 h 后频谱完全变样

现象：用 200 m 步长跑出来频谱毛糙，把步长改成 10 m 后边瓣数目和位置都变了，能量也不守恒。

原因：单个步长内的非线性相位旋转超过了安全范围。步长 h 内 SPM 相位增量为γ|u|²h，当它接近或超过 π，分裂近似就失效了。此时不是“步长越小越精确”这么简单，而是整个格式在小步长范围内才开始进入可信区。

解决：除了用min(LD, LNL)/20约束，再额外加一条判断：gamma * max(abs(u).^2) * h < 0.1。如果这条不满足，就缩步长。等脉冲峰值功率降下来之后，后半程可以适当放大步长，这个策略放到最后一章讲。

5.3 β2 单位和符号错乱，脉冲在反常色散区被“压缩”或相位翻转

现象：原本应该展宽的脉冲，模拟结果却出现脉宽变窄；或者把正色散光纤跑出负色散的特征。

原因：最常见的是 D 到 β2 的换算错误。我用一个最直接的公式避开黑匣子：beta2 = -lambda^2/(2*pi*c)*D，其中 D 必须是 s/m²、λ 必须是 m。很多人直接拿 17 代入并写成“ps²/km”，最后数值差了三个数量级。

解决：固定一套单位换算函数，在代码开头写清楚。每次新建仿真脚本都直接复制这段参数定义，不要手敲数值。同时先跑第 4.1 节的线性解析对照，β2 符号错了会在第一个对比图上立刻暴露。

5.4 FFT 频率轴索引没对齐，输出频谱左右镜像或高频毛刺

现象：SPM 效果越强，光谱图上高频部分越像“雪花”，或频谱左右不对称到离谱；有时候明明加了滤波器，结果完全不是那么回事。

原因：fft和fftshift的顺序、频率轴数组顺序没有同步。常见做法是omega = fftshift(2*pi*(-N/2:N/2-1)/(N*dt))，但这个写法很容易让零频位置和fftshift(fft(u))错位。

解决：用本章第 3.1 节里的固定组合：freq = (-N/2:N/2-1)/(N*dt)，uh = fftshift(fft(u))，反变换用ifft(ifftshift(uh))。这个组合确保频域数组从负频率到正频率连续排列，零频在正中间。确认方法是：不输入任何色散和非线性，直接把u0做一次正反 FFT，看能否原样返回。

5.5 加入损耗项后能量持续下降，但输出功率与光功率计对不上

现象：用公式计算 0.2 dB/km 损耗，跑了 20 km 应该掉 4 dB，模拟结果却掉了 8 dB。

原因：损耗项写法错了。方程里是∂u/∂z = -αu/2，所以离散迭代中每步乘的是exp(-αh/2)。如果误写成exp(-αh)，相当于每步算了两倍损耗，总损耗变成 2αL。

解决：记住损耗是对功率而言，幅度衰减系数是 α/2。更稳妥的做法是先忽略损耗跑通，确认波形正确后再加上，并且用一个固定 1 km 的测试距离对比手算的 dB 值。

6. 自适应步长与 GPU 化：一次仿真省掉三倍时间的习惯

固定步长实现起来简单，但计算效率并不好。脉冲在光纤里传播时峰值功率从 100 mW 慢慢下降，SPM 的影响随之变弱，非线性长度 LNL 越来越大，后半程根本不需要那么密的步长。我的习惯是改成自适应步长，用局部非线性相位增量来锁定 h，同时给一个上限防止步长超过色散尺度。

核心判据很简单：每一步走完，SPM 相位增量小于约定阈值，比如 0.05 rad。那h就取0.05 / (γ * max(|u|²))，并夹在[h_min, h_max]之间。代码骨架是：

%% 自适应步长定义 h_min = 1; % 最短 1 m，防止高功率时步长过小 h_max = 500; % 最长 500 m，防止超过色散尺度 alpha_step = 0.05; % 单步非线性相位阈值，rad % 每段循环开始前计算 peak_I = max(abs(u).^2); h_cur = alpha_step / (gamma * peak_I); h_cur = min(h_max, max(h_min, h_cur));

注意，h 变化时lin_op里的小 h 也要跟着变。所以色散算子不能放在循环外一次算完，而是每次更新h_cur后重新构造。虽然多了一次exp的向量运算，但比起整体步数减少带来的收益，这点开销完全可以忽略。记录输出距离时，也不要按固定索引取，改成“走到预定距离点就记录一次”。

这套自适应逻辑在 100 mW、20 km 的场景里，前 5 km 步长可能只有 30 m，后 10 km 能放到 500 m，总步数从 400 降到 150 左右，结果和固定 50 m 步长的频谱包络几乎完全重合。唯一的代价是没了均匀网格，分析瞬态时需要多写两行插值。

如果脉冲数不止一个，或者要做几百组参数扫描，还可以把主循环里的 FFT 全部扔到 GPU 上跑。MATLAB 里只需要把u初始化成gpuArray，lin_op也转成gpuArray，循环体一行都不用改。需要注意fftshift和ifftshift对 GPU 数组同样支持，但gather回 CPU 时要放在最终绘图之前。

这套流程走顺之后，我给自己定了一个死规矩：每次新建仿真脚本，第一件事是跑一小时前的固定步长脚本做基准，第二件事是关闭非线性验证解析展宽，第三件事才开 SPM。从那以后，我每次做传输仿真都会强制走一遍这三个环节，再让自适应步长去替我省时间。希望帮到你。

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