## 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 regardsHi,

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 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 matlab, fft, optics MATLAB Answers — New Questions