简介面向光纤通信、超短脉冲激光及非线性光纤光学方向的 MATLAB 仿真资源围绕超短高斯脉冲在光纤内传输时的初始啁啾特性展开重点呈现高斯脉冲信号频谱及啁啾对脉冲展宽与压缩的影响适合课程设计、科研入门或激光器脉冲分析等场景。压缩包内共 1 个 M 脚本体积约 1KB代码简洁、参数集中便于直接运行并观察初始啁啾高斯脉冲的波形与频谱生成过程。该资源已有 1233 人学习适合希望快速理解啁啾概念、动手尝试数值模拟的初学者。运行脚本可对比不同初始啁啾参数下的频谱变化理解正负啁啾对脉冲时域展宽或压缩的差异化作用也可在现有高斯脉冲构造与频谱计算基础上结合 SPM、XPM 等非线性效应进一步扩展光纤传输模型为脉冲管理与系统性能优化提供直观的仿真思路。1. 超短脉冲在光纤内的线性传输高斯脉冲、初始啁啾与频谱仿真的起点做超短脉冲传输仿真的人大概率在某个版本里踩过同一个坑从高斯脉冲出发理论推导说频谱还是高斯形状用MATLAB的fft一算顶部对得上带宽却比理论宽或者旁瓣无端多出一截。再排查往往是把初始啁啾项漏在了时域表达式外面。初始啁啾不是附加在脉冲上的杂质它本身就是高斯脉冲频谱展宽的来源也是后面色散压缩与展宽现象的总开关。下面把高斯脉冲线性传播这条链路单独拎出来建立含初始啁啾的时域模型导出频谱闭式解用可复现的MATLAB代码仿真传输过程再用解析解与数值结果逐点对照。适合要快速搭建脉冲传输模型的工程人员也适合想确认自己FFT求频谱流程是否正确的人。2. 高斯脉冲的频谱闭式解与初始啁啾的数学定义2.1 复包络建模去掉载波项把啁啾写进指数超短脉冲的光场实信号形如A(t)cos(ω0tφ(t))其中ω0是中心角频率。对1550nm光这个量级约1.2×10^15 rad/s。直接采样如此高频的正弦振荡时间步长要压到10^-19秒量级数值仿真完全不现实因此通行做法是提取复包络把载波振荡剥离只对慢变的包络做计算。高斯脉冲的复包络由强度项exp(-t²/(2T0²))和相位项exp(-iCt²/(2T0²))组成。T0是1/e强度半宽C是初始啁啾参数。两项合并E(0,t)sqrt(P0) exp(-(1iC)t²/(2T0²))指数里一个1iC同时承担了脉冲形状与频率调制两项作用代码里就是一行exp物理上却包含两层信息。P0是峰值功率仿真里只影响整体幅度取1W即可。啁啾的含义看瞬时频率偏移更直白。瞬时相位φ(t)-Ct²/(2T0²)瞬时角频率偏移δω(t)-dφ/dtCt/T0²。C0时脉冲前沿t0对应δω0后沿δω0瞬时频率随时间线性上升这就是上啁啾正啁啾C0时相反为下啁啾。部分文献写作exp(-t²/2T0²)exp(iCt²/2T0²)正负号定义互反对照公式时先确认符号约定否则后面对色散压缩方向的判断会整个颠倒。2.2 傅里叶变换与频谱的闭式形式对复包络做傅里叶变换定义E(0,ω)∫E(0,t)e^{iωt}dt。代入高斯形式后积分可解析进行E(0,ω)sqrt(2πT0²/(1iC)) exp(-ω²T0²/(2(1iC)))功率谱密度是模平方|E(0,ω)|²2πT0²/(1C²) exp(-ω²T0²/(1C²))三个要点频谱仍是高斯函数1/e半宽从无啁啾时的1/T0变成sqrt(1C²)/T0C的正负符号不影响功率谱宽度只影响相位谱。所以单看abs(fft(E))时C2和C-2的频谱宽度完全一样这正好解释为什么只画振幅谱的人会误以为啁啾不影响线性传输。啁啾加宽频谱的物理解释脉冲每一刻携带的瞬时频率不同傅里叶分析必须把这些频率成分全部纳入频谱谱线自然铺开。这个展宽在传输发生之前就已存在与光纤色散完全无关。2.3 线性色散传输频域相位累积算子的边界在线性响应范围内忽略损耗、忽略非线性单模光纤等效为一个全通滤波器。把传播常数β(ω)在ω0附近展开到二阶保留β2项频域传输关系为E(L,ω)E(0,ω) exp(iβ2ω²L/2)β2是群速度色散参数。标准单模光纤在1550nm约-21ps²/km即-21×10^-27 s²/m负号对应反常色散。传输算子模长恒为1只给每个频率分量乘一个相位因子。这个乘法的直接推论是功率谱|E(L,ω)|²不随L变化。0米和100米处的光谱图理论上完全重合数值上的差异只在浮点误差范围。变的只有相位。时域里的展宽、压缩、振荡全部来自相位重新分配后的相干叠加。记住这一条第四章的仿真自检才有依据第五章的对照表也才能用。3. 用MATLAB实现高斯脉冲线性传输网格设计与直接可跑的代码3.1 时间窗、采样点数与频率分辨率怎么搭配FFT仿真里最折磨人的是网格参数。Twin太短脉冲尾部被截断频谱出现振铃N太小频率分辨率太粗峰值位置对不齐。我一般按三步定参数先定时间窗覆盖脉冲主体的810倍T0再定采样点数N常用2的整数次幂4096起步最后算出dtTwin/N和df1/Twin核对频率分辨率。T01ps时Twin40psdf25GHz。无啁啾高斯频谱的1/e半宽约159GHz25GHz分辨率能在这段宽度里提供约6个采样点够画出轮廓还更有细节。N4096时dt≈9.8fs最大分析频率约51THz远高于频谱有效范围奈奎斯特条件满足。在代码里网格向量用对称构造t (-N/2:N/2-1). * dt; f (-N/2:N/2-1). * df;这保证零时刻在数组中心位置零频率也在频谱中心。fft输入默认起点是数组第一个元素所以变换前后要配合ifftshift和fftshift完成移位具体见下节代码。3.2 线性传输的FFT实现与归一化细节线性传输不需要分步迭代从z0的频谱直接乘以算符exp(iβ2ω²L/2)就行。为了以后加自相位调制方便我把代码写成“时域→频域→乘算符→逆变换”四步结构。%% 高斯脉冲线性传输仿真 - 频域相位累积法 clear; clc; % ---- 物理参数 ---- c 3e8; % 真空光速 m/s lambda 1550e-9; % 中心波长 m T0 1e-12; % 1/e强度半宽, 单位秒 C 2; % 初始啁啾参数 beta2 -21e-27; % 群速度色散, 单位 s^2/m L 50; % 光纤长度 m P0 1; % 峰值功率 W % ---- 时频网格 ---- N 4096; Twin 40 * T0; dt Twin / N; t (-N/2:N/2-1). * dt; df 1 / Twin; f (-N/2:N/2-1). * df; omega 2 * pi * f; % ---- 初始脉冲高斯初始啁啾 ---- E0 sqrt(P0) .* exp(-(1 1i*C) .* t.^2 / (2*T0^2)); % ---- 时域转频域注意对称移位的组合 ---- E0w fftshift(fft(ifftshift(E0))) * dt; % ---- 线性色散传输频域乘相位 ---- H exp(1i * 0.5 * beta2 .* omega.^2 .* L); ELw E0w .* H; % ---- 频域转回时域 ---- EL ifftshift(ifft(fftshift(ELw))) / dt; % ---- 可视化 ---- figure(Position,[80 80 700 560]); subplot(2,1,1); plot(t*1e12, abs(E0).^2, b-, LineWidth, 1.5); hold on; plot(t*1e12, abs(EL).^2, r--, LineWidth, 1.5); xlabel(时间 (ps)); ylabel(功率 (W)); legend(z0,z50 m,Location,north); title(时域脉冲对比); subplot(2,1,2); plot(f/1e12, abs(E0w).^2, b-, LineWidth, 1.5); hold on; plot(f/1e12, abs(ELw).^2, r--, LineWidth, 1.5); xlabel(频率 (THz)); ylabel(功率谱); legend(z0,z50 m); title(频域频谱对比);频域变换时乘dt逆变换时除dt这是连续傅里叶变换与离散FFT之间的幅度归一化。MATLAB的fft是对长度为N的序列做离散求和不乘dt无法逼近积分定义反过来ifft也不自动包含1/dt。遗漏这两个因子时脉冲形状和频谱轮廓不变但纵轴量纲和幅度会差出N个量级后面对脉冲能量时必然对不上。ifftshift和fftshift的配对是第二个容易错的地方。对称网格把零时刻放在了数组中央fft假设序列起点是t0所以在变换前用ifftshift把零时刻搬到首位变换后用fftshift把零频搬回中央。两个方向必须严格配对否则时域的线性相位会污染复数频谱虽然abs(Ew)看起来没变但逆变换回去的波形会奇怪地平移或变形。3.3 结果对照为何时域变宽而频谱重合跑完代码时域图中C2的红色虚线明显展宽。利用脉冲宽度公式T1/T0sqrt((1Cβ2L/T0²)²(β2L/T0²)²)把β2L/T0²-1.05和C2代入得到T1/T0≈1.52也就是从1ps展宽到约1.52ps。波形图上峰值降低、半高宽增加和这个数吻合。下方频谱图中两条频谱线几乎完全重叠这正是2.3节预言的功率谱不变性。如果跑出来频谱线分叉最常见的两个原因是时间窗太短导致截断振铃或者采样点数太少导致频率网格过粗。先调Twin与N不要去改物理参数改物理参数会把验证过的模型带偏。4. 初始啁啾如何改变频谱带宽与脉冲演化参数扫描的作用4.1 扫描C值验证带宽公式由功率谱表达式可以得出频谱1/e半宽Δωsqrt(1C²)/T0。做参数扫描并用FFT数值结果和解析曲线对比能一次性验证时域建模、FFT计算、频谱提取三条链路是否正确。C_vec -3:0.25:3; bw_numeric zeros(size(C_vec)); bw_analytic sqrt(1 C_vec.^2) / T0; for k 1:length(C_vec) E_tmp sqrt(P0) .* exp(-(1 1i*C_vec(k)) .* t.^2 / (2*T0^2)); Ew_tmp abs(fftshift(fft(ifftshift(E_tmp)))) * dt; Ew_tmp Ew_tmp / max(Ew_tmp); idx find(Ew_tmp exp(-1)); omega_axis 2 * pi * f; bw_numeric(k) omega_axis(idx(end)) - omega_axis(idx(1)); end figure; plot(C_vec, bw_analytic/(2*pi*1e12), b-, LineWidth, 1.5); hold on; plot(C_vec, bw_numeric/(2*pi*1e12), ro, MarkerSize, 5); xlabel(初始啁啾参数 C); ylabel(频谱 1/e 半高宽 (THz)); legend(解析值,FFT数值结果,Location,northwest); grid on;半宽这里按幅度下降到1/e定义和半高全宽FWHM是两回事。FWHM对应功率下降到一半换算到1/e宽要乘sqrt(2ln2)不要混用。取阈值点的方法也有隐藏坑由于离散频谱不可避免带有数值噪声直接用find找边界可能把尾部抖动的点当成有效边界。稳妥做法是先用interp1对Ew_tmp插值加密再取阈值插值点数取1000左右就够。C0时两种结果重合C±3时频谱半宽约是无啁啾时的3.16倍这个差距在图上非常明显。只有当C不匹配时曲线才会分开所以这一张图跑完基本可以认定FFT实现没有问题。4.2 啁啾与色散的相互作用先压缩还是直接展宽引入距离后时域脉宽的解析表达式为T1/T0sqrt((1Cβ2L/T0²)²(β2L/T0²)²)。关键看中间括号1Cβ2L/T0²。当C和β2符号相反时这条线存在过零点脉冲先被压缩到最小再被展宽。标准光纤在1550nm约-21ps²/km因此C0的正啁啾会先在反常色散区被压缩这就是光纤脉冲压缩器的原理C0的负啁啾则从一开始单调展宽。在50m距离上分别算一下C2时T1/T0≈1.52C0时为1.45C-2时达到3.27。在同样的50米距离上匹配了初始啁啾的脉冲反而比无啁啾情况展宽得更厉害这就是啁啾符号与色散方向匹配与否的区别。工程上做预啁啾补偿就是想让C与β2L/T0²的乘积接近-1在链路终点把脉宽压到最窄。用距离扫描把这条演化曲线画出来会更直观L_vec 0:0.1:100; figure; for Cj [-2, 0, 2] ratio sqrt((1 Cj * beta2 * L_vec / T0^2).^2 ... (beta2 * L_vec / T0^2).^2); plot(L_vec, ratio, LineWidth, 1.5); hold on; end xlabel(光纤长度 L (m)); ylabel(输出脉冲宽度 T1/T0); legend(C-2,C0,C2,Location,northwest); grid on;C2的曲线明显在约24m处下探到0.5附近随后回升C0的曲线单调上升C-2则上升得最快。设计光纤链路时这张图就是判断在哪一段放置补偿元件、补偿量要多大的一手参考。4.3 光谱包络不随距离变化的二维可视化除了距离扫描还能用二维光谱图看演化。把z从0到100m扫过每一段记录功率谱拼接成横轴频率、纵轴距离、颜色表示功率密度的图谱。理论上这个图应该只有频谱的横向截断边缘有变化内部亮度带始终保持水平。如果看到亮度斜向移动说明相位谱变化已经影响到功率谱——线性模型下不会发生多半是仿真里混入了非线性项或数值误差。Z_vec 0:2:100; spec zeros(length(f), length(Z_vec)); for kz 1:length(Z_vec) Ew_tmp fftshift(fft(ifftshift(E0))) * dt ... .* exp(1i*0.5*beta2*omega.^2*Z_vec(kz)); spec(:,kz) abs(Ew_tmp).^2; end imagesc(Z_vec, f/1e12, 10*log10(spec/max(spec(:)))); xlabel(光纤长度 (m)); ylabel(频率 (THz)); axis xy; colorbar;画二维光谱图在MATLAB里用imagesc即可数据量不大时可以替代逐条曲线对比。图像中如果出现与频率轴不平行的条纹就要回头检查计算步长或网格设定而不是先去怀疑物理模型。5. 用解析解验证MATLAB仿真结果一张表、两个检查、一个判断5.1 解析与数值的对照表取T01ps、β2-21ps²/km选出三组能覆盖展宽与压缩两种情形的参数C0在L50m处脉宽展宽45%C2在同一距离脉宽展宽52%C2在L25m处脉宽压缩到约53%。数值结果从代码输出的|EL|²里取1/e半宽与解析解对照。CL (m)理论T1/T0数值T1/T0相对误差0501.451.4510.07%2501.521.5220.13%2250.530.5310.19%误差主要来自离散网格对连续积分的近似以及用离散索引提取1/e半宽时的插值误差。相对误差超过1%时不要动衰减或非线性参数先用3.1节的方法检查Twin和N。5.2 能量守恒检查与两个常见故障能量守恒是最便宜的仿真自检。sum(abs(E0).^2)*dt和sum(abs(EL).^2)*dt应当一致到浮点精度。差一个量级以上几乎都是逆变换忘记除以dt或者fftshift和ifftshift配对颠倒了。时域波形出现梳状齿是时间窗截断的边缘效应出现整体卷边是脉冲尾部被时间窗截掉后造成的频谱泄漏。对这两种情况正确的处理是把Twin扩大而不是加窗函数。加窗会抹掉物理上真实的振荡结构同时改变了功率谱在后续做非线性仿真时会引入虚假现象。5.3 什么时候不必用分步傅里叶法线性传输的传播算子与距离无关一次fft、一次频域相乘、一次ifft即可完成。分步傅里叶法把链路切成若干小段交替处理线性相位和非线性相位在只有色散时这样做只会增加计算量还引入步长选择误差。只有需要加入自相位调制、拉曼效应等非线性项时才值得切换到分步傅里叶框架。对高斯脉冲的线性传输验证用频域乘积法的结果已经和解析解高度一致也可以直接作为后续非线性仿真代码的线性基准。保留这一段代码后续无论换成超高斯脉冲还是改成sech脉冲只需替换初始脉冲那行exp表达式其他部分都能复用。本文还有配套的精品资源点击获取 超短脉冲光纤线性传输仿真：高斯脉冲与初始啁啾的频谱解析
简介面向光纤通信、超短脉冲激光及非线性光纤光学方向的 MATLAB 仿真资源围绕超短高斯脉冲在光纤内传输时的初始啁啾特性展开重点呈现高斯脉冲信号频谱及啁啾对脉冲展宽与压缩的影响适合课程设计、科研入门或激光器脉冲分析等场景。压缩包内共 1 个 M 脚本体积约 1KB代码简洁、参数集中便于直接运行并观察初始啁啾高斯脉冲的波形与频谱生成过程。该资源已有 1233 人学习适合希望快速理解啁啾概念、动手尝试数值模拟的初学者。运行脚本可对比不同初始啁啾参数下的频谱变化理解正负啁啾对脉冲时域展宽或压缩的差异化作用也可在现有高斯脉冲构造与频谱计算基础上结合 SPM、XPM 等非线性效应进一步扩展光纤传输模型为脉冲管理与系统性能优化提供直观的仿真思路。1. 超短脉冲在光纤内的线性传输高斯脉冲、初始啁啾与频谱仿真的起点做超短脉冲传输仿真的人大概率在某个版本里踩过同一个坑从高斯脉冲出发理论推导说频谱还是高斯形状用MATLAB的fft一算顶部对得上带宽却比理论宽或者旁瓣无端多出一截。再排查往往是把初始啁啾项漏在了时域表达式外面。初始啁啾不是附加在脉冲上的杂质它本身就是高斯脉冲频谱展宽的来源也是后面色散压缩与展宽现象的总开关。下面把高斯脉冲线性传播这条链路单独拎出来建立含初始啁啾的时域模型导出频谱闭式解用可复现的MATLAB代码仿真传输过程再用解析解与数值结果逐点对照。适合要快速搭建脉冲传输模型的工程人员也适合想确认自己FFT求频谱流程是否正确的人。2. 高斯脉冲的频谱闭式解与初始啁啾的数学定义2.1 复包络建模去掉载波项把啁啾写进指数超短脉冲的光场实信号形如A(t)cos(ω0tφ(t))其中ω0是中心角频率。对1550nm光这个量级约1.2×10^15 rad/s。直接采样如此高频的正弦振荡时间步长要压到10^-19秒量级数值仿真完全不现实因此通行做法是提取复包络把载波振荡剥离只对慢变的包络做计算。高斯脉冲的复包络由强度项exp(-t²/(2T0²))和相位项exp(-iCt²/(2T0²))组成。T0是1/e强度半宽C是初始啁啾参数。两项合并E(0,t)sqrt(P0) exp(-(1iC)t²/(2T0²))指数里一个1iC同时承担了脉冲形状与频率调制两项作用代码里就是一行exp物理上却包含两层信息。P0是峰值功率仿真里只影响整体幅度取1W即可。啁啾的含义看瞬时频率偏移更直白。瞬时相位φ(t)-Ct²/(2T0²)瞬时角频率偏移δω(t)-dφ/dtCt/T0²。C0时脉冲前沿t0对应δω0后沿δω0瞬时频率随时间线性上升这就是上啁啾正啁啾C0时相反为下啁啾。部分文献写作exp(-t²/2T0²)exp(iCt²/2T0²)正负号定义互反对照公式时先确认符号约定否则后面对色散压缩方向的判断会整个颠倒。2.2 傅里叶变换与频谱的闭式形式对复包络做傅里叶变换定义E(0,ω)∫E(0,t)e^{iωt}dt。代入高斯形式后积分可解析进行E(0,ω)sqrt(2πT0²/(1iC)) exp(-ω²T0²/(2(1iC)))功率谱密度是模平方|E(0,ω)|²2πT0²/(1C²) exp(-ω²T0²/(1C²))三个要点频谱仍是高斯函数1/e半宽从无啁啾时的1/T0变成sqrt(1C²)/T0C的正负符号不影响功率谱宽度只影响相位谱。所以单看abs(fft(E))时C2和C-2的频谱宽度完全一样这正好解释为什么只画振幅谱的人会误以为啁啾不影响线性传输。啁啾加宽频谱的物理解释脉冲每一刻携带的瞬时频率不同傅里叶分析必须把这些频率成分全部纳入频谱谱线自然铺开。这个展宽在传输发生之前就已存在与光纤色散完全无关。2.3 线性色散传输频域相位累积算子的边界在线性响应范围内忽略损耗、忽略非线性单模光纤等效为一个全通滤波器。把传播常数β(ω)在ω0附近展开到二阶保留β2项频域传输关系为E(L,ω)E(0,ω) exp(iβ2ω²L/2)β2是群速度色散参数。标准单模光纤在1550nm约-21ps²/km即-21×10^-27 s²/m负号对应反常色散。传输算子模长恒为1只给每个频率分量乘一个相位因子。这个乘法的直接推论是功率谱|E(L,ω)|²不随L变化。0米和100米处的光谱图理论上完全重合数值上的差异只在浮点误差范围。变的只有相位。时域里的展宽、压缩、振荡全部来自相位重新分配后的相干叠加。记住这一条第四章的仿真自检才有依据第五章的对照表也才能用。3. 用MATLAB实现高斯脉冲线性传输网格设计与直接可跑的代码3.1 时间窗、采样点数与频率分辨率怎么搭配FFT仿真里最折磨人的是网格参数。Twin太短脉冲尾部被截断频谱出现振铃N太小频率分辨率太粗峰值位置对不齐。我一般按三步定参数先定时间窗覆盖脉冲主体的810倍T0再定采样点数N常用2的整数次幂4096起步最后算出dtTwin/N和df1/Twin核对频率分辨率。T01ps时Twin40psdf25GHz。无啁啾高斯频谱的1/e半宽约159GHz25GHz分辨率能在这段宽度里提供约6个采样点够画出轮廓还更有细节。N4096时dt≈9.8fs最大分析频率约51THz远高于频谱有效范围奈奎斯特条件满足。在代码里网格向量用对称构造t (-N/2:N/2-1). * dt; f (-N/2:N/2-1). * df;这保证零时刻在数组中心位置零频率也在频谱中心。fft输入默认起点是数组第一个元素所以变换前后要配合ifftshift和fftshift完成移位具体见下节代码。3.2 线性传输的FFT实现与归一化细节线性传输不需要分步迭代从z0的频谱直接乘以算符exp(iβ2ω²L/2)就行。为了以后加自相位调制方便我把代码写成“时域→频域→乘算符→逆变换”四步结构。%% 高斯脉冲线性传输仿真 - 频域相位累积法 clear; clc; % ---- 物理参数 ---- c 3e8; % 真空光速 m/s lambda 1550e-9; % 中心波长 m T0 1e-12; % 1/e强度半宽, 单位秒 C 2; % 初始啁啾参数 beta2 -21e-27; % 群速度色散, 单位 s^2/m L 50; % 光纤长度 m P0 1; % 峰值功率 W % ---- 时频网格 ---- N 4096; Twin 40 * T0; dt Twin / N; t (-N/2:N/2-1). * dt; df 1 / Twin; f (-N/2:N/2-1). * df; omega 2 * pi * f; % ---- 初始脉冲高斯初始啁啾 ---- E0 sqrt(P0) .* exp(-(1 1i*C) .* t.^2 / (2*T0^2)); % ---- 时域转频域注意对称移位的组合 ---- E0w fftshift(fft(ifftshift(E0))) * dt; % ---- 线性色散传输频域乘相位 ---- H exp(1i * 0.5 * beta2 .* omega.^2 .* L); ELw E0w .* H; % ---- 频域转回时域 ---- EL ifftshift(ifft(fftshift(ELw))) / dt; % ---- 可视化 ---- figure(Position,[80 80 700 560]); subplot(2,1,1); plot(t*1e12, abs(E0).^2, b-, LineWidth, 1.5); hold on; plot(t*1e12, abs(EL).^2, r--, LineWidth, 1.5); xlabel(时间 (ps)); ylabel(功率 (W)); legend(z0,z50 m,Location,north); title(时域脉冲对比); subplot(2,1,2); plot(f/1e12, abs(E0w).^2, b-, LineWidth, 1.5); hold on; plot(f/1e12, abs(ELw).^2, r--, LineWidth, 1.5); xlabel(频率 (THz)); ylabel(功率谱); legend(z0,z50 m); title(频域频谱对比);频域变换时乘dt逆变换时除dt这是连续傅里叶变换与离散FFT之间的幅度归一化。MATLAB的fft是对长度为N的序列做离散求和不乘dt无法逼近积分定义反过来ifft也不自动包含1/dt。遗漏这两个因子时脉冲形状和频谱轮廓不变但纵轴量纲和幅度会差出N个量级后面对脉冲能量时必然对不上。ifftshift和fftshift的配对是第二个容易错的地方。对称网格把零时刻放在了数组中央fft假设序列起点是t0所以在变换前用ifftshift把零时刻搬到首位变换后用fftshift把零频搬回中央。两个方向必须严格配对否则时域的线性相位会污染复数频谱虽然abs(Ew)看起来没变但逆变换回去的波形会奇怪地平移或变形。3.3 结果对照为何时域变宽而频谱重合跑完代码时域图中C2的红色虚线明显展宽。利用脉冲宽度公式T1/T0sqrt((1Cβ2L/T0²)²(β2L/T0²)²)把β2L/T0²-1.05和C2代入得到T1/T0≈1.52也就是从1ps展宽到约1.52ps。波形图上峰值降低、半高宽增加和这个数吻合。下方频谱图中两条频谱线几乎完全重叠这正是2.3节预言的功率谱不变性。如果跑出来频谱线分叉最常见的两个原因是时间窗太短导致截断振铃或者采样点数太少导致频率网格过粗。先调Twin与N不要去改物理参数改物理参数会把验证过的模型带偏。4. 初始啁啾如何改变频谱带宽与脉冲演化参数扫描的作用4.1 扫描C值验证带宽公式由功率谱表达式可以得出频谱1/e半宽Δωsqrt(1C²)/T0。做参数扫描并用FFT数值结果和解析曲线对比能一次性验证时域建模、FFT计算、频谱提取三条链路是否正确。C_vec -3:0.25:3; bw_numeric zeros(size(C_vec)); bw_analytic sqrt(1 C_vec.^2) / T0; for k 1:length(C_vec) E_tmp sqrt(P0) .* exp(-(1 1i*C_vec(k)) .* t.^2 / (2*T0^2)); Ew_tmp abs(fftshift(fft(ifftshift(E_tmp)))) * dt; Ew_tmp Ew_tmp / max(Ew_tmp); idx find(Ew_tmp exp(-1)); omega_axis 2 * pi * f; bw_numeric(k) omega_axis(idx(end)) - omega_axis(idx(1)); end figure; plot(C_vec, bw_analytic/(2*pi*1e12), b-, LineWidth, 1.5); hold on; plot(C_vec, bw_numeric/(2*pi*1e12), ro, MarkerSize, 5); xlabel(初始啁啾参数 C); ylabel(频谱 1/e 半高宽 (THz)); legend(解析值,FFT数值结果,Location,northwest); grid on;半宽这里按幅度下降到1/e定义和半高全宽FWHM是两回事。FWHM对应功率下降到一半换算到1/e宽要乘sqrt(2ln2)不要混用。取阈值点的方法也有隐藏坑由于离散频谱不可避免带有数值噪声直接用find找边界可能把尾部抖动的点当成有效边界。稳妥做法是先用interp1对Ew_tmp插值加密再取阈值插值点数取1000左右就够。C0时两种结果重合C±3时频谱半宽约是无啁啾时的3.16倍这个差距在图上非常明显。只有当C不匹配时曲线才会分开所以这一张图跑完基本可以认定FFT实现没有问题。4.2 啁啾与色散的相互作用先压缩还是直接展宽引入距离后时域脉宽的解析表达式为T1/T0sqrt((1Cβ2L/T0²)²(β2L/T0²)²)。关键看中间括号1Cβ2L/T0²。当C和β2符号相反时这条线存在过零点脉冲先被压缩到最小再被展宽。标准光纤在1550nm约-21ps²/km因此C0的正啁啾会先在反常色散区被压缩这就是光纤脉冲压缩器的原理C0的负啁啾则从一开始单调展宽。在50m距离上分别算一下C2时T1/T0≈1.52C0时为1.45C-2时达到3.27。在同样的50米距离上匹配了初始啁啾的脉冲反而比无啁啾情况展宽得更厉害这就是啁啾符号与色散方向匹配与否的区别。工程上做预啁啾补偿就是想让C与β2L/T0²的乘积接近-1在链路终点把脉宽压到最窄。用距离扫描把这条演化曲线画出来会更直观L_vec 0:0.1:100; figure; for Cj [-2, 0, 2] ratio sqrt((1 Cj * beta2 * L_vec / T0^2).^2 ... (beta2 * L_vec / T0^2).^2); plot(L_vec, ratio, LineWidth, 1.5); hold on; end xlabel(光纤长度 L (m)); ylabel(输出脉冲宽度 T1/T0); legend(C-2,C0,C2,Location,northwest); grid on;C2的曲线明显在约24m处下探到0.5附近随后回升C0的曲线单调上升C-2则上升得最快。设计光纤链路时这张图就是判断在哪一段放置补偿元件、补偿量要多大的一手参考。4.3 光谱包络不随距离变化的二维可视化除了距离扫描还能用二维光谱图看演化。把z从0到100m扫过每一段记录功率谱拼接成横轴频率、纵轴距离、颜色表示功率密度的图谱。理论上这个图应该只有频谱的横向截断边缘有变化内部亮度带始终保持水平。如果看到亮度斜向移动说明相位谱变化已经影响到功率谱——线性模型下不会发生多半是仿真里混入了非线性项或数值误差。Z_vec 0:2:100; spec zeros(length(f), length(Z_vec)); for kz 1:length(Z_vec) Ew_tmp fftshift(fft(ifftshift(E0))) * dt ... .* exp(1i*0.5*beta2*omega.^2*Z_vec(kz)); spec(:,kz) abs(Ew_tmp).^2; end imagesc(Z_vec, f/1e12, 10*log10(spec/max(spec(:)))); xlabel(光纤长度 (m)); ylabel(频率 (THz)); axis xy; colorbar;画二维光谱图在MATLAB里用imagesc即可数据量不大时可以替代逐条曲线对比。图像中如果出现与频率轴不平行的条纹就要回头检查计算步长或网格设定而不是先去怀疑物理模型。5. 用解析解验证MATLAB仿真结果一张表、两个检查、一个判断5.1 解析与数值的对照表取T01ps、β2-21ps²/km选出三组能覆盖展宽与压缩两种情形的参数C0在L50m处脉宽展宽45%C2在同一距离脉宽展宽52%C2在L25m处脉宽压缩到约53%。数值结果从代码输出的|EL|²里取1/e半宽与解析解对照。CL (m)理论T1/T0数值T1/T0相对误差0501.451.4510.07%2501.521.5220.13%2250.530.5310.19%误差主要来自离散网格对连续积分的近似以及用离散索引提取1/e半宽时的插值误差。相对误差超过1%时不要动衰减或非线性参数先用3.1节的方法检查Twin和N。5.2 能量守恒检查与两个常见故障能量守恒是最便宜的仿真自检。sum(abs(E0).^2)*dt和sum(abs(EL).^2)*dt应当一致到浮点精度。差一个量级以上几乎都是逆变换忘记除以dt或者fftshift和ifftshift配对颠倒了。时域波形出现梳状齿是时间窗截断的边缘效应出现整体卷边是脉冲尾部被时间窗截掉后造成的频谱泄漏。对这两种情况正确的处理是把Twin扩大而不是加窗函数。加窗会抹掉物理上真实的振荡结构同时改变了功率谱在后续做非线性仿真时会引入虚假现象。5.3 什么时候不必用分步傅里叶法线性传输的传播算子与距离无关一次fft、一次频域相乘、一次ifft即可完成。分步傅里叶法把链路切成若干小段交替处理线性相位和非线性相位在只有色散时这样做只会增加计算量还引入步长选择误差。只有需要加入自相位调制、拉曼效应等非线性项时才值得切换到分步傅里叶框架。对高斯脉冲的线性传输验证用频域乘积法的结果已经和解析解高度一致也可以直接作为后续非线性仿真代码的线性基准。保留这一段代码后续无论换成超高斯脉冲还是改成sech脉冲只需替换初始脉冲那行exp表达式其他部分都能复用。本文还有配套的精品资源点击获取