I have been trying to obtain a spectrum and a spectral phase of a Gaussian pulse using the Fast Fourier Transform provided with numpy library in Python. Here are the results:
It is known that the spectral phase of a Fourier-limited Gaussian pulse should be flat (i.e. equal to some constant across the whole spectrum). However, after having Fourier-transformed the pulse, I am getting a strange spectral phase as a result. I find this phase using the function angle provided with numpy which should calculate the phase of a complex number. Is there any idea what I a missing here?
Here is the code that I use.
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import speed_of_light
from scipy.fftpack import fft, fftshift, ifft
# construct Electric field of a Gaussian pulse in time domain
def buildPulse(delays, wavelength = 700, fwhm_fs = 20, time_position = 0, phi0 = np.pi/2):
frequency = speed_of_light/wavelength*10**9 # frequency in Hz
period = 1/frequency*10**15 # period in fs
sigma = fwhm_fs*(2*np.sqrt(2*np.log(2)))**(-1)
envelope = np.exp(-(delays-time_position)**2/2/sigma**2)
Ef = envelope*np.sin(2*np.pi*frequency/10**15*delays-phi0) # electric field
return Ef, envelope, frequency, period, fwhm_fs
# fast Fourier Transform
def FFT(Ef,t0,tf,t_step,delays,fwhm_fs):
# get frequency axis in Hz and wavelength axis in nm
frequencies = np.linspace(10**(-25), (2*((tf-t0)/((tf-t0)/t_step+1)))**(-1), delays.size//2)*10**15
clbr_const = (frequencies[3]-frequencies[2])*10**15 # freq axis calibration const
print("clbr_const = ",clbr_const)
wavelengths = speed_of_light*frequencies**(-1)*10**(9) # in nm
wavelengths = wavelengths[::-1] # reversing wvelength axis to make it increasing
# fast Fourier transform of Ef (electric field in time domain)
spectrum = fft(Ef)
# half of the spectrum-vector
spectrum_half = spectrum[0:len(spectrum)//2]
# normalization factor
A = 2*(tf-t0)/spectrum.size/fwhm_fs
return spectrum,spectrum_half,frequencies,wavelengths,A
# construct time-axis (delays in fs)
t0 = -200 # first time-point in fs
tf = 200 # last time-point in fs
t_step = 0.01 # time-step in fs
delays = np.linspace(t0, tf, int((tf-t0)/t_step)+1)
# construct the electric field of a Gaussian pulse
Ef, envelope, frequency, period, fwhm_fs = buildPulse(\
delays, wavelength = 700, fwhm_fs = 20, time_position = 0, phi0 = np.pi/4)
print("frequency = ",frequency," Hz\nperiod = ",period,' fs')
# fast Fourier transform
spectrum,spectrum_half,frequencies,wavelengths,A = FFT(Ef,t0,tf,t_step,delays,fwhm_fs)
print("len(spectrum) = ",len(spectrum))
# amplitude spectrum
amp_s = A*np.abs(spectrum_half)
# real and imaginary parts of the spectrum
real_s = A*np.real(spectrum_half)
imag_s = np.imag(spectrum_half)
arg_s = np.angle(spectrum_half)
# plotting the results
plt.figure()
plt.plot(delays,Ef,delays,envelope)
plt.xlim([-100, 100])
plt.legend(["Electric field (Ef)","Envelope"],loc="upper left")
axs[0,0].set_xlabel('time (fs)')
plt.figure()
plt.plot(frequencies/10**12,amp_s)
plt.xlim([0, 1100])
plt.legend(["abs(fft(Ef))"],loc="upper left")
axs[1,2].set_xlabel('frequency (THz)')
plt.figure()
plt.plot(frequencies/10**12,arg_s)
plt.xlim([0, 1100])
plt.legend(["angle(fft(Ef))"])
plt.xlabel("frequency, THz")
