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-frequencies)*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")