Why am I not getting a flat phase for a Gaussian pulse when doing a Fourier transform in Python? 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")
The group delay will be flat assuming the filter is linear phase (which means the phase is linear vs frequency). Group delay is the derivative of phase with respect to frequency. Either unwrap your phase and detrend it or use the group delay function that is part of scipy.signal to evaluate directly. Any FIR filter (not just the Gaussian) will be linear phase if (and only if) the filter coefficients are symmetric or anti-symmetric.