Fourier Transform Using Numpy

Question

I'm trying to calculate the Fourier Transform of the following Gaussian:

# sample spacing
dx = 1.0 / 1000.0

# Points
x1 = -5
x2 = 5

x = np.arange(x1, x2, dx)

def light_intensity():
    return 10*sp.stats.norm.pdf(x, 0, 1)+0.1*np.random.randn(x.size)

fig, ax = plt.subplots()
ax.plot(x,light_intensity())

I create a new array in the spacial frequency domain (Fourier Transform of Gaussian is a Gaussian so these values should be similar). I plot and get this:

fig, ax = plt.subplots()

xf = np.arange(x1,x2,dx)
yf= np.fft.fftshift(light_intensity())
ax.plot(xf,np.abs(yf))

Why is it splitting into two peaks?

OK, this took me a while, but basically, your x labels should start at 0. — Mateen Ulhaq, Sep 16 '18 at 04:13
[`fftshift`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftshift.html) does not perform a Fourier transform. To actually compute the FFT, use the function [`numpy.fft.fft`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html). — Warren Weckesser, Sep 16 '18 at 04:18
Using fft.fft gives an even more bizarre plot when I substitute it into the code above. I get two delta function peaks at each side. — Luke Polson, Sep 16 '18 at 04:27

Mateen Ulhaq · Accepted Answer · 2020-07-26T00:01:57.883

Advice:

use np.fft.fft
fft starts at 0 Hz
normalize/rescale

Complete example:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def norm_fft(y, T, max_freq=None):
    N = y.shape[0]
    Nf = N // 2 if max_freq is None else int(max_freq * T)
    xf = np.linspace(0.0, 0.5 * N / T, N // 2)
    yf = 2.0 / N * np.fft.fft(y)
    return xf[:Nf], yf[:Nf]

def generate_signal(x, signal_gain=10.0, noise_gain=0.0):
    signal = norm.pdf(x, 0, 1)
    noise = np.random.randn(x.size)
    return signal_gain * signal + noise_gain * noise

# Signal parameters
x1 = 0.0
x2 = 5.0
N = 10000
T = x2 - x1

# Generate signal data
x = np.linspace(x1, x2, N)
y = generate_signal(x)

# Apply FFT
xf, yf = norm_fft(y, T, T / np.pi)

# Plot
fig, ax = plt.subplots(2)
ax[0].plot(x, y)
ax[1].plot(xf, np.abs(yf))
plt.show()

Or, with noise:

Plots with symmetry

Alternatively, if you want to enjoy the symmetry in the frequency domain:

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

def norm_sym_fft(y, T, max_freq=None):
    N = y.shape[0]
    b = N if max_freq is None else int(max_freq * T + N // 2)
    a = N - b
    xf = np.linspace(-0.5 * N / T, 0.5 * N / T, N)
    yf = 2.0 / N * np.fft.fftshift(np.fft.fft(y))
    return xf[a:b], yf[a:b]

def generate_signal(x, signal_gain=10.0, noise_gain=0.0):
    signal = norm.pdf(x, 0, 1)
    noise = np.random.randn(x.size)
    return signal_gain * signal + noise_gain * noise

# Signal parameters
x1 = -10.0
x2 = 10.0
N = 10000
T = x2 - x1

# Generate signal data
x = np.linspace(x1, x2, N)
y = generate_signal(x)

# Apply FFT
xf, yf = norm_sym_fft(y, T, 4 / np.pi)

# Plot
fig, ax = plt.subplots(2)
ax[0].plot(x, y)
ax[1].plot(xf, np.abs(yf))
plt.show()

Or, with noise:

Quick question. The Fourier Transform of a Gaussian is supposed to be the same Gaussian function itself right? Take for example: http://www.wolframalpha.com/input/?i=fourier+transform+of+4e%5E(-x%5E2%2F2) where the Fourier transformed function is the same as the original function (in the frequency domain). Is there a reason the height of the fft function and the standard deviation are different in the plots you showed? — Luke Polson, Sep 16 '18 at 06:21
I believe it's because the FFT=DFT, which is different from a FT. — Mateen Ulhaq, Sep 16 '18 at 09:02

score 3 · Answer 2 · answered Sep 16 '18 at 05:26

First, use np.fft.fft to computes the Fourier Transform then use np.fft.fftshift to shift the zero-frequency component to the center of the spectrum.

Replace the second part of your code with:

xf = np.arange(x1,x2,dx)
yf = np.fft.fft(light_intensity())
yfft = np.fft.fftshift(np.abs(yf))
fig,ax = plt.subplots(1,2,figsize=(10,5))
ax[0].plot(xf,light_intensity())
ax[1].plot(xf,yfft)
ax[1].set_xlim(-0.05,0.05)
plt.show()

This is the result:

Fourier Transform Using Numpy

2 Answers2

Plots with symmetry