Question about the FFT version of the gradient of a function.

fast fourier transformfourier transformnumerical methodsnumerical-calculus

We know that for a sufficiently smooth function $f:\mathbb{R}^{3}\to\mathbb{R}$, its Fourier Transform $\hat{f}(\mathbf{k}) \colon= \mathcal{F}\{f\}$ should satisfy (using integration by parts):

$$\mathcal{F}\{\boldsymbol{\nabla}f\}(\mathbf{k}) = \int\limits_{\mathbb{R}^3} e^{-2\pi i\mathbf{k}\cdot\mathbf{x}}\boldsymbol{\nabla}f(\mathbf{x})\mathrm{d}\mathbf{x} = \mathbf{k}\int\limits_{\mathbb{R}^3} e^{-2\pi i\mathbf{k}\cdot\mathbf{x}}f(\mathbf{x})\mathrm{d}\mathbf{x}\equiv\mathbf{k}\hat{f}(\mathbf{k})$$

so that differential equations can be transformed into algebraic equations through the Fourier Transform.

As you can see from the above you only need to multiply by the wave vector, and for a Laplacian it reduces to $\Delta\to k^{2}=\mathbf{k}\cdot\mathbf{k}$. What I do not understand is the way this is translated to the discrete version in the Fast Fourier Transform libraries for Python and Matlab.

In these programs, for a direction of the grid with with $n$ points and spacing $h$, the wave vector component is given by (the square braces in Python are used for arrays):

\begin{equation}
\left[0, 1, 2, \dots, \frac{n}{2}-1, -\frac{n}{2}, \dots, -2, -1\right]\div hn, \ \ \text{for $n$ even}
\label{1}\tag{1}
\end{equation}

\begin{equation}
\left[0, 1, 2, \dots, \frac{n-1}{2}, -\frac{n-1}{2}, \dots, -2, -1\right]\div hn, \ \ \text{for $n$ odd}
\label{2}\tag{2}
\end{equation}

What I see from equations $\eqref{1}$ and $\eqref{2}$ is that the array seems to have been "cut in half" and "displaced". But the FFT algorithm should give an $n$ element array with the wavenumbers in order (when you plot the Fourier Trasform, you do it right from the array without cutting and moving parts).

Am I missing something from this?

Best Answer

The FFT algorithm naturally produces an array which is "cut in half and displaced". The function FFT-Shift is usually applied at the end to fix this. The reason the shift isn't applied by default is that many programs which use the FFT don't require the shift and they will run more efficiently by not calling it.

I once heard a professor comment on this by saying "you don't mess with a $O(n\log(n))$ algorithm by unnecessarily sticking FFT-shift on it"

You can read more about this here : https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fftshift.html

Related Question