★ FEATURED ARTICLE 简介面向光纤通信与MATLAB仿真学习者的色散研究代码包适合具备通信原理基础、希望用脚本复现光纤模式分布与色散补偿思路的读者。压缩包共14个文件均为.m脚本大小约14KB主要包含LP模式系列函数与三维极化绘图函数可用于搭建轻量级光纤信道模型、观察模式场分布及信号演化过程。已有2858人学习下载在课程实验和入门仿真中具备一定参考价值。读者可借助这些脚本理解模式色散、材料色散与波导色散的基本表现并进一步结合频域重叠保留补偿方法、衰减因子以及四波混频、自相位调制、交叉相位调制等非线性项比较补偿前后的信号质量差异量化评估不同策略的效果。整体结构精简便于按函数逐个调试、替换参数和扩展模块可作为光纤通信色散分析与Simulink建模之前的代码基础。1. 光纤通信链路里的色散预算先在 MATLAB 里算明白一条 10 Gbps 的链路发射端眼图睁得挺开传到 80 km 外的接收机就糊成一团误码率卡在 1e-6 下不去。现场第一反应往往是查光功率、换连接器、怀疑激光器其实很多时候问题出在光纤通信里最朴素的敌人身上——色散不同波长的光在光纤里跑的速度不一样脉冲还没衰减完就先被摊平了。这份基于 MATLAB 的色散研究压缩包干的就是把这件事从“凭经验猜”变成“先在电脑上算一遍”给定波长、色散系数、光纤长度和脉宽跑出脉冲展宽曲线、眼图和 Q 值让链路预算在开工前就有数。适合做光模块验证、链路预算、课程设计以及想弄明白 D、β2、色散长度之间换算关系的开发者。2. 环境与目录压缩包跑起来前先补齐的三件事拿到一个 MATLAB 工程压缩包最怕的是双击 main.m 之后满屏红字报的还是Undefined function。这多半不是代码坏了而是路径、版本、参数三者里有一个没对齐。先把运行环境摸清楚后面调参才不会被环境问题带偏。2.1 MATLAB 版本与工具箱边界核心算法只用到fft、ifft、fftshift、矩阵逐元素运算和plot这些从 R2016b 起就稳定存在所以版本门槛并不高。真正需要注意的是语法兼容性如果脚本里用了arguments参数块R2019b 之后才支持、字符串双引号R2017a 之后或者隐式扩展之外的新函数在老版本上会直接报解析错误。我一般建议用 R2020a 以上的版本跑兼容性和fft性能都更稳。工具箱方面Signal Processing Toolbox 主要用于可选的滤波器设计和眼图辅助函数Communications Toolbox 里的comm.EyeDiagram能画更漂亮的眼图但这两个都不是必需品。项目里如果没有comm.EyeDiagram用普通的plot叠加比特周期同样能看眼图只是样式朴素一点。检查方式很简单在命令行敲ver看输出里有没有 Signal Processing Toolbox 和 Communications Toolbox。没有也能跑但如果你打开eye_q.m发现第一行调用了comm.EyeDiagram要么装工具箱要么把那几行换成手写叠加逻辑后面会讲。还有一个容易忽略的点路径。MATLAB 不会自动递归加载子目录main.m里通常会有addpath(genpath(src))这类语句。如果你直接cd到src里运行某个函数config目录下的参数文件就找不到了报错信息会指向一个看起来莫名其妙的变量。运行前先确认当前目录是工程根目录这是血泪经验。2.2 目录结构与脚本职责这类色散仿真工程结构一般不会太复杂压缩包解开后大致是这么一套职责划分fiber_dispersion/ ├── main.m % 入口读参数、跑仿真、出图、存结果 ├── config/ │ └── fiber_params.m % 所有可调参数集中在这里 ├── src/ │ ├── D_S_to_beta.m % D、S 与 beta2、beta3 的单位换算 │ ├── gaussian_pulse.m % 生成高斯或超高斯输入脉冲 │ ├── prbs_seq.m % 伪随机比特序列用于眼图 │ ├── ssfm_propagate.m % 分步傅里叶法求解脉冲演化 │ └── eye_q.m % 眼图叠加与 Q 值估计 ├── results/ % 输出的 png、mat 文件落盘位置 └── README.md如果你手上的包文件名略有出入按职责定位就行凡是做单位换算的一定在D_S_to_beta这一类文件里凡是出现for k 1:nz循环、循环体里有fft的就是分步傅里叶法的主循环凡是出现reshape和plot叠加的多半是眼图。不要一上来就逐行读先按职责把文件分类再挑核心的三个看。config/fiber_params.m是整个工程的“控制面板”所有魔法数字都应该在这里而不是散落在各个函数里。判断一个工程是否好改就看这一点如果改一个波长要翻五个文件那二次开发成本会很高。2.3 一次基准运行的完整命令最省事的跑法是在 MATLAB 命令行里直接执行入口脚本cd(fiber_dispersion) % 进到工程根目录 main % 跑基准算例如果要在无图形界面的服务器上批量跑用批处理模式matlab -batch cd(/path/to/fiber_dispersion); main-batch会在执行完后自动退出适合挂到后台或者写进 shell 脚本里循环扫参。注意-batch模式下figure不会弹窗saveas和print仍然可用所以main.m里如果有figure却没配saveas跑完你什么都看不到只剩一个空的results目录。我第一次在服务器上跑就栽在这以为程序没执行其实是图被吞了。基准运行的预期输出是三样东西一张脉冲时域波形随距离演化的图、一张眼图、一个命令行打印的 Q 值dB。如果只出了前两张、Q 值是NaN八成是眼图采样点索引越界或者比特序列长度不够一个完整周期后面避坑章节会展开。参数文件的内容大致长这样每一项的含义值得先过一遍function p fiber_params() % 光纤与信号参数集中定义单位写在注释里 p.lambda 1550; % 波长nm p.D 17; % 色散系数ps/(nm·km)G.652 在 1550 nm 的典型值 p.S 0.057; % 色散斜率ps/(nm^2·km) p.alpha_dB 0.2; % 衰减系数dB/km p.gamma 1.3; % 非线性系数1/(W·km) p.L 80; % 光纤长度km p.P0_dBm 0; % 入纤功率dBm p.T0 25; % 高斯脉冲 1/e 半宽ps p.bitrate 10e9; % 比特率bps p.nz 200; % 分步傅里叶的步数 p.N 2^14; % 时间采样点数 p.tspan 3200; % 时间窗总长ps endD和S是工程单位T0是脉冲宽度nz和N是数值精度旋钮。N取 2 的幂是为了 FFT 效率tspan要覆盖脉冲展宽后的全部能量否则循环卷积会把尾巴绕回头部出现看起来像“反射”的假波形。这几个参数的互相牵制关系是后面调参的主要战场。3. 核心模块拆解从 D 到 β2再到脉冲在光纤里怎么走整个工程的计算链路可以概括成三步把工程单位换成非线性薛定谔方程能吃的单位、用分步傅里叶法把脉冲一段一段推过去、把输出波形切成比特周期叠成眼图。三步里任何一步的单位或栅格错位最后的曲线都会“看起来像那么回事但数值全错”这也是色散仿真最容易翻车的地方。3.1 色散系数换算D、S 与 β2、β3 的三角关系工程上习惯用色散系数D单位是 ps/(nm·km)意思是每公里、每纳米波长差带来多少皮秒的时延差。但非线性薛定谔方程里用的是β2单位 ps²/km描述的是群速度随频率的变化。两者描述的是同一件事换算公式是D -(2πc / λ²) · β2负号来自约定β2 0叫正常色散区D 0β2 0叫反常色散区D 0。1550 nm 处 G.652 光纤的D ≈ 17对应β2 ≈ -21.7 ps²/km正好落在反常色散区这也是为什么长距离链路要用色散补偿。代码实现时单位是最容易错的地方。D的 SI 量纲是 s/m²换算系数是 1e-6S的 SI 量纲是 s/m³系数是 1e3。这两个系数不一样混用就会差三个数量级function [beta2, beta3] D_S_to_beta(D, S, lambda_nm) % D : 色散系数ps/(nm·km) % S : 色散斜率ps/(nm^2·km) % lambda_nm: 波长nm % beta2 : 返回 ps^2/km % beta3 : 返回 ps^3/km c 299792458; % 光速m/s lam lambda_nm * 1e-9; % nm - m D_si D * 1e-6; % ps/(nm·km) - s/m^2 S_si S * 1e3; % ps/(nm^2·km) - s/m^3 beta2_si -D_si * lam^2 / (2*pi*c); % s^2/m % 由 S (2*pi*c)^2*beta3/lam^4 4*pi*c*beta2/lam^3 反解 beta3 beta3_si (S_si - 4*pi*c*beta2_si/lam^3) * lam^4 / (2*pi*c)^2; beta2 beta2_si * 1e27; % s^2/m - ps^2/km beta3 beta3_si * 1e39; % s^3/m - ps^3/km end这段代码的关键在最后两行的换算系数。1 ps²/km 1e-27 s²/m所以从 SI 换回工程单位要乘 1e271 ps³/km 1e-39 s³/m对应乘 1e39。这两个数字记不住没关系但一定要在代码里写清楚注释否则半年后回来看自己都不确定是乘还是除。拿基准参数验算一下D 17、λ 1550 nm算出来beta2 ≈ -21.7 ps²/kmbeta3 ≈ 0.13 ps³/km。如果算出来是-2.17e-26这种量级说明忘了乘 1e27单位还停在 SI 制如果符号是正的说明负号约定反了。这两个错误都会让后面的脉冲展宽方向完全反过来补偿模块会越补越差。β3的作用是描述色散斜率在窄带、单波长场景里影响很小但在宽带波分复用或者超短脉冲场景下不能忽略。基准算例里它主要用来验证S参数有没有被正确读入真正影响脉冲形状的还是β2。3.2 分步傅里叶法把 NLSE 拆成频域线性步和时域非线性步脉冲在光纤里的演化由非线性薛定谔方程描述∂A/∂z -α/2·A - i·β2/2·∂²A/∂T² i·γ·|A|²·A右边三项分别是衰减、色散、非线性。这个方程没有通用的解析解数值上最常用的是分步傅里叶法SSFM把光纤切成nz段每段里先假设只有色散和衰减频域里是乘一个相位因子再假设只有非线性时域里是乘一个相位因子。因为色散在频域是简单的二次相位非线性在时域是简单的强度相关相位两边交替做就能逼近真实解。function [A_out, omega] ssfm_propagate(A_in, t, L, beta2, gamma, alpha_dB, nz) % A_in : 输入复包络列向量单位 sqrt(W) % t : 时间栅格ps % L : 光纤长度km % beta2 : ps^2/km % gamma : 1/(W·km) % alpha_dB: dB/km % nz : 步数 dt t(2) - t(1); N numel(t); % 频域栅格单位 rad/ps用 fftshift 把零频挪到中间 omega 2*pi*fftshift((-(N/2):(N/2-1)) / (N*dt)); h L / nz; % 每步长度km alpha alpha_dB / (10*log10(exp(1))); % dB/km - 1/km功率衰减系数 A A_in; for k 1:nz % 线性半步频域乘色散与衰减相位 A_f fftshift(fft(ifftshift(A))); A_f A_f .* exp((-alpha/2 - 1i*beta2/2*omega.^2) * h); A fftshift(ifft(ifftshift(A_f))); % 非线性整步时域乘强度相关相位 A A .* exp(1i*gamma*abs(A).^2*h); % 线性半步与第一步对称 A_f fftshift(fft(ifftshift(A))); A_f A_f .* exp((-alpha/2 - 1i*beta2/2*omega.^2) * h); A fftshift(ifft(ifftshift(A_f))); end A_out A; end逐段说明。omega的计算里(-(N/2):(N/2-1))/(N*dt)给出的是以 Hz 为单位的频率栅格乘2*pi转成角频率 rad/ps。beta2 * omega.^2的量纲是ps²/km × 1/ps² 1/km和h相乘后无量纲指数项才合法。如果omega忘了乘2*pi色散量会差(2π)² ≈ 39.5倍脉冲展宽速度看起来“慢了很多”很容易被误判成参数填错。alpha的换算用10*log10(exp(1)) ≈ 4.343因为功率随距离按exp(-alpha*z)衰减而工程上的 dB 是10*log10。写成alpha_dB / 4.343和上面的表达式等价看个人习惯。混淆 dB/km 和 1/km 是最常见的错误之一0.2 dB/km 对应约 0.046 1/km差 4.3 倍80 km 下来功率差几十倍眼图直接闭合。对称分步半步线性、整步非线性、半步线性比非对称分步精度高一阶代价是每步多做两次 FFT。nz 200对 80 km 来说步长 0.4 km在 10 Gbps、几毫瓦功率下足够收敛。判断步长够不够不要凭感觉把nz翻倍再跑一遍如果 Q 值变化小于 0.1 dB就认为收敛了如果差很多说明非线性相位单步变化超过 0.05 rad得继续加步数。输入脉冲的生成一般用高斯function A gaussian_pulse(t, T0, P0_dBm) % t : 时间栅格ps % T0 : 1/e 半宽ps % P0_dBm: 峰值功率dBm P0 1e-3 * 10^(P0_dBm/10); % dBm - W A sqrt(P0) * exp(-t.^2 / (2*T0^2)); endT0 25 ps对应 10 Gbps 下大约 0.25 个比特周期是个比较典型的窄脉冲设置。P0_dBm 0对应 1 mWgamma 1.3时非线性长度L_NL 1/(gamma*P0) ≈ 770 km远大于 80 km所以基准算例里非线性很弱脉冲展宽主要由色散贡献。这个判断很重要如果L_NL和L同量级就不能用纯色散的直觉去解释结果得把自相位调制一起考虑进去。3.3 眼图与 Q 值怎么把一条曲线变成一个可判断的指标单条脉冲波形只能看展宽看不出对通信系统的影响。要判断“这个链路能不能用”得把连续的比特流叠成眼图再看眼高和眼宽。工程上的做法是生成一段 PRBS 比特序列每个比特按 NRZ 或 RZ 调制到光载波上过光纤然后在接收端按比特周期切段叠加。function bits prbs_seq(order, nbits) % 7 阶 PRBS本原多项式 x^7 x^6 1 len 2^order - 1; reg ones(1, order); bits zeros(1, nbits); for k 1:nbits bits(k) reg(end); fb xor(reg(end), reg(end-1)); reg [fb, reg(1:end-1)]; end bits bits(1:min(nbits, len)); endPRBS 的作用是保证比特序列里有足够多的 0/1 跳变和长连码眼图才有代表性。如果只用[1 0 1 0 ...]这种周期序列眼图会假性张开Q 值虚高看不出真实的码间干扰。Q 值的估计有两种路子。严谨的做法是统计眼图中心采样点的均值与标准差按Q (μ1 - μ0)/(σ1 σ0)算再转 dB简化做法是只取“1”电平的均值和标准差Q 20*log10(μ/σ)。后者的绝对值不能直接对标误码率门限但用来做相对比较、看趋势足够了。基准脚本一般用简化版速度快扫参时不容易把时间耗在统计上。function Q_dB eye_q(A, t, bitrate, nbits) % A : 接收端复包络 % t : 时间栅格ps % bitrate: 比特率bps % nbits : 参与统计的比特数 Tb 1e12 / bitrate; % 比特周期ps dt t(2) - t(1); idx round((0:nbits-1)*Tb/dt) 1; idx idx(idx 1 idx numel(A)); y abs(A(idx)).^2; % 采样点功率 mu mean(y); sig std(y); Q_dB 20*log10(mu / max(sig, eps)); end这段代码里max(sig, eps)是防止标准差为零导致除零在理想无噪声仿真里很常见。idx的边界检查也不能省时间窗如果没覆盖完整的nbits个比特索引会越界MATLAB 直接报错退出。更稳妥的做法是让nbits * Tb tspan在参数文件里就约束好。眼图叠加本身很简单把A按Tb切段、每段归一化时间轴、全部画在同一张图上透明度调低一点就能看出“眼睛”的张开程度。真正需要留意的是采样相位如果idx固定取比特起点而实际最佳采样点在比特中心Q 值会偏低看起来像链路性能差其实是采样位置没对齐。我一般会在一个比特周期内扫一遍采样相位取 Q 值最大的那个点作为报告值这样不同配置之间的比较才公平。4. 参数标定与实验设计色散长度、补偿比和扫参矩阵仿真跑通只是第一步真正决定这份工程有没有用的是参数怎么取。色散研究里最核心的两个量是色散长度L_D和色散补偿比前者决定“多远开始明显展宽”后者决定“补偿光纤要拉多长”。这两个数算对了后面的扫参才有意义。4.1 先算色散长度再决定光纤长度怎么取色散长度定义为L_D T0² / |β2|它描述的是脉冲宽度因为色散展宽到原来的 √2 倍所需要的距离。T0 25 ps、β2 21.7 ps²/km时L_D ≈ 28.8 km。这意味着 80 km 链路对应L/L_D ≈ 2.8脉冲会明显展宽眼图开始闭合正好是一个能看出问题的算例。如果只跑 10 kmL/L_D ≈ 0.35展宽不到 6%图上看不出差别容易误以为“色散影响不大”。function L_D dispersion_length(T0, beta2) % T0 : 脉冲 1/e 半宽ps % beta2: ps^2/km L_D T0^2 / abs(beta2); % km end这个函数虽然只有一行但建议单独成文件。原因是它会被多处调用主脚本用它决定光纤长度扫描范围补偿模块用它估算补偿后的残余色散画图时用它做归一化横轴。散落在各处手算改一个参数就要改五个地方早晚出错。扫光纤长度时合理的范围是0.2*L_D到5*L_D覆盖“几乎没展宽”到“完全糊掉”的全过程。低于 0.2 倍看不出趋势高于 5 倍脉冲已经展宽到时间窗边缘循环卷积开始污染结果。时间窗tspan也要跟着放大经验值是至少覆盖10*T0*sqrt(1(L/L_D)^2)否则展宽后的脉冲尾部会被截断能量不守恒Q 值算出来偏乐观。4.2 色散补偿比DCF 长度不是拍脑袋长距离链路要在接收端或中继处加色散补偿光纤DCF用负色散抵消正色散。补偿比的定义是CR (|D_DCF| · L_DCF) / (D_SMF · L_SMF)CR 1表示完全补偿残余色散为零CR 1是欠补偿CR 1是过补偿。实际工程里很少追求严格的 1因为完全补偿会让非线性效应尤其是四波混频在零色散点附近累积通常留 5% 到 10% 的残余色散更稳。按基准参数算SMF 段D 17、L 80 km总色散量1360 ps/nmDCF 取D -100 ps/(nm·km)完全补偿需要L_DCF 1360/100 13.6 km。这个长度不算短实际模块里会用更高负色散的 DCF 把长度压到几公里但仿真里用 -100 这个量级足够说明问题。function L_dcf dcf_length(D_smf, L_smf, D_dcf, CR) % CR: 目标补偿比1 表示完全补偿 L_dcf CR * D_smf * L_smf / abs(D_dcf); end补偿后的残余色散是D_smf*L_smf - |D_dcf|*L_dcf把它代回D_S_to_beta得到等效β2再跑一遍 SSFM就能看到眼图重新张开。这里有个容易翻车的点补偿光纤的β3符号和 SMF 相同还是相反取决于具体光纤设计仿真里如果直接把β2取反而不改β3宽带场景下补偿后的残余色散斜率会对不上眼图只在中心波长处好看。窄带单波长场景可以忽略波分复用场景必须把S一起补偿。4.3 扫参矩阵把单点仿真变成一张等高线图单点仿真只能回答“这组参数行不行”扫参才能回答“哪个区间行”。把D和L各取一组值双重循环跑 SSFM把 Q 值存成矩阵用imagesc画出来一眼就能看出可用区间。D_list 15:1:19; % ps/(nm·km) L_list 40:20:200; % km Qmap zeros(numel(D_list), numel(L_list)); for i 1:numel(D_list) for j 1:numel(L_list) p fiber_params(); p.D D_list(i); p.L L_list(j); [beta2, ~] D_S_to_beta(p.D, p.S, p.lambda); % ... 生成脉冲、过光纤、算 Q Qmap(i, j) Q_dB; end end imagesc(L_list, D_list, Qmap); colorbar; xlabel(Length (km)); ylabel(D (ps/nm/km));双重循环里每次都要重新生成脉冲和 PRBS开销不小。优化办法是把 PRBS 和输入脉冲提到循环外只把beta2和L传进 SSFM能省掉一半以上的时间。nz也要跟着L缩放保持步长恒定否则长距离那几组会因为步长过大而精度不够Q 值出现非物理的抖动。扫参结果一般会呈现出一个“高 Q 值走廊”在某个D·L乘积附近 Q 值最高偏离越远越低。这个走廊的宽度就是色散容限10 Gbps NRZ 大概在 ±800 ps/nm 量级和教科书上的经验值能对上说明仿真链路是可信的。如果扫出来的走廊位置明显偏移先别怀疑算法回去查D_S_to_beta的单位换算和alpha的 dB/km 转换这两个地方错了整张图会整体平移。5. 避坑与排查从量纲错位到 FFT 栅格混叠的五种翻车现场色散仿真的结果有一种很强的迷惑性曲线画出来总是“像那么回事”指数衰减、脉冲展宽、眼图闭合视觉上都合理但数值可能差了十几倍。下面这几条是调试时最常撞上的每条都按“现象、原因、解决”走一遍。5.1 脉冲展宽曲线指数爆炸或出现高频振荡现象随着L增大脉冲峰值不是平滑下降而是在某一段突然剧烈振荡甚至出现负值功率把nz调大后振荡位置会移动。原因FFT 栅格和时间窗不匹配导致频域混叠。时间窗tspan太短展宽后的脉冲尾部绕回到头部或者omega栅格没有用fftshift对齐线性相位作用到了错误的频率分量上。解决先按tspan ≥ 10*T0*sqrt(1(L/L_D)^2)检查时间窗不够就加倍N或tspan。然后确认fft和ifftshift的配对方式MATLAB 的fft输出零频在第一个点如果omega用fftshift排过序变换前后必须用ifftshift/fftshift成对处理少一个就会让频谱整体平移半个带宽表现为高频振荡。最后做步长收敛性检查nz翻倍结果不变才算过。5.2 加了色散补偿之后眼图反而更差现象补偿前 Q 值 12 dB按CR 1加了 DCF 之后掉到 8 dB检查补偿长度也没算错。原因β2的符号约定搞反了。D 17对应β2 -21.7是负值如果补偿模块里直接把D取反后套同一个换算公式得到的β2符号会错相当于把色散量加倍而不是抵消。解决不要对D取反再换算而是直接对β2取反或者在D_S_to_beta外面显式乘-1并加注释。验证方法很简单补偿后算残余β2完全补偿时应该接近 0如果算出来是-43.4说明符号反了。再退一步用脉冲前沿判断反常色散区β2 0里高频分量跑得慢脉冲波形会向前沿倾斜补偿正确时倾斜方向应该反转。5.3 Q 值恒为极大值或Inf现象不管光纤多长、色散多大Q 值都在 40 dB 以上眼图看起来完美张开。原因比特序列太短或者周期性太强采样点功率的标准差接近零。比如只用了[1 0 1 0]循环每个比特的采样值几乎相同std算出来是 1e-15 量级除下来 Q 值爆表。另一种情况是nbits设得比实际生成的比特数大idx越界后被截断统计样本只剩几个点。解决换成 7 阶以上的 PRBS长度至少 127 比特在eye_q里加一行断言检查numel(y) 0.8*nbits不够就报错而不是静默截断。还可以人为加一点ASE噪声让 Q 值落在合理区间10 到 20 dB这样不同配置之间的差异才看得出来。理想无噪声仿真里 Q 值本来就该是无穷大非要拿它做比较得先把噪声模型补上。5.4 功率随距离掉得离谱80 km 后几乎为零现象alpha_dB 0.2按理 80 km 衰减 16 dB功率剩 1/40实际跑出来功率只剩 1e-5差了三个数量级。原因dB/km 到 1/km 的换算用错了系数。把alpha_dB直接当成了 1/km 用等于把衰减系数放大了 4.3 倍或者把功率衰减和场衰减搞混exp(-alpha*z)里少除了 2。解决统一用alpha alpha_dB / (10*log10(exp(1)))并在代码注释里写明这是功率衰减系数。验证方法无非线性、无色散时输出功率应该等于P0*exp(-alpha*L)手算对比一下16 dB 衰减对应功率比 0.025对不上就是换算错了。这个错误不会让曲线形状变得离谱只会整体压低所以特别容易被忽略。5.5 换一个波长所有结果都偏移现象1550 nm 下调好的参数换到 1310 nm 后 Q 值、展宽曲线全部对不上但D和S看起来都按手册改了。原因β2和β3不是独立参数它们随波长变化如果只改了D而没同步改S或者D_S_to_beta里的lambda_nm忘了传新值换算出来的β2还是旧波长的。解决把波长作为参数结构体的一部分所有换算函数都从p.lambda取值不允许在函数内部写死 1550。换波长时同时更新D和SG.652 在 1310 nm 附近D ≈ 0、S ≈ 0.09和 1550 nm 完全不同。跑完之后用零色散波长做一次自检D接近 0 时β2也应该接近 0脉冲几乎不展宽如果还展宽得厉害说明换算链路里有一个环节没跟着波长走。6. 二次开发把单脚本改造成可复用的色散分析工具链基准脚本跑通之后大多数人下一步都会面临同一个问题每换一组参数就要改一次main.m改完还得手动记录结果跑上十几组就乱了。我一般会做三件事把它变成工具链参数结构体化、结果落盘、加一组能自动跑的校验。第一件事是把main.m拆成“配置、计算、出图”三层。配置层只留fiber_params计算层封装成run_case(p)输入一个参数结构体、输出一个结果结构体Q 值、展宽比、残余色散出图层单独一个函数。这样扫参时只需要在循环里构造不同的p调用run_case结果自动进数组不用碰任何画图代码。function r run_case(p) % 输入参数结构体输出结果结构体不画图、不落盘 [beta2, beta3] D_S_to_beta(p.D, p.S, p.lambda); A_in gaussian_pulse(p.t, p.T0, p.P0_dBm); A_out ssfm_propagate(A_in, p.t, p.L, beta2, p.gamma, p.alpha_dB, p.nz); r.Q_dB eye_q(A_out, p.t, p.bitrate, p.nbits); r.beta2 beta2; r.L_D p.T0^2 / abs(beta2); r.P_out sum(abs(A_out).^2) * (p.t(2)-p.t(1)); end第二件事是结果落盘。每次run_case返回的r加上参数p一起存成.mat文件名带上关键参数比如D17_L80_T025.mat。这样跑完几十组之后可以直接load回来画对比图不用重跑。文件名里带参数看着土但比在 MATLAB 里维护一张索引表可靠得多尤其是在多台机器上跑的时候。第三件事是加校验这也是最容易省掉、但后面最后悔的一步。色散仿真至少要有两个自检能量守恒和解析解对照。无衰减alpha_dB 0、无非线性gamma 0时SSFM 输出的总能量应该和输入一致误差在 1e-10 量级纯色散情况下高斯脉冲的展宽满足T1 T0*sqrt(1(L/L_D)^2)把仿真得到的脉宽和这个公式对一下偏差超过 5% 就说明步长或时间窗有问题。% 能量守恒自检 E_in sum(abs(A_in).^2); E_out sum(abs(A_out).^2); assert(abs(E_out/E_in - 1) 1e-6, 能量不守恒检查 alpha 与 FFT 栅格); % 解析解对照纯色散、无衰减 T1_sim sqrt(sum(p.t.^2 .* abs(A_out).^2) / sum(abs(A_out).^2)); T1_ref p.T0 * sqrt(1 (p.L / r.L_D)^2); assert(abs(T1_sim/T1_ref - 1) 0.05, 展宽与解析解偏差过大);这两个断言写完后面再改beta2换算、换波长、加补偿模块跑一遍就知道有没有碰坏东西。比起肉眼盯曲线断言的反馈快得多也不会因为“图看起来差不多”而放过一个数量级的错误。从那以后我每次动色散换算或者 FFT 栅格都强制先跑一遍能量守恒和解析解对照两个都过了才看眼图。这个习惯救过我好几次——有一次改omega的排序方式图上面看着一切正常能量却掉了 3%追下去发现是ifftshift少写了一处。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
