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