# 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-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")


## 2 Answers

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.

• Thnks, after unwrapping the phase I am still not getting the flat phase. What do you mean by detrending the unwrapped phase? I will try the group delay function, but I wanted to understand how to do it with it first. – Capo Pavel Mestre Feb 17 '20 at 16:57
• The phase will not be flat is my point but could certainly be linear which would result in a “flat” group delay, meaning the delay at all frequencies is the same. Just use the group-delay function; detrending means subtracting our the linear phase but that isn’t the easier path for you to see this. – Dan Boschen Feb 17 '20 at 17:04
• Could you please give an example of the usage of the group-delay function in my case? – Capo Pavel Mestre Feb 17 '20 at 21:47
• See this post where the results for the matlab command grpdelay are shown (result is delay in samples): dsp.stackexchange.com/questions/63141/… For python see lagrange.univ-lyon1.fr/docs/scipy/0.17.1/generated/… – Dan Boschen Feb 17 '20 at 21:54
• It would be scipy.signal.group_delay((b,a)) where b is your Gaussian pulse shape and a would be 1. – Dan Boschen Feb 17 '20 at 21:55

Do an fftshift (Rotate the input vector by n/2) before the FFT to place the peak of the Gaussian at index 0.