☰
光纤通信课程仿真:用Python实现分步傅里叶法求解非线性薛定谔方程
2026/10/15 15:52:24 网站建设 项目流程

简介:一套用于光纤通信课程教学与实验的MATLAB仿真代码包,面向通信工程专业学生、高校教师以及需要快速掌握光纤系统建模的初学者。资源以光发射-传输-接收为主线,系统覆盖光调制、光纤传播与色散、光电检测、误码率分析、前向纠错编码以及EDFA光放大器等关键主题,并包含PCM编码解码、NRZ/RZ码型时域频域波形、升余弦滚降系统、啁啾光纤光栅反射谱、布拉格光栅透射/反射谱、单模光纤色散等经典仿真场景,借助这些模块可将课堂理论转化为直观的波形曲线和性能指标,理解衰减、噪声、非线性效应及带宽受限对通信质量的影响。压缩包共含24个文件,核心为13个MATLAB的.m源文件,涵盖主程序、函数脚本与绘图辅助代码,另有光纤传光路径模拟、判决电路、非线性处理电路等辅助文件,资源总大小仅13KB,代码紧凑、模块独立,便于按需取用,适合课程作业、实验报告和自学复现。目前已有899人学习这份代码包。学习者可依据main入口脚本快速运行典型仿真,也可调整调制格式、放大器增益、光栅参数等变量,观察不同条件下的BER变化,为进一步研究波分复用、相干通信等方向提供可复用的实验基础。

1. 光纤通信课程仿真代码:为什么理论课上的公式要亲手跑一遍

我见过不少同学学光纤通信时卡在同一个地方:书上推导非线性薛定谔方程能推满三页纸,但问起「一个 1 mW 的高斯脉冲经过 20 km 标准单模光纤后变成什么样」,手里没有任何工具能回答。光纤通信课程仿真代码就是用来填这个空白的——把光信号在光纤中的衰减、色散、自相位调制这些抽象概念,变成你能亲眼看到、亲手调参的可运行程序。它能解决的不是“验证某个公式对不对”,而是帮你在做课程设计、期末项目或面试讲项目时,对“信号到底怎么变坏的、补偿为什么有效”有真实的体感。适合正在学光纤通信原理、做课程仿真作业、或者准备把光纤传输系统作为项目方向的同学。这一篇就沿着“选平台 → 建模型 → 写最小代码 → 跑典型场景 → 排坑”的顺序,把整套流程拆给你。

2. 仿真平台的选型与物理模型:先决定用什么写,再决定算什么

写光纤通信课程仿真,最难的不是写代码,而是在动手之前把平台和物理模型选对。选错了,后面每一步都在跟工具搏斗;模型简化过头,结果和教科书对不上,又会怀疑自己代码写错了。

2.1 Python 还是 MATLAB:课程仿真选型的三条标准

我自己的习惯是:先看三点,再决定平台。第一点是生态里有没有现成的 FFT、随机数和绘图工具,因为光纤仿真 90% 的运算都是傅里叶变换;第二点是交互式调试是否顺畅,能不能在跑完一次仿真后立刻改参数再看结果;第三点是代码能不能迁移到更大的场景,比如多信道 WDM 传输,而不是只能跑通一个演示。

在这三点上,Python 的 numpy、scipy、matplotlib 组合对课程仿真基本够用,而且免费、可迁移,适合写进课程报告甚至后续的毕业设计。MATLAB 胜在通信工具箱完备,眼图仪、调制解调函数开箱即用,但授权成本和代码可迁移性在团队协作里是麻烦。我的建议是:如果你只是单人完成课程仿真,用 Python 完全够;如果课程指定了 MATLAB 或者你需要和硬件平台联动,再切到 MATLAB。别在选平台上花太多时间,真正的耗时在物理模型。

2.2 物理模型拆解:从麦克斯韦方程组到非线性薛定谔方程

光纤通信里光信号的演化,严格讲要从麦克斯韦方程组出发,但课程仿真不需要这么重的起点。所有仿真代码共用的核心方程是非线性薛定谔方程,通常写作:

∂A/∂z = -i·β2/2·∂²A/∂T² - α/2·A + i·γ·|A|²·A

左边是光场包络 A 沿传播方向 z 的演化;右边第一项是群速度色散,β2 是色散参数,单位是 ps²/km;第二项是光纤损耗,α 是幅度损耗系数;第三项是自相位调制,γ 是非线性系数,|A|² 代表光功率。这个方程涵盖了课程里 80% 的仿真实验:不加第三项就是线性色散传输,去掉第一项只剩 SPM 频谱展宽,全部保留就是完整的光纤传输。

有个快速判断模型复杂度的经验:分别算出色散长度 LD = T0²/|β2| 和非线性长度 LNL = 1/(γ·P0),再和光纤长度 L 比较。如果 L 远小于 LD 和 LNL,衰减主导,用指数衰减公式就够了;如果 L 与 LD 同量级,色散主导,可以忽略非线性;如果 L 接近 LNL,SPM 必须保留。课程仿真里波分复用、四波混频这些课题才需要更完整的方程,做基础传输仿真时,上面的式子已经足够。

3. 用分步傅里叶法搭出第一个可运行仿真:最小代码与参数含义

方程有了,接下来就是数值求解。光纤通信课程仿真代码里最常用的求解方法是分步傅里叶法,几乎所有你能搜到的“光纤仿真源码”都在用它。把这部分理解透,后面换任何参数、加任何效应都是在它的框架上打补丁。

3.1 分步傅里叶法的核心思想:把方程拆成线性步和非线性步

非线性薛定谔方程的难点在于:色散项在频域里处理最方便,非线性项在时域里处理最方便,但两项混在一起时,没有一个域能同时让两项都简单。分步傅里叶法的做法是“耍赖”——把传输距离切成很多小段 h,在每一小段内,先假装只有色散和损耗,在频域里乘一个相位因子;再假装只有非线性,在时域里乘一个相位因子。这样交替进行,只要 h 足够小,误差就可控。

为了提高精度,常见做法是对称分步:每段先走半个线性步,再走一个非线性步,最后再走半个线性步。这样单步误差从 O(h²) 降到 O(h³),代码只多写一行,收益很划算。下面是完整的最小可运行例子,模拟一个高斯脉冲在 20 km 标准单模光纤里的传输。

3.2 最小仿真代码:一个高斯脉冲在光纤里的完整演化

import numpy as np from scipy.fft import fft, ifft import matplotlib.pyplot as plt # 物理参数:1550 nm 标准单模光纤的典型值 alpha_dB = 0.2 # 损耗,单位 dB/km alpha = alpha_dB / 4.343 # 换算成幅度损耗系数 /km beta2 = -20.0 # 色散参数 ps^2/km,反常色散区为负值 gamma = 1.3 # 非线性系数 1/(W·km) L = 20.0 # 光纤长度 km # 数值参数 T0 = 10.0 # 高斯脉冲 1/e 半宽,单位 ps P0 = 1.0e-3 # 峰值功率,单位 W(1 mW) N = 2**15 # 采样点数,2 的幂便于 FFT T = 600.0 # 时间窗口,单位 ps dt = T / N t = np.arange(-T/2, T/2, dt) omega = 2 * np.pi * np.fft.fftfreq(N, d=dt) # 角频率,单位 rad/ps # 输入高斯脉冲:注意幅值 = 根号功率 A = np.sqrt(P0) * np.exp(-0.5 * (t / T0)**2) # SSFM 参数 z_steps = 200 h = L / z_steps # 频域线性算子:色散 + 损耗合并 lin_half = np.exp((-1j * beta2 / 2 * omega**2 - alpha / 2) * h / 2) A_out = A.copy() for _ in range(z_steps): # 前半线性步,在频域完成 A_out = ifft(fft(A_out) * lin_half) # 非线性步,在时域完成 A_out = A_out * np.exp(1j * gamma * np.abs(A_out)**2 * h) # 后半线性步 A_out = ifft(fft(A_out) * lin_half) # 对照:同一段光纤,只走线性,不走非线性 A_lin = ifft(fft(A) * np.exp((-1j * beta2 / 2 * omega**2 - alpha / 2) * L)) plt.figure(figsize=(8, 3)) plt.plot(t, np.abs(A)**2 / P0, label='输入') plt.plot(t, np.abs(A_lin)**2 / P0, label='线性传输') plt.plot(t, np.abs(A_out)**2 / P0, label='SSFM 完整传输') plt.xlabel('时间 (ps)') plt.ylabel('归一化功率') plt.xlim(-100, 100) plt.legend() plt.grid(alpha=0.3) plt.show()

这段代码里最关键的是lin_half这个算子:它把色散项-i·β2/2·ω²和损耗项-α/2合并成频域乘法,每步只乘h/2,理由是对称分步需要前后各补一半。非线性步np.exp(1j * gamma * np.abs(A_out)**2 * h)直接在时域乘相位因子,功率越高、步长越大,相位旋转越快。你跑完应该看到:线性情况下脉冲只是对称展宽、峰值下降,而加上 SPM 后脉冲时域波形会略偏向一侧,频谱也明显展宽。

3.3 参数怎么定:步长、采样点数、时间窗口和参考功率

第一次跑通不算本事,知道每个参数为什么取这个值才算入门。我把最常调的几个参数整理成了一张表,课程仿真基本够用。

参数上面代码取值含义调参方向
N2**15采样点数太小则频谱混叠,建议不少于 2**14
T600 ps时间窗口要覆盖展宽后脉冲的 8 倍以上
z_steps200分段数非线性强时加大,步长 h 要小于 0.1/(γP0)
beta2-20 ps²/km色散参数符号不能错,负号对应反常色散
alpha0.2/4.343损耗系数dB 值除以 4.343 才是幅度系数
gamma1.3 /W/km非线性系数功率高时影响显著,弱光可忽略

时间窗口的坑最常见。上面代码里 T0=10 ps、beta2=-20 ps²/km 时,LD = T0²/|β2| = 5 km,传输 20 km 后脉冲展宽约 4 倍,窗口留到 600 ps 是安全的。如果你把 T 改成 60 ps 再跑,会看到波形从窗口边缘折回来,产生“伪脉冲”——这不是物理现象,是 FFT 天然周期性导致的。判断窗口够不够的土办法:先按 LD 估算展宽后的脉宽,再乘 8,取个整数当窗口。

步长同样有物理约束。非线性步的相位增量是 γ·P0·h,当这个值超过 0.1 rad 时,每步之间相互干涉会出现数值振荡。所以只要峰值功率高,就必须减小步长。课程场景里 P0 通常在 1 mW 量级,γP0≈0.0013 /km,200 步绰绰有余;如果你仿真几瓦的脉冲,z_steps 可能要加到几千。

4. 课程仿真常见场景:色散补偿、ASE 噪声与眼图测试

最小例子跑通后,就要往课程作业的要求靠拢了。我梳理了三个出现频率最高的课程仿真场景,每个场景对应一类报告图表,也对应一套容易抄错的细节。

4.1 场景一:色散补偿仿真,验证补偿光纤的效果

课程里最喜欢让你做的验证题,就是“一段 20 km 标准单模光纤后面接一段色散补偿光纤,证明脉冲宽度能被压回去”。原理很简单:传输链路累计色散是 β2·L,补偿光纤只要满足 β2_comp·L_comp + β2_fiber·L_fiber = 0 就能把频域相位赔回去。在仿真里,最简单的补偿就是在接收端对信号做一次反号色散运算,等效于理想 DCF。

# 接上节代码,A_out 是传输后的信号 # 理想补偿:频域乘一个反号色散相位 beta2_comp = -beta2 # 补偿光纤色散参数取反 L_comp = -beta2 * L / beta2_comp # 补偿长度,使累计色散为零 comp_op = np.exp((-1j * beta2_comp / 2 * omega**2 - alpha / 2) * L_comp) A_rec = ifft(fft(A_out) * comp_op) plt.figure(figsize=(8, 3)) plt.plot(t, np.abs(A)**2 / P0, label='原始脉冲') plt.plot(t, np.abs(A_rec)**2 / P0, label='补偿后') plt.xlim(-80, 80) plt.legend() plt.grid(alpha=0.3) plt.show()

这里有个课程报告里常被追问的点:为什么补偿后脉宽没有完全恢复?原因是链路里还有损耗和非线性,幅度变化不满足无损线性系统,色散补偿只能抵消相位,不能把损耗造成的能量损失补回来。如果想做更完整的补偿仿真,可以在链路中间加光放大器,两个仿真模块拼起来就是一个典型的“DCM+EDFA”链路模型。

4.2 场景二:ASE 噪声模型,看光放大器如何劣化信号

加了放大器之后,噪声就进来了。课程仿真里最常用的是 ASE 白噪声模型:掺铒光纤放大器增益为 G,自发辐射因子为 nsp,噪声功率谱密度近似为 2·nsp·h·ν·(G-1),其中 h 是普朗克常数,ν 是光频率。在离散仿真里加这类噪声的正确姿势是到频域加,因为 ASE 是加性白噪声,频域每个采样点的方差等于噪声谱密度乘以频谱分辨率。

# 放大器参数 G = 100.0 # 线性增益,对应 20 dB nsp = 2.0 # 自发辐射因子,典型值 1.5~3 h_planck = 6.626e-34 nu = 193.1e12 # 1550 nm 对应频率 # 在信号幅度上加入 ASE 噪声 df = 1.0 / (T * 1e-12) # 频率分辨率,单位 Hz S_ase = 2 * nsp * h_planck * nu * (G - 1) # 双偏振 ASE 谱密度 noise_power = S_ase * df # 每个频点上的噪声功率 noise = np.sqrt(noise_power / 2) * (np.random.randn(N) + 1j * np.random.randn(N)) A_ase = ifft(fft(A_out * np.sqrt(G)) + noise)

这段代码的坑在于单位换算:仿真里 t 的单位是 ps,而普朗克常数用的单位是 J·s,频率是 Hz,所以时间窗口必须从 ps 换算成 s,也就是再乘 1e-12,频谱分辨率才正确。噪声的实部和虚部各取一半功率,所以幅度要除以根号2,否则总噪声功率会翻倍。很多同学做“噪声系数越高 BER 越差”的曲线时发现趋势不对,基本都是这里差了 3 dB 或者少了归一化。

4.3 场景三:眼图与 Q 因子,把结果变成课程报告里的图

课程作业的“硬通货”是眼图和 BER 曲线。眼图的做法是:把接收信号按照符号周期切片,然后全部叠加画在一起,张开的“眼睛”越大说明信号质量越好。配套的量化指标是 Q 因子,它把眼图的上下电平均值和方差压缩成一个数,再换算成误码率。

# 假设 rec 是接收端采样后的 NRZ 信号,符号周期 Ts=50 ps Ts = 50.0 Ns = int(round(Ts / dt)) # 每个符号的采样点数 n_symbols = int(len(rec) // Ns) eye_segments = rec[:n_symbols * Ns].reshape(n_symbols, Ns) plt.figure(figsize=(6, 4)) for idx in range(min(n_symbols, 200)): plt.plot(eye_segments[idx], color='blue', alpha=0.03) plt.xlabel('符号内时间 (样本点)') plt.ylabel('幅度') plt.show() # Q 因子:对“1”电平和“0”电平分别统计 # 采样点取每个符号中间的 30% center = eye_segments[:, int(Ns*0.35):int(Ns*0.65)] levels = center.max(axis=0) # 简化:把中心区最大值当“1”电平分布 mu1, sigma1 = np.mean(levels), np.std(levels) mu0, sigma0 = 0.0, 0.05 # “0”电平按噪声底估算 Q = (mu1 - mu0) / (sigma1 + sigma0) BER = 0.5 * np.erfc(Q / np.sqrt(2))

注意,画眼图前最好加一个接收机低通滤波器,模拟光电探测器的带宽限制,否则眼图的过渡带会非常窄,张不开,看起来像一团乱麻。低通滤波可以直接在频域乘一个高斯形滤波器,带宽按 0.7~0.8 倍符号速率设置,这也是课程里容易漏掉的细节。

5. 光纤通信仿真避坑指南:从结果不对劲到参数玄学的 5 个排查方向

跑仿真最折磨人的是“代码没报错,但结果就是不对”。这类问题表面看是参数玄学,实际全是物理单位、符号约定和离散化细节。我把自己踩过和帮别人查过的五类高频问题列在这里,按“现象 → 原因 → 解决”写,你遇到类似情况可以直接对着查。

1. 波形从窗口边缘折回来,出现镜像脉冲。原因是时间窗口太小,展宽的脉冲碰到了窗口边界,FFT 的周期性把它当成新信号折回另一侧。解决方法是先估算展宽后的脉冲宽度,窗口至少给它留 8 倍余量,然后把 T 加大重跑。

2. 脉冲宽度不增反减,或者演化方向和教科书完全相反。多半是 β2 的符号写反了。1550 nm 的标准单模光纤处于反常色散区,β2 是负值;如果你按公式 β2 = -λ²/(2πc)·D 推算,D=17 ps/(nm·km) 算出来大约是 -21 ps²/km,符号别手滑。正负出错时,色散和非线性的竞争方向全反,结果自然对不上教科书。

3. 损耗衰减量明显不对,功率掉得比理论值快很多。alpha 的单位问题。0.2 dB/km 要先除以 4.343 换算成幅度系数,而不是直接拿来用;如果直接写alpha / 2,你把 dB 值当线性值用了。快速自检:20 km 光纤损耗应为 4 dB,功率衰减到 0.4 倍左右,如果代码里衰减成 0.01,一定是单位没换算。

4. 提高非线性强度后,波形出现锯齿状振荡,加密步长也没用。这是数值发散的特征。先检查 γ·P0·h 是否小于 0.1,如果峰值功率在传输过程中被 SPM 放大,局部功率可能远高于 P0,步长要按传输过程中的最大功率来算,而不是输入功率。建议在循环里加一个np.max(np.abs(A)**2)打印,观察功率突变的位置。

5. 计算出的 BER 好得离谱,或者加不加噪声结果没区别。大概率是 ASE 噪声功率没标定对,或者随机种子固定后恰好抽到一组好噪声。处理方法是:跑多次随机种子,取 BER 的平均值和方差;标定噪声时检查频率分辨率的单位,ps 域仿真一定要乘 1e-12 才能和 Hz 对接。

6. 让仿真代码变成课程项目里最耐用的部分:三个值得留下的习惯

走到这一步,你手里已经有能跑的最小 SSFM 代码和三个典型场景了。最后分享三个我后来做课程项目和毕业设计时觉得最值得留下的习惯,它们能帮你把“能跑的脚本”升级成“能反复使用的仿真工具”。

第一个习惯是把所有物理参数收进一个配置文件,而不是散落在代码里。我通常用 dataclass 把光纤参数、脉冲参数、数值参数装成三个对象,跑新场景时只改配置,不动核心传输函数。这样课程报告的“参数表”可以直接读出来,而且换光纤型号、换波长时,改动范围一目了然。第二个习惯是在核心传输函数里加一个“验证开关”:设 gamma=0 时,仿真结果要和纯色散理论解完全重合;设 beta2=0、弱光时,频谱展宽应该符合 SPM 的理论值。每次改完代码先跑这两个验证,通过再继续加功能,能省掉大半排查时间。第三个习惯是刻意把一个场景做成参数扫描,比如固定光纤长度,扫描输入功率从 0.1 mW 到 50 mW,画出输出脉宽随功率变化的曲线。这条曲线能同时讲出色散、SPM 和两者竞争三个知识点,比单张眼图更能展示你建模的能力。

我后来做课程项目时,习惯先写一个能打印所有参数和验证量的最小脚本,确认跑出来的第一个数是对的,再往里面加模块、加图表。这比一次性写完再回头找 bug 省力得多。希望帮到你。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询