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:

enter image description here enter image description here enter image description here

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.xlim([-100, 100])
plt.legend(["Electric field (Ef)","Envelope"],loc="upper left")
axs[0,0].set_xlabel('time (fs)')
plt.xlim([0, 1100])
plt.legend(["abs(fft(Ef))"],loc="upper left")
axs[1,2].set_xlabel('frequency (THz)')    
plt.xlim([0, 1100])
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.

  • $\begingroup$ 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. $\endgroup$ Feb 17 '20 at 16:57
  • 1
    $\begingroup$ 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. $\endgroup$ Feb 17 '20 at 17:04
  • $\begingroup$ Could you please give an example of the usage of the group-delay function in my case? $\endgroup$ Feb 17 '20 at 21:47
  • $\begingroup$ 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/… $\endgroup$ Feb 17 '20 at 21:54
  • $\begingroup$ It would be scipy.signal.group_delay((b,a)) where b is your Gaussian pulse shape and a would be 1. $\endgroup$ 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.