fenbufuliye.m

利用分步傅里叶方法模拟光纤激光器的输出脉冲与传输距离的关系

clc;
clear all;                   % close all;

wavelength=1550*(10^-9);     %wavelength 单位m
c=3*(10^8);                 %光速 单位:m/s
w0=2*pi*c/wavelength;        %入射波长中心
alpha=9/4.343*(10^-3);        %attenuation coefficient 单位:/m
D=-0.7*(10^-6);              %色散参量   单位s/m^2%-0.7
beta2=0.72*(10^-27);
beta3=0.02*(10^-39);          %三阶色散参数  单位s^3/m
gama=11*(10^-3);             %单位：/w/m
TR=3*(10^-15);               %拉曼增益，单位：s   3fs 
TFWHM=1.6*(10^-12);         %脉冲半高宽度，单位为s   1.6ps
T0=TFWHM/1.763;
P=26.12*(10^-3);               %输入光平均功率，单位：w
P=(2^0.5)*P/TFWHM/50/(10^6);
N=2^15;                     %取样点数
n=0:1:(N-1);  
n1=0:1:N/2-1;
n2=-N/2:1:-1;
h=1e-8;                      %步长，单位:m
Nh=1e2;                     %模拟步长数          
NN=10;                     %记录波形数
Detat=3*10^-15;              %抽样点时间间隔
T=N*Detat;                  %总时间长度，单位为s
t=(n-N/2)*Detat;              %时间序列
t1=n1*Detat;
t2=n2*Detat;
Detaw=2*pi/(T);              %频谱间隔
w=(n-N/2)*Detaw;            %频谱序列
wave=2*pi*c./(w+w0)*10^6;    %波长序列um
D_w=-alpha/2+j*beta2*(w.^2)/2-j*beta3*(w.^3)/6;
UU1=D_w(N/2+1:N);               %改变频率显示顺序
UU2=D_w(1:N/2);
D_w=[UU1,UU2];
U_in=(P)^0.5*sech(t/T0);        %输入时间脉冲的包络,Hyperbolic-Secant Pulse
U_mh=ones(NN+1,N);
U_mhw=ones(NN+1,N);
U=U_in;
Uw=ifft(U,N)*N;              %傅立叶变换
ABS_Uw=Uw.*conj(Uw);
U_mh(1,:)=abs(U).^2;
U_mhw(1,:)=fftshift(ABS_Uw);
for k=1:Nh 
    ABS_U=abs(U).^2;
    ABS_Uw=Uw.*conj(Uw);
   if(mod(k,Nh/NN)==0)
        U_mh(k*NN/Nh+1,:)=abs(U).^2;
        U_mhw(k*NN/Nh+1,:)=fftshift(ABS_Uw);
   end
      Selfcraggedness=i*[ABS_U(1)*U(1)-ABS_U(N)*U(N),diff(ABS_U.*U)]./(w0*Detat*U); 
        ABS_UU=[ABS_U(1)-ABS_U(N),diff(ABS_U)]./Detat;
        Raman_effect=TR*ABS_UU;
        NL=i*gama*(ABS_U+Selfcraggedness-Raman_effect);
        U=ifft(exp(h*NL).*U)*N;   %傅立叶变换
        U=fft(exp(h*D_w).*U)/N;%反傅立叶变换
        Uw=ifft(U)*N;
end
    Maxw=max(U_mhw(1,:));
    Max=max(U_mh(1,:));
    Max11=max(U_mhw(11,:));
%图形输出         

figure   %输出波长谱
plot(wave*1000,10*log10(U_mhw(11,:)/Max11));
title('出射脉冲波长谱');
xlabel('\fontsize{16} Wavelength \lambda(nm)');
ylabel('\fontsize{16} Nomorlized Intensity(dB)');
figure
waterfall(w*T0,[0:NN]*Nh*h/NN,U_mhw/Maxw)
xlabel('(\omega-{\omega}_{0})T_0')
ylabel('Z(m)')
zlabel('相对强度')
title('频谱波形')