Split step Fourier and shannon sampling theorem issue

Hi,

I am building a code to model a laser with multiple step optical pulse amplifiers, however when I increase the gain I tend to get spectral aliasing. Little explanation of how the code is supposed to work : I get a pulse of a few ps travelling through optical fibers, due to the non linear contribution of the Kerr effect we observe spectral broadening. And this spectral broadening depends on the energy of the pulse and on the length of fiber travelled. Meaning when I get it through an amplifier stage I get more spectrum.

Here is what I do not get, I defined my sampling parameters this way :

Nsamples=2^16;
TimeWindow=39.5;

Timewin=tau0*TimeWINDOW;    %%time window defined as TimeWINDOW*the width of the pulse
dtau=Timewin/Nsamples;      %%time steps
Fs=1/dtau;                  %%Sampling frequency -> Must be superior to 2*Fmax
dF=1/Timewin;               %%Frequency window
domega=2*pi*dF;             %%angular freq

It works fine until I increase the gain in my amplifier SSF code. However as I understand it my Fmax should be the highest frequency of my window so the lowest wavelength of my spectrum? With I am way under the criteria so I guess a little spectral broadening is enough to get me some spectral aliasing. The issue is I can increase the sampling frequency by reducing the time window but when I reach the criteria my pulse in the time domain can not be described accurately due to its width being larger than the window. Also the closer I get to the frequency the faster my results get damaged.

For TimeWINDOW=10

For TimeWINDOW=5

So I need to increase the amount of samples... I get Nsamples to 2^18 and TimeWINDOW to 15 It should satisfy the criteria yet I still get the following error

And if I increase the gain it will still come faster.

I am kind of lost as to where is my mistake ? I am pretty sure it is an obvious one. Here is how I defined my functions :

%%Pulse building%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[u_in,s_in,Domega,tau,dtau,omega0,Lambda,P00,Dfreq]=Pulse_generation(pulse_form,Pmoy,freq_pulse,Chirp,tau0,Nsamples,TimeWINDOW,lambda0)
Timewin=tau0*TimeWINDOW;    %sqrt(Nsamples)*
dtau=Timewin/Nsamples;              %time sep samples
Fs=1/dtau;                            %frequence d'echantillonage
dF=1/Timewin;                         %sep freq d'echantillonage
domega=2*pi*dF;                     %angular freq

%%time space
tau=(-Nsamples/2:(Nsamples/2)-1)*dtau*1e12;

%%frequency space
Dfreq=(-Nsamples/2:1:(Nsamples/2)-1)*dF;

%%spectral space
omega0=2*pi*3e8/lambda0;
dL=2*pi*3e8/omega0^2*domega*1e9;
Domega=(-Nsamples/2:Nsamples/2-1)*domega;
Lambda=((-Nsamples/2:Nsamples/2-1)*dL+lambda0*1e9);


tauFWHM=tau0;
tau0=tauFWHM/(2*sqrt(log(2)));
%%pulse defined
u_gauss=exp(-(1+1i*Chirp).*(((tau.*1e-12).^2)./(2*tau0^2))); %Impulsion gaussienne 
u_sech=exp(-(1i*Chirp).*(tau.^2)./(2*tau0^2)).*sech(tau/tau0); %%sech hyp
if(pulse_form == 0)
    u_in=u_gauss;
elseif(pulse_form == 1)
    u_in=u_sech;
end

%%Normal pulse
%%P0=40e-6;%%probleme de val de P0
%%P0
I00=u_in.^2;
E01=sum(u_in)*dtau;
E0=Pmoy/freq_pulse;
I01=I00./E01*E0;
P00=E0./I00;
u_in=sqrt(E0/E01)*I01;
P00=(freq_pulse)^-1/(sum(u_in.^2)*dtau)*Pmoy;%%erreur dim??
u_in=sqrt(P00)*u_in;

s_in=fftshift(fft(u_in))/Nsamples;
end
%%%SSF%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u_out,s_out,L_u,L_s]=SSF(u_in,L,beta2,gamma,StepSpatial,Domega,att,g,MapDataSpatial)
    [u_out,s_out]=DISP(u_in,beta2,StepSpatial/2,Domega,g);
    [u_out,s_out]=KERR(u_out,gamma,StepSpatial);
    [u_out,s_out]=DISP(u_out,beta2,StepSpatial/2,Domega,g);
    if(MapDataSpatial==1)
        L_u=u_out;
        L_s=s_out;
    end
    iii=1;
    for ii=2*StepSpatial:StepSpatial:L
        iii=iii+1;

        [u_out,s_out]=DISP(u_out,beta2,StepSpatial/2,Domega,g);
        [u_out,s_out]=KERR(u_out,gamma,StepSpatial);
        [u_out,s_out]=DISP(u_out,beta2,StepSpatial/2,Domega,g);
        if(MapDataSpatial==1)
            L_u=[L_u;u_out];
            L_s=[L_s;s_out];
        end
    end
    if(MapDataSpatial~=1)
            L_u=0;
            L_s=0;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u_out,s_out]=DISP(u_in,beta2,L,Domega,g)
   % mask=zeros(1,length(Domega));
   % mask(1000:length(Domega)-1000)=1;
   Phidisp=exp(1i*((beta2)/2)*Domega.^2*L+g*L);
   s_out=(fftshift(fft(u_in)))/length(u_in).*Phidisp;
   %%s_out=s_out.*mask;
   u_out=(ifft(ifftshift(s_out)))*length(u_in);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Function non linearity only%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [u_out,s_out]=KERR(u_in,gamma,L)
    PhiNL=gamma.*abs(u_in).^2*L;
    u_out=u_in.*exp(1i.*PhiNL);%%(g-att)*L
    s_out=(fftshift(fft(u_out)))/length(u_out);
end
%%%%%%%

g is the gain, gamma the kerr nonlinear coefficient, beta2 the dispersion coefficient.

Also I have to keep a large frequency bandwidth as I have to include Raman effects later on (Raman effect builds up energy at another wavelength around 13.2THz away so around 1081 nm for me I work at 1030nm).

Best regards

[dpb]

Nyquist of Fmax/2 is the theoretical minimum.

From a practical standpoint, one should sample at 3-4X the Fmax. With such high frequency signals, you may well run into memory issues if you need more than very short time pulses; I didn't try to read the whole novel...

[Hugo]

Thanks for the input,

It works now, I lost a lot in computational time but whatever.

Also if anyone ever face the same gain issues, the reason why the gain increase was so impactfull on the overall error is because the non linear contribution to the phase is dependant on the gain leading to my spatial steps no longer being adapted (that was quite stupid of me) so I switched to a variable spatial step.