# Numpy FFT gives me a pulse shorter than it should be. Not sure what I am doing wrong

I've created a code (Python, numpy) that defines an ultrashort laser pulse in the frequency domain (pulse duration should be 4 fs), but when I perform the Fourier Transform using DFT, my pulse in the time domain is actually shorter than it should be.

Here's my code:

import numpy as np
import matplotlib.pyplot as plt

#pulse duration [fs] FWHM -- what I should get after the FT
T0  = 4

#speed of light [nm/fs]
c = 299792458*10**(-6)
#central wavelength [nm]
wl0 = 800
w0 = (2*np.pi*c)/wl0
delta_w = 2*np.pi*0.441/T0
#bandwidth in [nm]
delta_wl = delta_w*wl0**2/(2*np.pi*c)

w = np.linspace(w0-delta_w*8, w0+delta_w*8, 2**9)
#wavelengths [nm]
wl = c/w
#frequencies [PHz]
f = w/(2*np.pi)

#to make the spectrum centered around the carrier frequency
diff_w = w-w0
sigma_w = delta_w/(np.sqrt(8*np.log(2)))
spectrum_w = np.exp(-(diff_w**2)/(2*sigma_w**2))

#phase terms (not relevant here)
phi_w = 0
def phase(phi_0,phi_1,phi_2,phi_3):
phi_w = phi_0 + phi_1*(w-w0) + (phi_2*(w-w0)**2)/math.factorial(2) + (phi_3*(w-w0)**3)/math.factorial(3)
return phi_w

#field in the frequency domain
E_w = np.exp(1j*phase(0,0,0,0))*spectrum_w

#FT:
n = len(E_w)
timestep = 0.01
fa = 1.0/timestep

t_1 = np.fft.fftfreq(n,d = timestep)
t = np.fft.fftshift(t_1)

field_ft = np.fft.ifft(E_w)
plt.plot(t,field_ft)
plt.show()

new = np.fft.fftshift(field_ft)
plt.xlim(-5,5)
plt.plot(t,new,t,np.abs(new))
plt.show()


The output I get is a pulse that is shorter than it should be. I've found a pretty similar question posted here almost 4 years ago, but it received no responses. I am being hopeful today!

As the previous poster said, this should be a pretty straightforward and simple code, but it hasn't been that simple due to this issue I keep running into.

• Have you tried adding the negative part of your spectrum? – nicoguaro Jun 28 '19 at 14:32
• I think I see the problem, and after looking at the fftshift code here github.com/numpy/numpy/blob/… I think I see the solution. If you post the code that shows what the expected graphs before and after are supposed to look like, I'll dig in a little deeper. – mathew gunther Jun 28 '19 at 14:47

Running your code, it seems like your pulse looks kinda like this:

(sorry for not adding units to the plots, I used the same as you, i.e. t is in fs and w in rad/fs)

So, the FWHM is not correct (should be 4 fs) and the angular carrier/central frequency is messed up (should be ca. 2.3 / ca. 0.37 per fs). I think there is two things not quite right here:

1) the definition of your w-axis

2) the conversion of frequency to time

As by definition of the FFT, where there is already a $$2\pi$$ in the exponents, we get $$\sigma_t = \frac{1}{2\sigma_w}$$ where $$\sigma_t$$ and $$\sigma_w$$ are the standard deviations of the intensity in t- and w-space respectively.

Now, according to wikipedia and with $$\Delta t$$ and $$\Delta w$$ as the FWHM in t- and w- space respectively $$\text{fac} \cdot \sigma_t = \Delta t$$ and respectively $$\text{fac} \cdot \sigma_w = \Delta w$$ with $$\text{fac} = \sqrt{8\cdot \text{log}2}$$. Plugging all these in you get $$\Delta w = \text{fac}^2 \cdot \frac{1}{2\Delta t} \approx 2\pi \cdot 0.441 \cdot \frac{1}{\Delta t}$$ as you did correctly. However, in laser physics, people use the FWHM of the intensity of the field, rather than the FWHM of the amplitude. This is why we use $$\sigma_t \sigma_w = \frac{1}{2}$$. The $$\frac{1}{2}$$ originates from the fact that the intensity is the amplitude squared. This is a small but important detail, as we have to pay attention to put factors of $$2$$ and $$\sqrt{2}$$ in the appropriate places.

Now for 1), the FFT-implementation is going to consider the spectrum to start at w=0. The way you defined the w-axis, what the algorithm saw in E_w, was a Gaussian centered half the way up your w-axis, because you had the axis symmetrically spread around your peak. That means, after the FFT, your central frequency was significantly higher than it should have been. This is fixed by letting w start at zero and go to some high enough value (you chose 16 bandwidths, so I stuck with that too):

#angular frequencies [rad/fs]
w = np.linspace(0, delta_w*16, 2**9)


As above, for the Gaussian you need the $$\sigma_w$$ from the FWHM. Don't forget a factor $$\sqrt{2}$$ because of the whole intensity-amplitude-issue:

sigma_w = delta_w/fac # standard deviation of intensity
sigma_w_field = np.sqrt(2)*sigma_w # standard deviation of amplitude

#to make the spectrum centered around the carrier frequency
diff_w = w-w0


I removed the phase function, as it didn't seem to do anything here and also make sure to use sigma_w_field:

spectrum_w = np.exp(-(diff_w**2)/(2*sigma_w_field**2))

#field in the frequency domain
E_w = spectrum_w
plt.figure()
plt.xlabel("w")
plt.ylabel("E_w")
plt.plot(w, E_w)
plt.grid()
plt.show()


Technically, the above plot is not 100% correct, because the left shoulder of the Gaussian might reappear in the high end of the spectrum. This is not an issue here, because at w=0 the Gaussian has decayed enough, but if you choose a wider band width or a lower central frequency, you should take care of that somehow (I have no idea how to do this elegantly though). This fixes Item 1).

For Item 2), let's look at the variable timestep. timestep is supposed to be the inverse of the sampling frequency, however, I don't think 0.01 was the correct value for that. The inverse of the sampling frequency, i.e. the sampling time, is the length of your signal in f-space divided by the number of samples. Here this is 16 bandwidths over $$2^9$$ (here called w_s). This is however in w-space and fftfreq wants f-space, so one division by $$2\pi$$ is requiered. And again, delta_w was computed for the intensity, so we need another $$\sqrt{2}$$ to turn it into the band-width of the spectrum:

#FT:
n = len(E_w)
w_s =  np.sqrt(2)*16*delta_w/2**9
timestep = w_s/(2*np.pi)
#fa = 1.0/timestep

t_1 = np.fft.fftfreq(n,d = timestep)
t = np.fft.fftshift(t_1)

field_ft = np.fft.ifft(E_w)

new = np.fft.fftshift(field_ft)
plt.figure()
plt.xlabel("t")
plt.ylabel("field")
plt.xlim(-10,10)
plt.plot(t,new, label="new")
plt.plot(t,np.abs(new), label="abs(new)")
plt.plot(t, np.ones(len(new))*0.5*np.max(np.abs(new)), label="half maximum")
plt.grid()
plt.legend()
plt.show()


The carrier-wave-plot can be improved by making the range of w larger. I hope this is the pulse you were looking for ;)

• This is exactly what I was looking for! This was a very clear explanation, thank you! – python_enthusiast Jun 29 '19 at 15:32

I don't think the sqrt(2) is necessary in: w_s = np.sqrt(2)*16*delta_w/2**9 When it is removed, the correct fwhm pulse duration for the intensity is obtained.