fbfly.m

利用分步傅里叶方法对光线传输长度对脉冲波形的影响进行比较

clear all
close all
% [e-mail]
%光纤参数          超连续谱的话 N=2^15  detat=10^-15 步长取小一点

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=6.89*(10^-13);   %脉冲半高宽度，单位为s  
T0=TFWHM/1.763;
%T0=2*(10^-13);              %1/e 点宽度0.20皮秒
% P_ave=26.12*10^-3;             %输入光平均功率，单位：w
% P=2^0.5*P_ave/TFWHM/50/10^6;
P=50;    %此处可以改参数  输入峰值脉冲功率 单位：w
%%%%%%%%%%% 入射脉冲参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
P=P^0.5;
N=2^12;              %取样点数 有可能需要改 但必须是2的偶数次方
n=0:(N-1);  
n1=0:N/2-1;
n2=-N/2:1:-1;
h=1e-3;              %步长，单位:m 这个地方可以改的
Nh=2e2;            %模拟步长数 这个地方可以改的
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*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);
% T_rms=sum(t1.*U_mh(1,(1:N/2)))./sum(U_mh(1,(1:N/2)));%用来计算均方差脉宽的 
% T2_rms=sum(t1.^2.*U_mh(1,(1:N/2)))./sum(U_mh(1,(1:N/2)));
% sigma_T(1)=2*sqrt(T2_rms-T_rms.^2);
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);
%         T_rms=sum(t1.*U_mh(k*NN/Nh+1,(1:N/2)))./sum(Umh(11,(1:N/2)));%用来计算均方差脉宽的
%         T2_rms=sum(t1.^2.*U_mh(k*NN/Nh+1,(1:N/2)))./sum(Umh(11,(1:N/2)));
%         sigma_T(k*NN/Nh+1)=2*sqrt(T2_rms-T_rms.^2);
    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,10*log10(U_mhw(11,:)/Max11));
 title('出射脉冲波长谱');
 xlabel('\lambda({\mu}m)');
 ylabel('相对强度');
  figure
 
  waterfall(w*T0,[0:NN]*Nh*h/NN,U_mhw/Maxw)
  xlabel('(\omega-{\omega}_{0})T_0')
  ylabel('Z(m)')
  zlabel('相对强度')
  title('频谱波形')

%   figure
% plot(t/T0, U_mh(11,:));
% xlabel('T/T0');
% ylabel('|U(z,t)|^2');
% title('出射脉冲');
% figure
%plot(sigma_T/sigma_T(1));
%ylabel('sigma_T');