Tutorials
Blog
About

Tutorials > Numerical Computing in Julia

FFT Derivative

by Martin D. Maas, Ph.D
@MartinDMaas

Last updated: 2021-12-03

Introduction

In many applications, we want to take advantage that the derivative operator is transformed into a multiplication operator in Fourier space. By applying a subsequent inverse transform, we can expect to obtain the derivative of a function using FFTs, using the well-known relation:

$$ \mathcal{F}[g’](\xi) = \frac{2\pi}{L} \xi \mathcal{F}[g](\xi) $$

However, this is not so simple to translate into the discrete setting.

Common Gotcha: Alising

For example, let’s compute the FFT-derivative of a periodic function defined in \( (0,L) \). The following code yields an incorrect result:

using FFTW
using plots

N = 20;
L = 1;
xj = (0:N-1)*L/N;
f = sin.(2π*xj)
df = 2π*cos.(2π*xj)

k = 1:N
df_fft = ifft( 2π/L * k.* fft(f) )

plot(xj,real(df),label="Exact derivative")
plot!(xj,real(df_fft),label="incorrect FFT derivative",markershape=:circle)

Wrong FFT-derivative

The reason this doesn’t work is related to aliasing, which we discussed in the previous section. In detail, what’s going on under the hood is that, even though in our grid, we have the equality

$$ \sin(5 x_j) = \frac{e^{-5 i x_j} + e^{5 i x_j}}{2i} + \frac{e^{15 i x_j} + e^{5 i x_j}}{2i}, $$

clearly, the frequency corresponding to k=15 should not be employed when computing the derivative.

The correct derivative algorithm uses the negative frequencies mentioned in the previous section.

Correct FFT-Derivative Algorithm

In order to obtain the correct result, we to correct these two lines:

k = fftfreq(N)*N;
df_fft = ifft( 2π*im/L * k.* fft(f) );

FFT-derivative

Don't miss any updates!

I'm writing a newsletter about once a month, with links to new posts and tutorials.

In particular, one of my goals for 2022 is to write a book about scientific software development using Julia. Make sure to subscribe if you are interested!