分步傅里叶算法——计算光脉冲在光纤中的传输

function fiberfourier（AD1550Ds1550DpmdraL）
%A为输入光脉冲序列
%D1550为1550nm处的色散系数，单位是ps/nm/km
%Ds1550为1550nm处色散斜率，但是是ps/nm^2/km
%Dpmd为偏振模色散系数，单位是s/km^0.5
%r为光纤的非线性系数，单位是1/（W*km）
%a为光纤损耗，单位是1/km
%L为光纤长度，单位是km
disp（‘调用函数：fiberfourier（AD1550Ds1550DpmdraL）‘）
disp（‘         A为输入光脉冲序列‘）
disp（‘         D1550为1550nm处的色散系数，单位是ps/nm/km‘）
disp（‘         Ds1550为1550nm处色散斜率，但是是ps/nm^2/km‘）
disp（‘         Dpmd为偏振模色散系数，单位是s/km^0.5‘）
disp（‘         r为光纤的非线性系数，单位是1/（W*km）‘）
disp（‘         a为光纤损耗，单位是1/km‘）
disp（‘         L为光纤长度，单位是km‘）
disp（‘程序运行中……‘）
global dt t N f w f0 freqspace numchannel freq
%////////////////////////////////////////////////////////////////
%-------------------------------------脉冲参数
%-----------------C=0;            %啁啾参数
%-----------------num=2;                           %超高斯脉冲阶数
%-----------------T0=100e-12/2/（（log（2））^（0.5/num））;%;1.655;        %脉冲半宽度
%-----------------peak=0.0449;         %脉冲峰值
c=3e8;                                          %speed of optical unit is m/s
wave=c./freq*1e9;        %信道波长%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ch=numchannel;
%------------------------------------时间设置
%-----------------dt=4e-12;
%-----------------t=-200e-12:dt:200e-12;
%-----------------N=length（t）;
%-----------------f=（（（-N/2+0.5）:1:（N/2-0.5））/（dt*N））; %
%-----------------w=2*pi*f;
w0=2*pi*freq;        %1550nm
%-------------------------------------光纤参数
%%-----------------D1550=3.7;                            %1550nm处色散系数ps/nm.km
%%-----------------Ds1550=0.07;                          %1550nm处色散斜率ps/nm2.km
Dwave=D1550+Ds1550*（wave-1550）;
b2=-wave.^2*（1e-9）/（2*pi*3e8）.*Dwave*（1e-12）;         % 二阶色散因子 （s^2/km）
b31=Ds1550*（1e-12）.*wave.^4*（1e-18）/（4*pi^2*9e16）;
b32=-b2.*wave*（1e-9）/（pi*3e8）;    % 三阶色散因子 （s^3/km）
b3=b31+b32;
%%-----------------Dpmd=0.2e-12;%0;              %偏振模色散系数（s/km^1/2）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%-----------------r=2;%0;                % 非线性参数 （1/（w*km））%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%-----------------a=0.275/4.343;        % 损耗参数 （1/km）
%%-----------------L=80;                % 光纤长度 （km）
h=0.2;              % 分成每小段的长度 （km）
m=L/h;              % 段数 
%t0=Dpmd*（3*pi*h/8）^（1/2）;      %每一段光纤主态的延时
b=2*pi*Dpmd*（（3*pi/8/h）^（1/2））*3e14/1.55;
b1=Dpmd/（h^（1/2））;
delay=（b+b1.*（w0-1.2161e15））*h./w0;
%///////////////////////////////////////////////////////////////////
%--------------------------------%信道存在概率
pos=poissrnd（6）;              %通信时间内信道数
posch=[ones（1（pos-1））zeros（1（ch-pos））];
poscho=randperm（ch-1）;
mis=int16（-（25e-12）/dt+（50e-12）/dt*rand（1ch））;    %相邻信道彼特重叠差%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------------
%%-----------------------------------------------A=peak*exp（-（1+i*C）/2*（t）.^（2*num）/T0^（2*num））;    %脉冲波型
%A0=fft（A）;
%A0=abs（fftshift（A0））;
Ay=A/（2^0.5）;
Ax=A/（2^0.5）;
Ax={Ax};
Ay={Ay};
Ax=repmat（Ax1ch）;
Ay=repmat（Ay1ch）;
Aresult=repmat（{}1ch）;
R0=rand（1m）*2*pi;       %偏振旋转方位角
for n=1:m
    f