简介面向光纤通信课程学习与实验教学的MATLAB仿真代码包围绕光信号产生、光调制、光纤传播、光电检测与误码率分析等核心环节展开借助数值仿真帮助本科生及研究生理解光调制、光纤衰减、色散、信道噪声等抽象概念适合课程设计、实验辅助与科研入门。压缩包共24个文件以13个.m脚本为主体涵盖PCM编解码、光纤传光路径模拟、TE01/HE11等模式场形、单模光纤色散、啁啾光纤光栅反射谱、布拉格光栅反射/透射谱、NRZ/RZ码时频域波形、升余弦滚降系统等模块并附有非线性处理电路、D触发器定时判决电路、幅度判决电路等辅助模型整体体积仅13KB轻量易用。目前已有899人学习下载。代码模块划分清晰、便于直接运行与参数修改既能辅助掌握脉冲编码调制、光纤模式与光栅滤波器设计等原理也可作为进一步开展误码率分析、信道编码仿真和传输性能优化的二次开发基础。1. 光纤通信课程仿真代码为什么能跑通和能讲清是两回事光纤通信课程里真正卡住学生的往往不是傅里叶变换而是“公式我都懂但让我把 OOK 和 QPSK 的误码率曲线画出来我连从哪下手都不知道”。这门课的仿真代码核心诉求从来不是写出一个能出图的脚本而是把发射机、信道、接收机三段链路拆成可独立验证的模块每一段都对应一个物理过程每一行代码都能映射回课本公式。你需要的是这样一套代码跑通只是底线能改参数、能复现理论曲线、能定位“为什么我的 BER 比理论值差一个数量级”才是价值所在。本文从零搭一套以 Python 为主的光纤通信仿真框架覆盖从光源调制、光纤传输到相干接收与误码率统计的完整闭环并指出那些课本不会写、但仿真里一定会踩的坑。2. 仿真框架选型为什么选择 Python 而非 MATLAB以及最小可运行架构2.1 选型理由MATLAB 工具箱虽全但 Python 更适合课程迭代与二次开发常见做法是 MATLAB 加通信工具箱封装度高但黑匣子问题严重。课程仿真恰恰需要你敢打开每一层看数据长什么样。Python 的 numpy 和 scipy 足以覆盖光纤通信仿真 90% 的数值计算需求而且光通信领域开源生态正快速补齐从信号生成到信道建模都有可参考的实现。Python 版本迭代快调试手段丰富出了问题能直接打印中间变量这对课程学习来讲比工具箱的“一键出图”重要得多。若你所在课程强制用 MATLAB本文的模块划分思路同样可直接迁移区别只在函数名和绘图语法。2.2 最小可运行架构发射机-信道-接收机三段式整个仿真链路划分为三个相互独立的模块。发射机负责把比特流映射为光信号波形信道负责施加衰减、色散与噪声接收机负责将光信号还原为比特并统计误码。三段之间通过“波形数组 采样率 符号率”这三个约定接口衔接这样你可以单独替换任何一段而不影响其余部分。以下是最小骨架跑通后你的代码就有了可扩展的基地import numpy as np # 全局参数所有模块共享 fs 80e9 # 采样率 80 GHz R 10e9 # 符号率 10 GBaud sps int(fs / R) # 每符号采样点数 8 # 随机比特流 bits np.random.randint(0, 2, 10000) # 发射机OOK 映射 symbols bits # OOK 下比特直接作为幅度符号 # 信道先占位后续填充色散与噪声 rx_signal symbols.copy() # 接收机硬判决 rx_bits (rx_signal 0.5).astype(int) # 误码统计 ber np.mean(rx_bits ! bits) print(f理想信道 BER {ber:.2e})这段代码的逻辑很直接参数区定义采样率、符号率和每符号采样点这是全仿真最核心的三个数字三者关系必须满足 fs R * sps。发射机部分是 OOK 的最简实现比特 0 和 1 直接映射为幅度 0 和 1。信道暂时直通接收机按 0.5 阈值做硬判决。这个骨架能跑但它不包含脉冲成形、光纤色散、噪声和匹配滤波距离“光纤通信仿真”还有距离下一步逐层填充。3. 光发射机仿真脉冲成形、马赫-曾德尔调制与 IQ 调制器建模3.1 为什么不能直接发送矩形脉冲脉冲成形是仿真第一课光纤带宽有限矩形脉冲在频域是 sinc 函数旁瓣会干扰相邻符号这就是码间串扰的根源。课程仿真必须加入脉冲成形滤波最常用的是升余弦和根升余弦滤波器。根升余弦在发射端和接收端各用一次合成升余弦响应匹配滤波后信噪比最优。滚降因子决定带宽占用与峰值功率的折中0.2 到 0.5 是工程常用区间滚降越小频谱效率越高但对时钟相位越敏感。以下是发射机完整建模包含比特映射、脉冲成形和 IQ 调制:from scipy.signal import upfirdn def rrc_filter(sps, beta, span11): 生成根升余弦滤波器系数 n np.arange(-span*sps//2, span*sps//2 1) / sps h np.zeros_like(n) for i, t in enumerate(n): if t 0: h[i] 1 - beta 4*beta/np.pi elif abs(t) 1/(4*beta): h[i] beta/np.sqrt(2) * ((12/np.pi)*np.sin(np.pi/(4*beta)) (1-2/np.pi)*np.cos(np.pi/(4*beta))) else: num np.sin(np.pi*t*(1-beta)) 4*beta*t*np.cos(np.pi*t*(1beta)) den np.pi*t*(1-(4*beta*t)**2) h[i] num/den return h / np.sqrt(np.sum(h**2)) beta 0.3 rrc_taps rrc_filter(sps, beta) # QPSK 映射每 2 比特映射一个符号 bits_qpsk np.random.randint(0, 2, 20000) symbols_qpsk (bits_qpsk[0::2]*2 bits_qpsk[1::2]) # 0,1,2,3 symbols_i np.cos(np.pi/2 * symbols_qpsk) # I 分量: 1,0,-1,0 symbols_q np.sin(np.pi/2 * symbols_qpsk) # Q 分量: 0,1,0,-1 # 上采样 脉冲成形 tx_i upfirdn(rrc_taps, symbols_i, sps) tx_q upfirdn(rrc_taps, symbols_q, sps) tx_signal tx_i 1j*tx_q # 复包络这段代码里最需要注意的是 rrc_filter 的系数归一化直接除以能量开根号保证滤波器不改变信号功率。QPSK 映射有多个等价写法上面的角度映射方式最贴近课本。upfirdn 一次完成上采样和滤波注意它默认是“先插值后滤波”输出长度会比你预想的多出滤波器时延后续对齐时要做群时延补偿。3.2 马赫-曾德尔调制器:从理想乘性模型到偏置点设置真实光发射机里电信号不会直接变成光场而是通过马赫-曾德尔调制器加载到光载波上。课程仿真通常用理想模型光场等于直流分量乘以电信号的某种转移函数。但有一个参数你必须设置那就是偏置点。偏置在正交点用于线性调制偏置在最小点用于产生载波抑制信号。仿真里最常见的错误是忘记考虑调制器插入损耗和非线性 transfer 曲线导致接收端幅度和理论预期对不上。# 光载波参数 fc 193.1e12 # 光频率 193.1 THzC band 中心 P0 1e-3 # 1 mW 入纤光功率 # MZM 传递函数正交偏置线性区 V_pi 5.0 V_bias V_pi/2 # 正交偏置点 def mzm_transfer(v_signal): MZM 强度调制传递函数v_signal 为归一化驱动电压 phi np.pi * (v_signal V_bias) / V_pi tf np.cos(phi) return P0 * tf**2 # 将 QPSK 电信号转换为光强度 v_drive tx_i[:len(tx_signal)//2] # 注意长度对齐 optical_power mzm_transfer(v_drive) optical_field np.sqrt(optical_power) * np.exp(1j*2*np.pi*fc*np.arange(len(optical_power))/fs)光场表达式里的 fs 是采样率如果你直接拿光频率除以采样率会发现相位变化太快无法可视化。实际仿真里更推荐用复包络表示即去掉光载波项只在接收端加入本振和频偏。这段代码展示了 MZM 非线性 transfer 的含义cos 函数把输入电压映射成光场幅度平方后才是光功率。仿真时如果发现波形形状畸变先检查偏置点是不是真的在正交点。另外注意 v_drive 长度对齐问题前一级滤波器输出的长度不等于原始符号数乘 sps必须做截断或延迟补偿。4. 光纤信道仿真衰减、色散与非线性的建模顺序与参数陷阱4.1 衰减模型系数、距离与放大器的位置设在哪里光纤损耗是最容易建模的部分α 0.2 dB/km 是标准单模光纤的值1550 nm 窗口下每公里衰减约 4.6%。课程仿真里如果只有衰减误码率曲线就是一条随距离单调恶化的线这就够验证链路预算计算了。放大器的位置在仿真里往往被简化但作为做课程设计的人你需要知道实际系统每隔 80 km 放一个掺铒光纤放大器补偿衰减的同时引入自发辐射噪声。仿真里放大器的位置直接决定噪声的注入点这会显著影响接收端的信噪比。# 光纤衰减 alpha_db_km 0.2 fiber_length_km 100 alpha_linear 10**(-alpha_db_km * fiber_length_km / 10) # 施加衰减 signal_att optical_field * np.sqrt(alpha_linear) # 噪声注入简化EDFA模型放大同时加ASE噪声 G 1/alpha_linear # 增益恰好补偿衰减 ASE_noise np.random.randn(len(signal_att)) 1j*np.random.randn(len(signal_att)) ASE_noise ASE_noise * np.sqrt((G-1) * 0.1) # 噪声系数 NF 相关 rx_after_amp signal_att * np.sqrt(G) ASE_noise这段代码里最容易翻车的是功率单位。光纤衰减在 dB 域是乘法转线性域要除 10 再取指数。ASE 噪声功率的准确计算需要噪声系数这里 0.1 是占位参数真实仿真时要从系统指标换算。放大器增益恰好补偿衰减是理想化设置实际中增益平坦度会造成信道间功率差这是波分复用系统要额外考虑的问题。本段代码适用于单信道课程仿真意义在于建立“衰减-补偿-噪声”的因果关系。4.2 色散仿真群速度色散如何用频域滤波器实现色度色散是单模光纤最主要的线性损伤脉冲不同频率分量传播速度不同导致时域展宽。仿真里最稳妥的实现是频域相位滤波器在频域乘一个相位因子 exp(j * D * λ² * (2πf)² * L / (4πc))。这里的关键是频率轴要对应基带即 f 是相对于光载波的频率偏移而不是绝对频率。很多课程项目的翻车点就在这里。def chromatic_dispersion(signal, fs, D, lam, L): 频域色散滤波器 D: 色散系数 ps/(nm·km)典型值 17 lam: 波长 nm典型值 1550 L: 光纤长度 km N len(signal) freq np.fft.fftfreq(N, 1/fs) # 基带频率 # 离散频点上的色散相位 beta2 -D * lam**2 / (2*np.pi*3e8) # s^2/m phase np.exp(1j * beta2 * (2*np.pi*freq)**2 * L * 1e3 / 2) return np.fft.ifft(np.fft.fft(signal) * phase) # 应用色散 D 17 # ps/(nm·km) lam 1550 # nm L 80 # km signal_disp chromatic_dispersion(rx_after_amp, fs, D, lam, L)这段代码里有三个参数需要特别留意。fs 必须和前面发射机保持一致不能 phonon 频率变了采样率也跟着变。L 的单位是 km但 beta2 算出来是 s²/m所以乘 L 时要乘 1e3 转成米。色散相位里的 2πff 单位是 Hz整个指数因子必须无量纲。如果你发现色散仿真后信号频谱出现折叠或时域波形出现周期性的边带通常原因只有一个频率轴没对齐频谱分辨率不够导致相位差点不连续。4.3 自相位调制什么时候必须仿真什么时候可以忽略非线性效应里课程最常提的是自相位调制光功率通过克尔效应改变折射率进而改变相位。仿真它的标准做法是分步傅里叶法把一段光纤分成很多小步每步先走线性的色散再走非线性的相位旋转。这个算法理解不难但实现有讲究。步长选太大非线性相位会积累误差步长太小计算量暴增。经验值是每步非线性相移控制在 0.01 弧度以内。课程仿真里如果入纤功率小于 0 dBm自相位调制的影响几乎可以忽略仿真里不做非线性部分结果不会有大偏差。但如果功率超过 10 dBm必须加非线性否则误码率曲线会比真实情况乐观很多。有经验的工程师在项目里会先跑一版线性仿真确认链路预算再叠加非线性看恶化程度这个习惯在课程设计里也很适用。5. 相干接收机与误码率统计本振频偏、时钟恢复和完整链路5.1 接收机模型相干检测、90度混频器与平衡探测器相干接收的核心是本振激光器与信号光混频90度混频器分出 I/Q 四路平衡探测器把光信号转为电流。课程仿真里这个模块通常简化为已知发射信号的复包络加上本振频偏和相位噪声再叠加接收机热噪声。你在课本上见过的“接收机灵敏度”和“本振功率”概念在仿真中的体现就是本振相位噪声方差和探测器响应度这两个参数。# 本振与混频 f_lo_offset 100e6 # 本振频偏 100 MHz phase_noise_std 0.01 # 相位噪声标准差 # 本振光场 t np.arange(len(signal_disp)) / fs LO np.exp(1j*(2*np.pi*f_lo_offset*t phase_noise_std*np.random.randn(len(signal_disp)))) # 相干混频简化模型 signal_coh signal_disp * np.conj(LO) # 匹配滤波接收端RRC rx_i upfirdn(rrc_taps, signal_coh.real, 1)[:len(signal_coh)] rx_q upfirdn(rrc_taps, signal_coh.imag, 1)[:len(signal_coh)] # 采样每符号取一个点 rx_sym rx_i[::sps][:len(symbols_qpsk)] 1j*rx_q[::sps][:len(symbols_qpsk)]本振频偏 100 MHz 是故意设置的非理想参数目的是让你看到星座图旋转。如果没有频偏模型星座图就是一个标准正方形看起来完美但和实际脱节。相位噪声标准差 0.01 弧度是低噪声激光器的水平如果调到 0.1你就会看到星座点拖尾连成环。这段代码里的匹配滤波用 upfirdn 时采样率参数设为 1含义是不做上采样只做滤波然后手动抽取每符号一个点。5.2 误码率统计理论曲线对比与蒙特卡洛次数选择误码率统计看起来是数错误比特数除以总比特数但这里有个隐藏的大坑序列对齐。发射端滤波器和接收端滤波器都有时延信道色散也会引入相位旋转直接逐位比较必然得到接近 0.5 的错误率。正确的做法是做互相关对齐或使用已知的训练序列。课程仿真里最省事的方案是计算接收符号和发送符号之间的相位旋转角把星座图转正再做判决。# 星座图相位校正 # 用已知符号估计相位误差 phase_est np.angle(np.sum(rx_sym * np.conj(symbols_qpsk))) print(f估计相位误差: {phase_est:.4f} rad) rx_corrected rx_sym * np.exp(-1j*phase_est) # QPSK 判决 rx_bits_hard np.zeros(len(rx_corrected)*2, dtypeint) for k, sym in enumerate(rx_corrected): # 象限判决 if sym.real 0 and sym.imag 0: rx_bits_hard[2*k:2*k2] [0, 0] elif sym.real 0 and sym.imag 0: rx_bits_hard[2*k:2*k2] [1, 0] elif sym.real 0 and sym.imag 0: rx_bits_hard[2*k:2*k2] [1, 1] else: rx_bits_hard[2*k:2*k2] [0, 1] # BER 计算 ber_measured np.mean(rx_bits_hard ! bits_qpsk) print(f测量 BER {ber_measured:.2e}) # 理论 QPSK BERAWGN 信道 def ber_qpsk_theory(ebn0_db): eta 10**(ebn0_db/10) from scipy.special import erfc return 0.5 * erfc(np.sqrt(eta)) print(f理论 BER (5dB) {ber_qpsk_theory(5):.2e})相位校正用到了所有符号做估计这在课程仿真里是可接受的简化。实际系统里会用前导序列做数据辅助估计但原理一致。蒙特卡洛仿真次数同样影响统计稳定性实测 1e5 个比特时 BER 能在 1e-3 到 1e-4 量级给出三位有效数字1e6 个比特才能把置信区间压到可发表水平。如果你发现 BER 曲线抖动厉害先加比特数不要先怀疑算法。5.3 完整链路脚本参数统一管理与循环扫描把前面所有模块串起来后下一步是扫参。课程作业通常要求画 BER 随信噪比或光纤距离的曲线。这里最忌讳的就是复制粘贴大段代码后改一个数字再跑一次参数统一管理的做法是写配置字典或参数类循环里只改目标参数。def run_link(ebn0_db, L_km): 完整链路返回 BER # 内部使用全局参数返回测量 BER # ... 依次调用发射机、信道、色散、接收机 # 这里用占位返回值实际实现按上文模块拼装 return ber_qpsk_theory(ebn0_db) # 示例返回理论值 # 扫描 EBN0 ebn0_list np.arange(0, 11, 1) ber_list [run_link(e, 80) for e in ebn0_list] import matplotlib.pyplot as plt plt.semilogy(ebn0_list, ber_list, o-) plt.xlabel(Eb/N0 (dB)) plt.ylabel(BER) plt.grid(True) plt.show()注意 run_link 函数体内必须把上文所有模块按顺序调用这里为了行文简洁省略了实际操作时长链路脚本会有 200 到 400 行。保存两套脚本一套是单次运行的链路带星座图和眼图输出另一套是扫描脚本批量改参数并收集结果。这两个脚本互相独立调试时只动第一个出结果时只跑第二个。这是做课程项目时效率最高的组织方式。6. 光纤通信仿真避坑指南五条血泪经验与参数排查清单6.1 采样率与符号率的关系错了全链路静默失真现象星座图看起来就是一团随机点怎么调滤波器都没用。原因采样率不是符号率的整数倍上采样和抽取环节丢错了位置序列错位。解决在全局参数区断言 sps 为整数且 fs/R 大于等于 2最好不小于 4。加一行 assert 比事后调试节省两个小时。assert fs % R 0, 采样率必须是符号率的整数倍 sps fs // R6.2 色散相位方向反了脉冲被压缩而非展宽现象色散后信号波形比发射端更“尖”或者频谱出现不对称畸变。原因beta2 的符号没搞清楚正常单模光纤在 1550 nm 处 beta2 为负即反常色散。计算公式里丢掉负号会直接导致相位方向反转。解决记住一个检查方法看脉冲展宽方向仿真里故意把色散系数设成 -17 对比一次如果波形变化完全对称说明公式写错了。6.3 滤波器时延未补偿BER 卡在 0.25 不动现象QPSK 误码率无论如何优化信噪比都稳定在 0.25 左右。原因发射端和接收端各用了一个 span11 的根升余弦滤波器总时延是 11 个符号周期采样时候没把这段延迟去掉。符号间错位一个符号误码率接近 0.25。解决计算滤波器群时延 span * sps // 2在采样前把信号向后平移这么多点。这是课程仿真里最常见也最隐蔽的问题强烈建议打印发射符号和接收符号的前 20 个值做肉眼对齐。6.4 蒙特卡洛次数不足误码率曲线像锯齿现象BER 曲线理论上该平滑下降实际每次运行都不同加信噪比后反而上升。原因仿真比特数太少统计波动大于信噪比改善带来的差异。解决扫描最低信噪比时用 1e6 比特确保最低误码率点至少有 100 个错误比特。如果跑 1e6 比特耗时太长把符号率降到 1 GBaud 以下计算量主要和采样点数相关和符号率数值没有关系。6.5 相位噪声功率与频偏耦合星座图拖尾成环现象QPSK 星座图上四个点变成四个弧形严重时环形化。原因本振相位噪声方差设得太大或者频偏没有和相位噪声分开建模。激光器线宽与相位噪声方差的关系公式在课本上有仿真时按激光器线宽参数换算不要直接随手填 0.1。解决频偏只做星座图定标时加做 BER 曲线时先补偿频偏再统计否则误码率和理论值的偏差没有物理意义。7. 进阶验证技巧用 EVM 和眼图定位链路瓶颈而不是只盯 BERBER 是最终指标但定位瓶颈时它太迟钝。一个更有效的工作习惯是把误差矢量幅度作为中间监测量在每个模块输出端算一次 EVM这样能瞬间看出哪一级引入的损伤最大。EVM 的定义是误差向量功率与理想符号功率之比的平方根换成 dB 后直观反映“这个模块让信号质量恶化了多少”。这个习惯会让你的课程报告比别人的有深度得多。def evm_db(rx_sym, tx_sym): 计算误差矢量幅度dB err rx_sym - tx_sym evm np.sqrt(np.mean(np.abs(err)**2) / np.mean(np.abs(tx_sym)**2)) return 20 * np.log10(evm) # 链路各节点 EVM 监测 print(f发射端 EVM: {evm_db(rx_sym, symbols_qpsk):.2f} dB) print(f加噪声后 EVM: {evm_db(rx_corrected, symbols_qpsk):.2f} dB)另一个验证技巧是眼图。把接收端信号按符号周期的整数倍切片叠加眼图睁开程度直接反映码间串扰和噪声的叠加效果。课程报告里放一张清晰的眼图比放十个表格都有说服力。做眼图要注意横轴不要画点要画线用短线段连接相邻采样点否则眼图会变成一坨密集点阵。跑完眼图后把眼图开启度与 EVM 做交叉对比你会发现二者高度相关但 EVM 能定量眼图能定性。我的个人习惯是白天调仿真参数晚上统一跑一轮全参数扫描把 EVM 和 BER 曲线放在同一张图上对比看哪条曲线先饱和哪条先劣化据此判断系统是受噪声限制还是受失真限制。这个习惯帮我多次在课程验收前发现了“看起来能跑但对不上理论”的隐患。希望这篇笔记能帮到你少踩几个我踩过的坑。本文还有配套的精品资源点击获取