computing dFT at the frequencies of the FFT

Question

I am computing the dFT of a function f(x) sampled at x_i, i=0,1,...,N (with a known dx) at the frequencies u_j, j=0,1,...,N where u_j are frequencies that np.fft.fftfreq(N, dx) generates and compare it with the result of np.fft.fft(f(x)). I find that the two do not agree...

Am I missing something? Shouldn't they by definition be the same? (The difference is even worse when I am looking at the imag parts of the dFT/FFT).

I am attaching the script that I used, which generates this plot that compares the real and imag parts of the dFT and FFT.

import numpy as np
import matplotlib.pyplot as plt
from astropy import units

def func_1D(x, sigma_x):
    return np.exp(-(x**2.0 / (2.0 * sigma_x**2)))


n_pixels = int(2**5.0)
pixel_scale = 0.05 # units of arcsec

x_rad = np.linspace(
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)


sigma_x = 0.5 # in units of arcsec
image = func_1D(
    x=x_rad,
    sigma_x=sigma_x * units.arcsec.to(units.rad),
)
image_FFT = np.fft.fftshift(np.fft.fft(np.fft.fftshift(image)))
u_grid = np.fft.fftshift(np.fft.fftfreq(n_pixels, d=pixel_scale * units.arcsec.to(units.rad)))

image_dFT = np.zeros(shape=n_pixels, dtype="complex")
for i in range(u_grid.shape[0]):
    for j in range(n_pixels):
        image_dFT[i] += image[j] * np.exp(
            -2.0
            * np.pi
            * 1j
            * (u_grid[i] * x_rad[j])
        )

value = 0.23

figure, axes = plt.subplots(nrows=1,ncols=3,figsize=(14,6))
axes[0].plot(x_rad * 10**6.0, image, marker="o")
for x_i in x_rad:
    axes[0].axvline(x_i * 10**6.0, linestyle="--", color="black")
axes[0].set_xlabel(r"x ($\times 10^6$; rad)")
axes[0].set_title("x-plane")

for u_grid_i in u_grid:
    axes[1].axvline(u_grid_i / 10**6.0, linestyle="--", color="black")
axes[1].plot(u_grid / 10**6.0, image_FFT.real, color="b")
axes[1].plot(u_grid / 10**6.0, image_dFT.real, color="r", linestyle="None", marker="o")
axes[1].set_title("u-plane (real)")
axes[1].set_xlabel(r"u ($\times 10^{-6}$; rad$^{-1}$)")
axes[1].plot(u_grid / 10**6.0, image_FFT.real - image_dFT.real, color="black", label="difference")

axes[2].plot(u_grid / 10**6.0, image_FFT.imag, color="b")
axes[2].plot(u_grid / 10**6.0, image_dFT.imag, color="r", linestyle="None", marker="o")
axes[2].set_title("u-plane (imag)")
axes[2].set_xlabel(r"u ($\times 10^{-6}$; rad$^{-1}$)")
#axes[2].plot(u_grid / 10**6.0, image_FFT.imag - image_dFT.imag, color="black", label="difference")
axes[1].legend()
plt.show()

If your input is real (I would allow the Gaussian to go to low values) then you should look for the absolute value in the spectrum. Please do read up on minimal example. I would leave out the units, some of these manipulations are not easy to follow. — roadrunner66, Feb 08 '20 at 07:49

score 1 · Accepted Answer · answered Feb 08 '20 at 08:09

I made a minimal example (I hope). I get essentially the same numbers for FFT and the naive Fourier integral (computed for the same frequency values).

import numpy as np
import matplotlib.pyplot as p 
%matplotlib inline

def signal(x, sigma_x):
    return np.exp(-(x**2.0 / (2.0 * sigma_x**2)))

t=np.linspace(-10,10,1000)
sigma=.3
sig=np.exp(-(t**2.0 / (2.0 * sigma **2)))

p.subplot(311)
p.plot(t,sig);

ft=np.fft.fftshift(np.fft.fft(sig))
freq=np.fft.fftshift(np.fft.fftfreq(1000,0.02))
p.subplot(312)
p.plot(freq,np.abs(ft))
print(np.abs(ft)[500:505])
# naive fourier integral 
fi=[]
for f in freq: 
  i=np.sum( sig* np.exp(- 1j* 2 *np.pi*f*t ))
  fi.append(np.abs(i))

p.subplot(313)
p.plot(freq,fi)

print(np.abs(fi)[500:505])

You need to plot both the real and imaginary components. Taking the magnitude hides problems caused by shifting of the signal (e.g. if the use of `fftshift` is incorrect). — Cris Luengo, Feb 08 '20 at 14:25
@roadrunner66 thank you for your example, it was very helpful. I updated it below to show instead the real and imaginary parts of the FFT. — Sketos, Feb 09 '20 at 15:50

score 0 · Answer 2 · answered Feb 09 '20 at 15:48

I updated the example by @roadrunner66 to show instead the real and imaginary parts of the FT rather than the magnitude, since the application I want to use it for involves dealing with the real and imaginary parts of the FT (commonly referred to as visibilities in interferometry).

Below is the slightly updated example.

import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(-10,10,1000)
sigma=.3
sig=np.exp(-(t**2.0 / (2.0 * sigma **2)))

ft=np.fft.fftshift(np.fft.fft(sig))
freq=np.fft.fftshift(np.fft.fftfreq(len(t),abs(t[0] - t[1])))

# naive fourier integral
fi_real=[]
fi_imag=[]
for f in freq:
  i=np.sum( sig* np.exp(- 1j* 2 *np.pi*f*t ))
  fi_real.append(i.real)
  fi_imag.append(i.imag)

figure, axes = plt.subplots(nrows=1,ncols=2)
axes[0].plot(freq,ft.real, color="b", label="np.fft.fft")
axes[0].plot(freq,fi_real, color="r", label="exact")
axes[0].set_xlim(-5.0, 5.0)
axes[0].set_title("real")
axes[0].legend()
axes[1].plot(freq,ft.imag, color="b", label="np.fft.fft")
axes[1].plot(freq,fi_imag, color="r", label="exact")
axes[1].set_xlim(-5.0, 5.0)
axes[1].set_title("imag")
axes[1].legend()
plt.show()

Looking at the output figure I think it's clear that the np.fft is not appropriate when you want to work with the real and imaginary parts of the FFT.

That is interesting. Could you provide a link? There are of course cases, where just the real or imaginary part are the interesting physical quantities (depends on the definition only), but I wasn't aware of it for the real or imaginary part of the electric field, where the applications I"m aware of are a mix of real and imaginary parts after any FFTs (field, power, FT(linear autocorrelation in time) = intensity spectrum , nonlinear autocorrelations in time, spatial fringe visibility <-> coherence length and many others). — roadrunner66, Feb 09 '20 at 20:29
The only thing you'd have to change from my example was to plot `real` and `imag` from `abs`. There is no need to generate purely real or purely imaginary vectors in numpy until you need to. The very idea of the complex FT is to avoid having to write everything twice. — roadrunner66, Feb 09 '20 at 20:33

computing dFT at the frequencies of the FFT

2 Answers2