3.4. Fourier Transform¶

The 1D Fourier transform is:

$F[f(x)] \equiv \tilde f(\omega) = \int_{-\infty}^{\infty} f(x) e^{-i\omega x}\,\d x F^{-1}[\tilde f(\omega)] = f(x) = {1\over2\pi}\int_{-\infty}^{\infty} \tilde f(\omega) e^{+i\omega x}\,\d \omega$

To show that it works:

$F^{-1} F [f(x)] = {1\over2\pi}\int_{-\infty}^{\infty} \left[\int_{-\infty}^{\infty} f(x) e^{-i\omega x}\,\d x\right] e^{+i\omega x}\,\d \omega = {1\over2\pi}\int_{-\infty}^{\infty} \left[\int_{-\infty}^{\infty} f(x') e^{-i\omega x'}\,\d x'\right] e^{+i\omega x}\,\d \omega = = \int_{-\infty}^{\infty} f(x') \left[{1\over2\pi}\int_{-\infty}^{\infty} e^{i\omega (x- x')}\,\d \omega \right] \,\d x' = \int_{-\infty}^{\infty} f(x') \delta(x-x') \,\d x' =f(x)$

If $x$ is time (unit $[\mathrm{s}]$ ), then $\omega$ is angular frequency (unit $[\mathrm{rad}/\mathrm{s}]$ ). One can express the Fourier transform in terms of ordinary frequency $\nu$ (unit $[1/\mathrm{s}] = [\mathrm{Hz}]$ ) by substituting $\omega = 2\pi \nu$ :

$\tilde f(\omega) = \tilde f(2\pi \nu) \equiv \tilde f'(\nu) = \int_{-\infty}^{\infty} f(x) e^{-2\pi i\nu x}\,\d x f(x) = \int_{-\infty}^{\infty} \tilde f'(\nu) e^{+2\pi i\nu x}\,\d \nu$

Both transformations are equivalent and only differ in whether we express the transform in terms of $\omega$ or $\nu$ , the conversion being given by $\tilde f(\omega) = \tilde f(2\pi \nu) \equiv \tilde f'(\nu)$ . Third frequently used convention that is however not equivalent to the above is:

$\tilde f(k) = {1\over\sqrt{2\pi}} \int_{-\infty}^{\infty} f(x) e^{-ik x}\,\d x f(x) = {1\over\sqrt{2\pi}} \int_{-\infty}^{\infty} \tilde f(k) e^{+ik x}\,\d k$

The 3D Fourier transform is:

(3.4.1)¶ $F[f(\mathbf{x})] \equiv \tilde f(\bomega) = \int_{-\infty}^{\infty} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x F^{-1}[\tilde f(\bomega)] = f(\mathbf{x}) = {1\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega$

With obvious analogs for other conventions and dimensions.

The sign convention in the exponentials $e^{\pm i\omega x}$ is arbitrary, one can as well flip the sign of the direct and inverse transforms. In particular, one often uses both sign conventions in the same equation. Consider a spacetime plane-wave $e^{i k\cdot x} = e^{i(\omega t - \mathbf{k}\cdot\mathbf{x})}$ . Then we obtain (using plus sign convention in the $e^{ik \cdot x}$ exponential for the direct transformation):

$F[f(x)] \equiv \tilde f(k) = \int_{-\infty}^{\infty} f(x) e^{ik \cdot x}\,\d^4 x = \int_{-\infty}^{\infty} f(x) e^{i(\omega t - \mathbf{k}\cdot\mathbf{x})}\,\d^4 x F^{-1}[f(k)] \equiv f(x) = {1\over(2\pi)^4} \int_{-\infty}^{\infty} \tilde f(k) e^{-ik \cdot x}\,\d^4 k = {1\over(2\pi)^4} \int_{-\infty}^{\infty} \tilde f(k) e^{-i(\omega t - \mathbf{k}\cdot\mathbf{x})}\,\d^4 k$

Finally, the equation $k\cdot x = \omega t - \mathbf{k}\cdot\mathbf{x}$ depends on the metric signature, in this case $\diag({1, -1, -1, -1})$ . For a signature $\diag({-1, 1, 1, 1})$ we would get $k\cdot x = -\omega t + \mathbf{k}\cdot\mathbf{x}$ .

Unlike the normalization convention, where one has to be very careful, the sign convention in Fourier transform is not a problem, one just has to remember to flip the sign for the inverse transform.

3.4.1. Shift Theorem¶

The Fourier transform of a shifted function, in 3D:

$F[f(\mathbf{x}+\mathbf{b})] = \int_{-\infty}^{\infty} f(\mathbf{x}+\mathbf{b}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = \int_{-\infty}^{\infty} f(\mathbf{x}) e^{-i\bomega \cdot (\mathbf{x}-\mathbf{b})}\,\d^3 x = = e^{i\bomega\cdot \mathbf{b}} \int_{-\infty}^{\infty} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = e^{i\bomega\cdot \mathbf{b}} F[f(\mathbf{x})]$

3.4.2. Scaling¶

For $a>0$ :

$F[f(ax)](\omega) = \int_{-\infty}^{\infty} f(ax) e^{-i\omega x}\,\d x = = {1\over a}\int_{-\infty}^{\infty} f(y) e^{-i{\omega\over a} y}\,\d y = = {1\over a}F[f(x)]\left({\omega\over a}\right)$

3.4.3. Derivative¶

The Fourier transform of a derivative, in 3D:

$F[\partial_i f(\mathbf{x})] = \int_{-\infty}^{\infty} (\partial_i f(\mathbf{x})) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = \left[f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\right]_{-\infty}^{\infty} -\int_{-\infty}^{\infty} f(\mathbf{x}) \partial_i e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = -\int_{-\infty}^{\infty} f(\mathbf{x}) \partial_i e^{-i\omega_j x^j}\,\d^3 x = = -(-i\omega_i)\int_{-\infty}^{\infty} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = i\omega_i F[f(\mathbf{x})]\,.$

An alternative derivation is to start from:

$f(\mathbf{x}) = F^{-1}[\tilde f(\bomega)] = {1\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega$

and differentiate both sides:

$\partial_i f(\mathbf{x}) = {1\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\bomega) \partial_i e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega \partial_i f(\mathbf{x}) = {1\over(2\pi)^3}\int_{-\infty}^{\infty} i\omega_i \tilde f(\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega \,,$

from which:

$F[\partial_i f(\mathbf{x})] = i\omega_i \tilde f(\bomega) = i\omega_i F[f(\mathbf{x})]\,.$

3.4.4. Convolution¶

The convolution of two functions $f(x)$ and $g(x)$ is defined as:

$f(x) * g(x) = \int_{-\infty}^\infty f(y) g(x-y) \,\d y$

The Fourier transform of a convolution is:

$F[f(x) * g(x)](\omega) = \int_{-\infty}^\infty \int_{-\infty}^\infty f(y) g(x-y) \,\d y\, e^{-i\omega x}\d x = = \int_{-\infty}^\infty \int_{-\infty}^\infty g(x-y) e^{-i\omega x} \d x\, f(y) \,\d y = = \int_{-\infty}^\infty \int_{-\infty}^\infty g(u) e^{-i\omega (u+y)} \d u\, f(y) \,\d y = = \int_{-\infty}^\infty g(u) e^{-i\omega u} \d u \int_{-\infty}^\infty f(y) e^{-i\omega y} \d y = F[f(x)](\omega)\ F[g(x)](\omega)$

And for the inverse transform:

$F^{-1}[f(\omega) * g(\omega)](x) = {1\over 2\pi} \int_{-\infty}^\infty \int_{-\infty}^\infty f(y) g(\omega-y) \,\d y\, e^{i\omega x}\d \omega = = {1\over 2\pi} \int_{-\infty}^\infty \int_{-\infty}^\infty g(\omega-y) e^{i\omega x} \d \omega\, f(y) \,\d y = = {1\over 2\pi} \int_{-\infty}^\infty \int_{-\infty}^\infty g(u) e^{ix (u+y)} \d u\, f(y) \,\d y = = 2\pi {1\over 2\pi} \int_{-\infty}^\infty g(u) e^{ix u} \d u {1\over 2\pi} \int_{-\infty}^\infty f(y) e^{ix y} \d y = 2\pi F^{-1}[f(\omega)](x)\ F^{-1}[g(\omega)](x)$

Fourier transform of a function multiplication is:

$F [ f g ] = F [\ F^{-1}[ F[f] ]\quad F^{-1}[ F[g] ]\ ] = {1\over 2\pi} F [ F^{-1} [ F[f] * F[g] ]] = {1\over 2\pi} F[f] * F[g]$

and for the inverse transform:

$F^{-1} [ f g ] = F^{-1} [\ F[ F^{-1}[f] ]\quad F[ F^{-1}[g] ]\ ] = F^{-1} [ F [ F^{-1}[f] * F^{-1}[g] ]] = F^{-1}[f] * F^{-1}[g]$

3.4.5. Radial Fourier Transform¶

As a special case when the function $f(\mathbf{x})=f(r)$ is spherically symmetric, we introduce spherical coordinates such that the $z$ -axis is along the $\bomega$ vector and calculate (we use $r=|\mathbf{x}|$ and $\omega=|\bomega|$ ):

$F[f(\mathbf{x})] \equiv \tilde f(\bomega) = \int_{-\infty}^{\infty} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = \int_{-\infty}^{\infty} f(r) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = \int_0^\infty\d r \int_0^\pi\d\theta \int_0^{2\pi}\d\phi f(r) e^{-i \omega r \cos\theta}\,r^2\sin\theta = = 2\pi \int_0^\infty\d r \int_0^\pi\d\theta f(r) e^{-i \omega r \cos\theta}\,r^2\sin\theta = = 4\pi \int_0^\infty f(r) \sinc(\omega r) \,r^2 \d r = = 4\pi \int_0^\infty f(r) {\sin\omega r \over \omega r}\,r^2 \d r = = {4\pi\over\omega} \int_0^\infty r\sin(\omega r) f(r) \,\d r\,,$

where we used:

$\int_0^\pi e^{-i \omega r \cos\theta}\,\sin\theta \d\theta = \int_{-1}^1 e^{i\omega r u} \d u = \left[e^{i\omega r u} \over i\omega r\right]_{-1}^1 = {e^{i\omega r} - e^{-i\omega r} \over i \omega r} = = 2 {\sin(\omega r) \over \omega r} = 2 \sinc(\omega r) = 2 j_0(\omega r)\,.$

So the transform is real and spherically symmetric, since the result only depends on $\omega$ .

Similarly, for the inverse transform:

$F^{-1}[\tilde f(\bomega)] = f(\mathbf{x}) = {1\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega = = {1\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\omega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega = = {1\over(2\pi)^3} {4\pi\over r} \int_0^\infty \omega\sin(\omega r) f(\omega) \,\d \omega = = {1\over 2\pi^2 r} \int_0^\infty \omega\sin(\omega r) f(\omega) \,\d \omega$

3.4.6. Examples¶

Rectangular Function¶

The rectangular function is defined as:

$\Pi(x) = H(x+\half) - H(x-\half)$

The Fourier transform is:

$F[\Pi(x)] \equiv \tilde \Pi(\omega) = \int_{-\infty}^{\infty} \Pi(x) e^{-i\omega x}\,\d x = \int_{-\half}^{\half} e^{-i\omega x}\,\d x = = \left[e^{-i\omega x} \over -i\omega x\right]_{x=-\half}^\half = {e^{i{\omega\over 2}} - e^{-i{\omega \over 2}} \over i {\omega\over2}} = 2{\sin({\omega \over 2}) \over {\omega \over 2}} = 2\sinc\left({\omega\over 2}\right)\,.$

Dirichlet Kernel¶

The Dirichlet kernel $D_N(x)$ is a partial sum of complex exponentials:

$D_N(x) = {1\over 2\pi}\sum_{n=-N}^N e^{inx} = = {1\over 2\pi}\left(1+2\sum_{n=1}^N \cos(nx)\right) = = {1\over 2\pi \sin\left(x\over2\right)}\left(\sin\left(x\over2\right) +2\sum_{n=1}^N \cos(nx)\sin\left(x\over2\right)\right) = = {1\over 2\pi \sin\left(x\over2\right)}\left(\sin\left(x\over2\right) +\sum_{n=1}^N\left( \sin\left(\left(n+\half\right)x\right) -\sin\left(\left(n-\half\right)x\right) \right)\right) = = {\sin\left(\left(N+\half\right)x\right) \over 2\pi \sin\left(x\over2\right)}$

From the definition, it is a periodic function with period $2\pi$ .

Integral of it is equal to one:

$\int_{-\pi}^\pi D_N(x) \d x = \int_{-\pi}^\pi {1\over 2\pi}\left(1+2\sum_{n=1}^N \cos(nx)\right) \d x = = 1 + {1\over \pi}\sum_{n=1}^N \int_{-\pi}^\pi \cos(nx)\, \d x = 1$

also

$\int_{-\pi}^\pi D_N(x-y) \d y = 1$

The Dirichlet kernel $D_N(x)$ converges towards a train of delta functions (called Dirac comb, see the equation (3.4.6.2) in the next section):

(3.4.6.1)¶ ${1\over 2\pi}\sum_{n=-\infty}^\infty e^{inx} = \lim_{N\to\infty} {1\over 2\pi}\sum_{n=-N}^N e^{inx} = \lim_{N\to\infty} D_N(x) = = \lim_{N\to\infty} {\sin\left(\left(N+\half\right)x\right) \over 2\pi \sin\left(x\over2\right)} = \sum_{n=-\infty}^\infty \delta(x-2\pi n)$

Let us do the crucial step in more details using distributions:

$\int_{-\infty}^\infty \lim_{N\to\infty} {\sin\left(\left(N+\half\right)x\right) \over 2\pi \sin\left(x\over2\right)} \varphi(x) \,\d x = = \sum_{n=-\infty}^\infty \lim_{N\to\infty} \int_{-\pi}^\pi {\sin\left(\left(N+\half\right)(x+2\pi n)\right) \over 2\pi \sin\left(x+2\pi n\over2\right)} \varphi(x+2\pi n) \,\d x = = \sum_{n=-\infty}^\infty \varphi(2\pi n) = =\int_{-\infty}^\infty \sum_{n=-\infty}^\infty \delta(x-2\pi n) \varphi(x)\,\d x$

Where we used the fact that

$\left[\lim_{N\to\infty} \int_{-\pi}^\pi {\sin\left(\left(N+\half\right)(x+2\pi n)\right) \over 2\pi \sin\left(x+2\pi n\over2\right)} \varphi(x+2\pi n) \,\d x \right] - \varphi(2\pi n) = = \left[\lim_{N\to\infty} \int_{-\pi}^\pi D_N(x+2\pi n) \varphi(x+2\pi n) \,\d x \right] - \varphi(2\pi n) = = \lim_{N\to\infty} \int_{-\pi}^\pi D_N(x+2\pi n) \left(\varphi(x+2\pi n)-\varphi(2\pi n)\right) \,\d x = = \lim_{N\to\infty} \int_{-\pi}^\pi {\varphi(x+2\pi n)-\varphi(2\pi n) \over 2\pi\sin\left(x+2\pi n\over 2\right)} \sin\left(\left(N+\half\right)(x+2\pi n)\right) \,\d x = = 0$

Dirac Comb (Shah) Function¶

The Dirac comb function, also called the Shah function, is defined as:

$\Sh(x) = \sum_{n=-\infty}^\infty \delta(x-n)$

It has the following scaling property:

$\Sh(ax) = \sum_{n=-\infty}^\infty \delta(ax-n) = \sum_{n=-\infty}^\infty \delta\left(a\left(x-{n\over a}\right)\right) = \sum_{n=-\infty}^\infty {1\over|a|}\delta\left(x-{n\over a}\right)$

and for $a={1\over L}$ with $L>0$ :

$\Sh\left({x\over L}\right) = \sum_{n=-\infty}^\infty L\delta(x-nL)$

From which a train of delta functions $L$ distance apart is expressed using a Dirac comb as:

$\sum_{n=-\infty}^\infty \delta(x-nL) = {1\over L}\Sh\left({x\over L}\right)$

Using the identity (3.4.6.1), the infinite sum of complex exponentials is also equal to a Dirac comb:

(3.4.6.2)¶ ${1\over 2\pi}\sum_{n=-\infty}^\infty e^{inx} = \lim_{N\to\infty} {1\over 2\pi}\sum_{n=-N}^N e^{inx} = \lim_{N\to\infty} D_N(x) = = \lim_{N\to\infty} {\sin\left(\left(N+\half\right)x\right) \over 2\pi \sin\left(x\over2\right)} = \sum_{n=-\infty}^\infty \delta(x-2\pi n) = {1\over 2\pi}\Sh\left({x\over 2\pi}\right)$

Using (3.4.6.2) we can now calculate the Fourier transform:

$F[\Sh(x)](\omega) \equiv \tilde \Sh(\omega) = \int_{-\infty}^{\infty} \Sh(x) e^{-i\omega x}\,\d x = = \int_{-\infty}^{\infty} \sum_{n=-\infty}^\infty \delta(x-n) e^{-i\omega x}\,\d x = = \sum_{n=-\infty}^\infty \int_{-\infty}^{\infty} \delta(x-n) e^{-i\omega x}\,\d x = = \sum_{n=-\infty}^\infty e^{-i\omega n} = = 2\pi \sum_{n=-\infty}^\infty \delta(\omega-2\pi n) = = \Sh\left({\omega\over 2\pi}\right)$

For the inverse Fourier transform we get (using the previous result):

$F^{-1}[\Sh(\omega)](x) = F^{-1}\left[\ \Sh\left({(2\pi\omega)\over2\pi}\right) \ \right](x) = F^{-1}[\ F[\Sh(x)](2\pi\omega) \ ](x) = F^{-1}\left[\ F\left[{1\over2\pi}\Sh\left({x\over2\pi}\right)\right](\omega) \ \right](x) = {1\over2\pi}\Sh\left({x\over2\pi}\right)$

The following Fourier transform is also useful:

$F\left[\sum_{n=-\infty}^\infty \delta(x-nL)\right](\omega) = F\left[{1\over L}\Sh\left({x\over L}\right)\right](\omega) = F[\Sh(x)](L\omega) = = \Sh\left({L\over 2\pi}\omega\right) = \sum_{n=-\infty}^\infty {2\pi\over L} \delta\left(x-{2\pi n\over L}\right)$

Periodic Summation¶

The convolution $f(x) * g(x) = \int_{-\infty}^\infty f(y) g(x-y) \,\d y$ of a Dirac comb $\Sh(x)$ and an arbitrary function $f(x)$ is called a periodic summation:

$f(x) * \Sh(x) = \int_{-\infty}^\infty f(y) \Sh(x-y) \,\d y = \int_{-\infty}^\infty f(y) \sum_{n=-\infty}^\infty \delta(x-y-n) \,\d y = = \sum_{n=-\infty}^\infty f(x-n) = \sum_{n=-\infty}^\infty f(x+n)$

because the result is a periodic function with period 1:

$(f * \Sh)(x+1) = \sum_{n=-\infty}^\infty f(x+n+1) = \sum_{m=-\infty}^\infty f(x+m) = (f * \Sh)(x)$

Poisson Summation Formula¶

The Poisson summation formula:

(3.4.6.3)¶ $\sum_{n=-\infty}^\infty f(2\pi n) = {1\over 2\pi} \sum_{n=-\infty}^\infty \tilde f(n)$

can be derived using a Dirac comb:

$\sum_{n=-\infty}^\infty f(2\pi n) = \int_{-\infty}^\infty f(x) \sum_{n=-\infty}^\infty \delta(x-2\pi n) \,\d x = = \int_{-\infty}^\infty f(x) {1\over 2\pi} \Sh\left(x\over2\pi\right) \,\d x = = {1\over 2\pi} \int_{-\infty}^\infty f(x) \cdot F[\Sh(\omega)](x) \,\d x = = {1\over 2\pi} \int_{-\infty}^\infty F[f(x)](\omega) \cdot \Sh(\omega) \,\d \omega = = {1\over 2\pi} \int_{-\infty}^\infty \tilde f(\omega) \cdot \sum_{n=-\infty}^\infty \delta(\omega-n) \,\d \omega = = {1\over 2\pi} \sum_{n=-\infty}^\infty \tilde f(x)$

An alternative derivation using Fourier series (see next sections):

$\sum_{n=-\infty}^\infty f(x+2\pi n) = g(x) = \sum_{n=-\infty}^\infty {1\over 2\pi}\int_{-\pi}^\pi g(y) e^{-iny} \d y\, e^{inx} = = \sum_{n=-\infty}^\infty {1\over 2\pi}\int_{-\pi}^\pi \sum_{m=-\infty}^\infty f(y+2\pi m) e^{-iny} \d y\, e^{inx} = = \sum_{n=-\infty}^\infty {1\over 2\pi} \sum_{m=-\infty}^\infty \int_{-\pi}^\pi f(y+2\pi m) e^{-in(y+2\pi m)} \d y\, e^{inx} = = \sum_{n=-\infty}^\infty {1\over 2\pi} \int_{-\infty}^\infty f(y) e^{-iny} \d y\, e^{inx} = = \sum_{n=-\infty}^\infty {1\over 2\pi} \tilde f(n) \, e^{inx}$

And setting $x=0$ we get the Poisson summation formula (3.4.6.3).

The last derivation can actually also be done using a Dirac comb function as follows:

$\sum_{n=-\infty}^\infty f(x+2\pi n) = f(x) * {1\over2\pi}\Sh\left(x \over 2\pi \right ) = = F^{-1}\left[\ F\left[f(x) * {1\over2\pi}\Sh\left(x \over 2\pi \right ) \right](\omega)\ \right](x) = = F^{-1}\left[\ F[f(x)](\omega)\ F\left[{1\over2\pi}\Sh\left(x \over 2\pi \right )\right](\omega)\ \right](x) = = F^{-1}[\ F[f(x)](\omega)\ \Sh(\omega)\ ](x) = = {1\over 2\pi} \int_{-\infty}^\infty F[f(x)](\omega) \Sh(\omega) e^{i\omega x} \d x = = \sum_{n=-\infty}^\infty {1\over 2\pi} F[f(x)](n) \, e^{inx}$

Fourier Series¶

Consider a periodic function $f(x)$ with a period $L$ and let us calculate the Fourier transform of it. We define a new function $f_0(x)=f(x)$ in the $[0, L]$ interval and zero otherwise. Then:

$f(x) = f_0(x) * {1\over L}\Sh\left(x\over L\right) = \sum_{n=-\infty}^\infty f_0(x+Ln)$

Apply Fourier transform:

(3.4.6.4)¶ $F[f(x)](\omega) = F\left[f_0(x) * {1\over L}\Sh\left(x\over L\right)\right](\omega) = = F[f_0(x)](\omega)\ F\left[{1\over L}\Sh\left(x\over L\right)\ \right](\omega) = F[f_0(x)](\omega)\ \Sh\left({L\over2\pi}\omega\right) = = \sum_{n=-\infty}^\infty F[f_0(x)]\left({2\pi n\over L}\right) \ {2\pi\over L}\delta\left(\omega-{2\pi n\over L}\right)= = \sum_{n=-\infty}^\infty \int_0^L f(x) e^{-i 2\pi n x / L} \d x \ {2\pi\over L}\delta\left(\omega-{2\pi n\over L}\right)= = 2\pi \sum_{n=-\infty}^\infty f_n \delta\left(\omega-{2\pi n\over L}\right)$

where $f_n$ are called Fourier coefficients:

$f_n = {1\over L} \int_0^L f(x) e^{-i2\pi n x / L} \d x$

We can see that the Fourier transform is zero for $\omega \neq {2\pi n\over L}$ . For $\omega={2\pi n\over L}$ it is equal to a delta function times a $2\pi$ multiple of a Fourier series coefficient. The delta functions structure is given by the period $L$ of the function $f(x)$ . All the information that is stored in the answer is inside the $f_n$ coefficients, so those are the only ones that we need to calculate and store.

The function $f(x)$ is calculated from the $f_n$ coefficients by applying the inverse Fourier transform to the final result of (3.4.6.4) as follows:

(3.4.6.5)¶ $f(x) = F^{-1} [F[f(x)](\omega)](x) = = F^{-1}\left[ 2\pi \sum_{n=-\infty}^\infty f_n \delta\left(\omega-{2\pi n\over L}\right) \right](x) = = {1\over 2\pi} \int_{-\infty}^\infty 2\pi \sum_{n=-\infty}^\infty f_n \delta\left(\omega-{2\pi n\over L}\right) e^{i\omega x} \d \omega = = \sum_{n=-\infty}^\infty f_n e^{i 2\pi n x / L}$

The expansion (3.4.6.5) is called a Fourier series. It is given by the Fourier coefficients $f_n$ . The equation (3.4.6.4) provides the relation between a Fourier transform and a Fourier series.

For example for $f(x) = \sin(x)$ , the only nonzero Fourier coefficients for $L=2\pi$ are $f_{-1} = {i\over2}$ and $f_1 = -{i\over2}$ . The Fourier transform then is:

$F[\sin(x)](\omega) = 2\pi\left(f_{-1}\delta(\omega-(-1)) + f_1\delta(\omega-1)\right) = = 2\pi\left({i\over 2}\delta(\omega+1)) -{i\over2}\delta(\omega-1)\right) = i\pi\delta(\omega+1) - i\pi\delta(\omega-1)$

For $f(x) = 1$ the only nonzero Fourier coefficient is $f_0=1$ , the Fourier transform then is:

$F[1](\omega) = 2\pi f_0 \delta(\omega-0) = 2\pi \delta(\omega)$

For $f(x) = e^{3ix}$ the only nonzero Fourier coefficient for $L=2\pi$ is $f_3=1$ , the Fourier transform then is:

$F[e^{3ix}](\omega) = 2\pi f_3 \delta(\omega-3) = 2\pi \delta(\omega-3)$

For $f(x) = \sum_{n=-\infty}^\infty \delta(x-2\pi n)$ the Fourier coefficients for $L=2\pi$ are all equal to $f_n={1\over 2\pi}$ and the Fourier transform is:

$F[f(x)](\omega) = 2\pi \sum_{n=-\infty}^\infty f_n \delta(\omega-n) = \sum_{n=-\infty}^\infty \delta(\omega-n)$

Note: if we start from (3.4.6.5), for simplicity on an interval $[-\pi, \pi]$ :

(3.4.6.6)¶ $f(x) = \sum_{n=-\infty}^\infty f_n e^{inx}$

To calculate the Fourier coefficients $f_n$ , we can just multiply both sides of (3.4.6.6) by $e^{-imx}$ and integrate:

$\int_{-\pi}^\pi f(x) e^{-imx} \d x = \int_{-\pi}^\pi \sum_{n=-\infty}^\infty f_n e^{inx} e^{-imx} \d x = = \sum_{n=-\infty}^\infty f_n \int_{-\pi}^\pi e^{i(n-m)x} \d x = = \sum_{n=-\infty}^\infty f_n 2\pi\delta_{nm} = = 2\pi f_m \,,$

so

(3.4.6.7)¶ $f_n = {1\over2\pi} \int_{-\pi}^\pi f(x) e^{-inx} \d x$

Convergence of Fourier Series¶

To see what conditions the function $f(x)$ must satisfy in order for the Fourier series to converge towards it, we can do the following analysis. Substituting (3.4.6.7) into (3.4.6.6) yields:

$f(x) = \sum_{n=-\infty}^\infty {1\over 2\pi}\int_{-\pi}^\pi f(y) e^{-iny} \d y\, e^{inx} = = \lim_{N\to\infty} \int_{-\pi}^\pi {1\over 2\pi}\sum_{n=-N}^N e^{in(x-y)} f(y) \d y = = \lim_{N\to\infty} \int_{-\pi}^\pi D_N(x-y) f(y) \d y$

We can now calculate the difference between the Fourier series and the function value:

$\lim_{N\to\infty} \int_{-\pi}^\pi D_N(x-y) f(y) \d y - f(x) = = \lim_{N\to\infty} \int_{-\pi}^\pi D_N(x-y) \left(f(y)-f(x)\right) \d y = = \lim_{N\to\infty} \int_{-\pi}^\pi {\sin\left(\left(N+\half\right)(x-y)\right) \over 2\pi \sin\left(x-y\over2\right)} \left(f(y)-f(x)\right) \d y = = \lim_{N\to\infty} \int_{-\pi}^\pi {f(y)-f(x)\over 2\pi\sin\left(x-y\over 2\right)} \sin\left(\left(N+\half\right)(x-y)\right) \d y = = \lim_{N\to\infty} \int_{x-\pi}^{x+\pi} {f(x-u)-f(x)\over 2\pi\sin\left(u\over 2\right)} \sin\left(\left(N+\half\right)u\right) \d u = = \lim_{N\to\infty} \int_{x-\pi}^{x+\pi} h(u) \sin\left(\left(N+\half\right)u\right) \d u = 0$

where $h(u)$ is finite and well behaved at the origin $u=0$ :

$h(u) = {f(x-u)-f(x)\over 2\pi\sin\left(u\over 2\right)} = - {f'(x)\over \pi} + {f''(x)\over 2\pi} u + O(u^2)$

The integral is zero because the more and more oscillating $\sin$ function cancels the contributions of positive and negative parts of the integrand. This can be proven explicitly as follows using the fact that $h(x)$ , $h'(x)$ and $\cos(Nx)$ is bounded as $N\to\infty$ :

$\lim_{N\to\infty} \int_a^b h(x) \sin(Nx)\,\d x = = \lim_{N\to\infty}{1\over N} \left(\left[-h(x)\cos(Nx)\right]_a^b + \int_a^b h'(x) \cos(Nx)\,\d x \right) = 0$

The conditions that we used are that the function $h(u)$ can be integrated, which is satisfied if e.g. $f(x)$ has derivatives. These conditions can be loosened in various ways.

3.5. Fourier Transform of a Periodic Function (e.g. in a Crystal)¶

The Fourier transform in (3.4.1) requires the function $f(\mathbf{x})$ to be decaying fast enough in order to converge. In an infinite crystal, on the other hand, the function $f(\mathbf{x})$ is typically periodic (and thus not decaying):

$f(\mathbf{x}+\mathbf{T}(n_1, n_2, n_3)) = f(\mathbf{x})$

where $\mathbf{T}(\mathbf{n})=\mathbf{T}(n_1, n_2, n_3)=n_1\mathbf{a}_1+n_2\mathbf{a}_2+n_3\mathbf{a}_3$ are the crystal translation vectors. As such, the Fourier transform in (3.4.1) is infinite, but it can be made finite by the following definition:

(3.5.1)¶ $F[f(\mathbf{x})] \equiv \tilde f(\bomega) = {1\over\Omega_\mathrm{crystal}}\int_{\Omega_\mathrm{crystal}} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = {1\over\Omega_\mathrm{crystal}} \sum_\mathbf{n} \int_{\Omega_\mathrm{cell}} f(\mathbf{x}+\mathbf{T}(\mathbf{n})) e^{-i\bomega \cdot (\mathbf{x}+\mathbf{T}(\mathbf{n}))}\,\d^3 x = = {1\over\Omega_\mathrm{crystal}} \sum_\mathbf{n} \int_{\Omega_\mathrm{cell}} f(\mathbf{x}) e^{-i\bomega \cdot (\mathbf{x}+\mathbf{T}(\mathbf{n}))}\,\d^3 x = = {1\over\Omega_\mathrm{crystal}} \sum_\mathbf{n} e^{-i\bomega \cdot \mathbf{T}(\mathbf{n})} \int_{\Omega_\mathrm{cell}} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = {1\over\Omega_\mathrm{crystal}} N_\mathrm{cell} \int_{\Omega_\mathrm{cell}} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x = = {1\over\Omega_\mathrm{cell}} \int_{\Omega_\mathrm{cell}} f(\mathbf{x}) e^{-i\bomega \cdot \mathbf{x}}\,\d^3 x$

This assumes that the wave vector $\bomega=\mathbf{G}$ is equal to the reciprocal space vectors $\mathbf{G}$ , defined by

(3.5.2)¶ $e^{i\mathbf{G} \cdot \mathbf{T}(\mathbf{n})} = 1\,,$

because then $\sum_\mathbf{n} e^{-i\bomega \cdot \mathbf{T}(\mathbf{n})} = \sum_\mathbf{n} 1 = N_\mathrm{cell}$ .

For $\bomega\neq\mathbf{G}$ , the expression ${1\over\Omega_\mathrm{crystal}} \sum_\mathbf{n} e^{-i\bomega \cdot \mathbf{T}(\mathbf{n})} = 0$ vanishes, because the sum is bounded, and so dividing by the (infinite) crystal volume makes the expression vanish, and so $\tilde f(\bomega)=0$ . In other words, the only non-zero Fourier components $\tilde f(\bomega)$ of any periodic function $f(\mathbf{x})$ are those with $\bomega=\mathbf{G}$ . Equivalently said, if the Fourier components of a given function are non-zero for some $\bomega\neq\mathbf{G}$ , then the function is not periodic.

Summary: the only difference between the crystal Fourier transform (3.5.1) and the usual Fourier transform (3.4.1) is the $\Omega_\mathrm{crystal}$ factor. The Fourier transform (3.5.1) of a periodic function is nonzero only for $\omega=\mathbf{G}$ and is equal to:

(3.5.3)¶ $F[f(\mathbf{x})] \equiv \tilde f(\mathbf{G}) = {1\over\Omega_\mathrm{cell}} \int_{\Omega_\mathrm{cell}} f(\mathbf{x}) e^{-i\mathbf{G} \cdot \mathbf{x}}\,\d^3 x$

Note: the fact that the sum is bounded follows from:

$\left| \sum_{n=-\infty}^\infty e^{ikn} \right| = \left| \lim_{N\to\infty} \sum_{n=-N}^N e^{ikn} \right| = \left| \lim_{N\to\infty} \left(1+2\sum_{n=1}^N \cos kn\right) \right|= = \left| \lim_{N\to\infty} {\cos kN - \cos k(N+1) \over 1-\cos k} \right| < {2 \over | 1-\cos k | }$

Because $| \cos kN - \cos k(N+1) | < 2$ . So for $k\neq 2\pi$ (i.e. the denominator is non-zero), the sum is bounded (to be precise, the infinite sum does not converge, because it oscillates, but the point is that the partial sum is always bounded). For $k=2\pi$ , the sum is infinite, because $e^{i2\pi n} = 1$ .

Since we divided the direct Fourier transform in (3.4.1) by $\Omega_\mathrm{crystal}$ to obtain (3.5.1), we need to multiply the inverse transform in (3.4.1) by $\Omega_\mathrm{crystal}$ :

(3.5.4)¶ $F^{-1}[\tilde f(\bomega)] = f(\mathbf{x}) = {\Omega_\mathrm{crystal}\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega = = {\Omega_\mathrm{cell}N_\mathrm{cell}\over(2\pi)^3}\int_{-\infty}^{\infty} \tilde f(\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega = = {N_\mathrm{cell}\over\Omega_\mathrm{BZ}} \sum_{\mathbf{G}} \int_{\Omega_\mathrm{BZ}} \tilde f(\mathbf{G}+\bomega) e^{+i(\mathbf{G}+\bomega) \cdot \mathbf{x}}\,\d^3 \omega = = {N_\mathrm{cell}\over\Omega_\mathrm{BZ}} \sum_{\mathbf{G}} e^{+i\mathbf{G} \cdot \mathbf{x}} \int_{\Omega_\mathrm{BZ}} \tilde f(\mathbf{G}+\bomega) e^{+i\bomega \cdot \mathbf{x}}\,\d^3 \omega = = \sum_{\mathbf{G}} \tilde f(\mathbf{G}) e^{+i\mathbf{G} \cdot \mathbf{x}} \int_{\Omega_\mathrm{BZ}} \delta(\boldsymbol\omega) e^{+i\boldsymbol\omega \cdot \mathbf{x}}\,d^3 \omega = = \sum_{\mathbf{G}} \tilde f(\mathbf{G}) e^{+i\mathbf{G} \cdot \mathbf{x}}$

where we used the fact that:

${N_\mathrm{cell}\over\Omega_\mathrm{BZ}}\tilde f(\mathbf{G}+\boldsymbol\omega) =\tilde f(\mathbf{G})\delta(\boldsymbol\omega) \,.$

Alternatively, if one is only interested to show that the inverse transformation works, one can directly substitute the direct formula (3.5.3) into (3.5.4) as follows:

$F^{-1}[\tilde f(\mathbf{G})] = \sum_{\mathbf{G}} \tilde f(\mathbf{G}) e^{+i\mathbf{G} \cdot \mathbf{x}} = = \sum_{\mathbf{G}} \left({1\over\Omega_\mathrm{cell}} \int_{\Omega_\mathrm{cell}} f(\mathbf{x'}) e^{-i\mathbf{G} \cdot \mathbf{x'}}\,d^3 x'\right) e^{+i\mathbf{G} \cdot \mathbf{x}} = = {1\over\Omega_\mathrm{cell}} \int_{\Omega_\mathrm{cell}} f(\mathbf{x'}) \sum_{\mathbf{G}} e^{i\mathbf{G} \cdot (\mathbf{x}-\mathbf{x'})}\,d^3 x' = = {1\over\Omega_\mathrm{cell}} \int_{\Omega_\mathrm{cell}} f(\mathbf{x'}) (2\pi)^3 \delta\left({(2\pi)^3\over\Omega_\mathrm{cell}} (\mathbf{x}-\mathbf{x'})\right) \,d^3 x' = = {1\over\Omega_\mathrm{cell}} \int_{\Omega_\mathrm{cell}} f(\mathbf{x'}) (2\pi)^3 {\Omega_\mathrm{cell}\over (2\pi)^3} \delta(\mathbf{x}-\mathbf{x'}) \,d^3 x' = =f(\mathbf{x})\,,$

where we used the fact that:

$\sum_{n=-\infty}^\infty e^{inx} = 2\pi\delta(x)\,.$

Thus we have shown that $F^{-1}[\tilde f(\mathbf{G})] = f(\mathbf{x})$ .

3.5.1. One Dimension (Fourier Series)¶

In one dimension with a periodic function $f(x+L)=f(x)$ , the volume of a unit cell is $\Omega_\mathrm{cell}=L$ and the reciprocal space vectors $G$ are defined using $e^{iGL}=1$ from which $G_k = {2\pi\over L} k$ . The equation (3.5.3) then becomes:

(3.5.1.1)¶ $F[f(x)] \equiv \tilde f(G_k) \equiv c_k = {1\over L} \int_{-{L\over2}}^{L\over2} f(x) e^{-i G_k x}\,\d x = {1\over L} \int_{-{L\over2}}^{L\over2} f(x) e^{-i(2\pi k x/L)}\,\d x$

This is exactly the definition of a Fourier series ( $c_k$ are the Fourier coefficients). The inverse transform follows from (3.5.4):

(3.5.1.2)¶ $f(x) = \sum_{k=-\infty}^\infty \tilde f(G_k) e^{i G_k x} = \sum_{k=-\infty}^\infty c_k e^{i(2\pi kx/L)}$

3.6. Discrete Fourier Transform¶

In the discrete case, we only have a finite number $N$ of reciprocal points:

$k=0, 1, \dots, N/2-1, -N/2, -N/2+1, \dots, -1 \quad\mbox{if $N$ is even} k=0, 1, \dots, (N-1)/2, -(N-1)/2, -(N-1)/2+1, \dots, -1 \quad\mbox{if $N$ is odd}$

E.g. for:

$N=8 \quad \mbox{we get} \quad k=0, 1, 2, 3, -4, -3, -2, -1 N=9 \quad \mbox{we get} \quad k=0, 1, 2, 3, 4, -4, -3, -2, -1$

The real space function $f(x)$ is sampled at points $x_n={L\over N}n$ for $n=-N/2,\dots,N/2-1$ and the equation (3.5.1.1) becomes:

$c_k = {1\over L} \int_{-{L\over2}}^{L\over2} f(x) e^{-i(2\pi k x/L)}\,\d x = = \lim_{N\to\infty} {1\over L}\sum_{n=-N/2}^{N/2-1} f(x_n) e^{-i(2\pi k x_n/L)}\,{L\over N} = = \lim_{N\to\infty} {1\over N}\sum_{n=-N/2}^{N/2-1} f(x_n) e^{-2\pi i {k\over N} n}$

The equation (3.5.1.2) becomes:

$f(x_n) = \sum_{k=-\infty}^\infty c_k e^{i(2\pi kx_n/L)} = = \lim_{N\to\infty} \sum_{k=-N/2}^{N/2-1} c_k e^{i(2\pi kx_n/L)} = = \lim_{N\to\infty} \sum_{k=-N/2}^{N/2-1} c_k e^{2\pi i {k\over N} n}$

Using the fact

$x_n + L = {L\over N}n + L = {L\over N}(n + N) = x_{n+N}\,,$

we can express the periodicity $f(x_n+L)=f(x_n)$ as $f(x_{n+N})=f(x_n)$ . The sums can then be rearranged:

$c_k = \lim_{N\to\infty} {1\over N}\sum_{n=-N/2}^{N/2-1} f(x_n) e^{-2\pi i {k\over N} n} = = \lim_{N\to\infty} {1\over N} \left( \sum_{n=-N/2}^{-1} f(x_n) e^{-2\pi i {k\over N} n} + \sum_{n=0}^{N/2-1} f(x_n) e^{-2\pi i {k\over N} n} \right) = = \lim_{N\to\infty} {1\over N} \left( \sum_{n=N/2}^{N-1} f(x_{n-N}) e^{-2\pi i {k\over N} (n-N)} + \sum_{n=0}^{N/2-1} f(x_n) e^{-2\pi i {k\over N} n} \right) = = \lim_{N\to\infty} {1\over N} \sum_{n=0}^{N-1} f(x_n) e^{-2\pi i {k\over N} n}$

and if we drop the limit and consider a finite $N$ only:

$f(x_n) = \sum_{k=-N/2}^{N/2-1} c_k e^{2\pi i {k\over N} n} = = \left( \sum_{k=-N/2}^{-1} c_k e^{2\pi i {k\over N} n} + \sum_{k=0}^{N/2-1} c_k e^{2\pi i {k\over N} n} \right) = = \left( \sum_{k=N/2}^{N-1} c_{k-N} e^{2\pi i {(k-N)\over N} n} + \sum_{k=0}^{N/2-1} c_k e^{2\pi i {k\over N} n} \right) = = \sum_{k=0}^{N-1} c_k e^{2\pi i {k\over N} n}$

Summary, the direct transform:

(3.6.1)¶ $c_k = {1\over N} \sum_{n=0}^{N-1} f(x_n) e^{-2\pi i {k\over N} n}$

and inverse transform:

(3.6.2)¶ $f(x_n) = \sum_{k=0}^{N-1} c_k e^{2\pi i {k\over N} n}\,,$

with $x_n={L\over N}n$ . In the limit $N\to\infty$ , the equation (3.6.1) becomes (3.5.1.1) and equation (3.6.2) becomes (3.5.1.2) and as we increase $N$ , the discrete Fourier transform numerically converges towards the Fourier series results.

The ${1\over N}$ factor is sometimes moved from the direct to the inverse transform, but then the correspondence with Fourier series is broken (one has to divide and multiply by $N$ appropriately to recover it).

3.7. Fast Fourier Transform (FFT)¶

We write the discrete Fourier transform (3.6.1) using a notation more commonly used for FFTs:

$X(k) = \sum_{n=0}^{N-1} x(n) W_N^{kn}$

where:

$W_N = e^{-2\pi i / N}$

Similarly, the inverse discrete Fourier transform (3.6.2) becomes:

$x(n) = {1\over N} \sum_{k=0}^{N-1} X(k) W_N^{-kn}$

3.7.1. Decimation In Frequency (DIF)¶

We start with radix-4:

$X(k) = \sum_{n=0}^{N-1} x(n) W_N^{kn} = =\sum_{n=0}^{{N\over4}-1} x(n) W_N^{kn} +\sum_{n={N\over4}}^{{2N\over4}-1} x(n) W_N^{kn} +\sum_{n={2N\over4}}^{{3N\over4}-1} x(n) W_N^{kn} +\sum_{n={3N\over4}}^{{4N\over4}-1} x(n) W_N^{kn} = =\sum_{n=0}^{{N\over4}-1} \left[ x(n) W_N^{kn} +x\left(n+{N\over4}\right) W_N^{k\left(n+{N\over4}\right)} +x\left(n+{2N\over4}\right) W_N^{k\left(n+{2N\over4}\right)} +x\left(n+{3N\over4}\right) W_N^{k\left(n+{3N\over4}\right)} \right] = =\sum_{n=0}^{{N\over4}-1} \left[ x(n) +x\left(n+{N\over4}\right) W_N^{kN\over4} +x\left(n+{2N\over4}\right) W_N^{2kN\over4} +x\left(n+{3N\over4}\right) W_N^{3kN\over4} \right] W_N^{kn} = =\sum_{n=0}^{{N\over4}-1} \left[ x(n) +x\left(n+{N\over4}\right) (-i)^k +x\left(n+{2N\over4}\right) (-1)^k +x\left(n+{3N\over4}\right) i^k \right] W_N^{kn}$

Now we subdivide the $X(k)$ sequence into 4 subsequences:

$X(4k) = \sum_{n=0}^{{N\over4}-1} \left[ x(n) +x\left(n+{N\over4}\right) (-i)^{4k} +x\left(n+{2N\over4}\right) (-1)^{4k} +x\left(n+{3N\over4}\right) i^{4k} \right] W_N^{4kn} = = \sum_{n=0}^{{N\over4}-1} \left[ x(n) +x\left(n+{N\over4}\right) +x\left(n+{2N\over4}\right) +x\left(n+{3N\over4}\right) \right] W_{N\over4}^{kn}$

Similarly:

$X(4k+1) = \sum_{n=0}^{{N\over4}-1} \left[ x(n) -i x\left(n+{N\over4}\right) -x\left(n+{2N\over4}\right) +i x\left(n+{3N\over4}\right) \right] W_N^{n} W_{N\over4}^{kn} X(4k+2) = \sum_{n=0}^{{N\over4}-1} \left[ x(n) -x\left(n+{N\over4}\right) +x\left(n+{2N\over4}\right) -x\left(n+{3N\over4}\right) \right] W_N^{2n} W_{N\over4}^{kn} X(4k+3) = \sum_{n=0}^{{N\over4}-1} \left[ x(n) +i x\left(n+{N\over4}\right) -x\left(n+{2N\over4}\right) -i x\left(n+{3N\over4}\right) \right] W_N^{3n} W_{N\over4}^{kn}$

This has a form of a DFT of length ${N\over4}$ :

$X(4k) = \sum_{n=0}^{{N\over4}-1} F_0(n) W_{N\over4}^{kn} X(4k+1) = \sum_{n=0}^{{N\over4}-1} F_1(n) W_{N\over4}^{kn} X(4k+2) = \sum_{n=0}^{{N\over4}-1} F_2(n) W_{N\over4}^{kn} X(4k+3) = \sum_{n=0}^{{N\over4}-1} F_3(n) W_{N\over4}^{kn}$

where

$\begin{pmatrix} F_0(n) \\ F_1(n) \\ F_2(n) \\ F_3(n) \\ \end{pmatrix} = \begin{pmatrix} x(n) +x\left(n+{N\over4}\right) +x\left(n+{2N\over4}\right) +x\left(n+{3N\over4}\right) \\ x(n) -i x\left(n+{N\over4}\right) -x\left(n+{2N\over4}\right) +i x\left(n+{3N\over4}\right) \\ x(n) -x\left(n+{N\over4}\right) +x\left(n+{2N\over4}\right) -x\left(n+{3N\over4}\right) \\ x(n) +i x\left(n+{N\over4}\right) -x\left(n+{2N\over4}\right) -i x\left(n+{3N\over4}\right) \end{pmatrix} = =\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix} \begin{pmatrix} x(n) \\ x\left(n+{N\over4}\right) \\ x\left(n+{2N\over4}\right) \\ x\left(n+{3N\over4}\right) \end{pmatrix}$

This coefficient matrix for various radix-n schemes can be generated by:

>>> from sympy import exp, I, pi, pprint, Matrix
>>> n = 2
>>> Matrix(n, n, lambda i, j: exp(-2*pi*I*i*j/n))
[1  1]
[1 -1]
>>> n = 3
>>> Matrix(n, n, lambda i, j: exp(-2*pi*I*(i*j % n)/n))
[1,              1,              1]
[1, exp(-2*I*pi/3), exp(-4*I*pi/3)]
[1, exp(-4*I*pi/3), exp(-2*I*pi/3)]
>>> n = 4
>>> Matrix(n, n, lambda i, j: exp(-2*pi*I*i*j/n))
[1  1  1  1]
[1 -I -1  I]
[1 -1  1 -1]
[1  I -1 -I]
>>> n = 5
>>> Matrix(n, n, lambda i, j: exp(-2*pi*I*(i*j % n)/n))
[1,              1,              1,              1,              1]
[1, exp(-2*I*pi/5), exp(-4*I*pi/5), exp(-6*I*pi/5), exp(-8*I*pi/5)]
[1, exp(-4*I*pi/5), exp(-8*I*pi/5), exp(-2*I*pi/5), exp(-6*I*pi/5)]
[1, exp(-6*I*pi/5), exp(-2*I*pi/5), exp(-8*I*pi/5), exp(-4*I*pi/5)]
[1, exp(-8*I*pi/5), exp(-6*I*pi/5), exp(-4*I*pi/5), exp(-2*I*pi/5)]
>>> n = 8
>>> Matrix(n, n, lambda i, j: exp(-2*pi*I*(i*j % n)/n))
[1,              1,  1,              1,  1,              1,  1,              1]
[1,   exp(-I*pi/4), -I, exp(-3*I*pi/4), -1, exp(-5*I*pi/4),  I, exp(-7*I*pi/4)]
[1,             -I, -1,              I,  1,             -I, -1,              I]
[1, exp(-3*I*pi/4),  I,   exp(-I*pi/4), -1, exp(-7*I*pi/4), -I, exp(-5*I*pi/4)]
[1,             -1,  1,             -1,  1,             -1,  1,             -1]
[1, exp(-5*I*pi/4), -I, exp(-7*I*pi/4), -1,   exp(-I*pi/4),  I, exp(-3*I*pi/4)]
[1,              I, -1,             -I,  1,              I, -1,             -I]
[1, exp(-7*I*pi/4),  I, exp(-5*I*pi/4), -1, exp(-3*I*pi/4), -I,   exp(-I*pi/4)]

One then recursively solves the smaller problems. This approach is used for example in FFTPACK. There are also other approaches how to decompose the DFT, used in various other libraries.

3.8. Laplace Transform¶

Laplace transform of $f(x)$ is:

$L[f(x)] = \int_0^{\infty} f(x) e^{-s x}\,\d x L^{-1}[\bar f(s)] = {1\over2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty} \bar f(s) e^{s x}\,\d s = \sum_{s_0} \res_{s=s_0} (\bar f(s) e^{s x})$

The contour integration is over the vertical line $\sigma+i\omega$ and $\sigma$ is chosen large enough so that all residues are to the left of the line (that’s because the Laplace transform $\bar f(s)$ is only defined for $s$ larger than the residues, so we have to integrate in this range as well). It can be shown that the integral over the left semicircle goes to zero:

$\left|\int_\Omega e^{sx}g(s) \d s \right| =\left|\int_{\pi\over2}^{3\pi\over2} e^{(\sigma + Re^{i\varphi})x} g(\sigma+Re^{i\varphi})iRe^{i\varphi}\d\varphi\right| \le \le R \max_\Omega |g(z)| e^{\sigma x} \int_{\pi\over2}^{3\pi\over2}\left| e^{xRe^{i\varphi}} \right|\d\varphi = = R \max_\Omega |g(z)| e^{\sigma x} \int_{\pi\over2}^{3\pi\over2}e^{xR \cos \varphi} \d\varphi = = R \max_\Omega |g(z)| e^{\sigma x} \int_0^{\pi}e^{-xR \sin \varphi} \d\varphi = < {\pi e^{\sigma x}\over x} \max_\Omega |g(z)|$

so the complex integral is equal to the sum of all residues of $\bar f(s)e^{sx}$ in the complex plane.

To show that it works:

$L^{-1} L [f(x)] = {1\over2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty} \left[\int_0^{\infty} f(x) e^{-s x}\,\d x\right] e^{s x}\,\d s = {1\over2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty} \left[\int_0^{\infty} f(x') e^{-s x'}\,\d x'\right] e^{s x}\,\d s = = \int_0^{\infty} f(x') \left[{1\over2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} e^{s (x- x')}\,\d s \right] \,\d x' = \int_0^{\infty} f(x') \delta(x-x') \,\d x' =f(x)$

where we used:

${1\over2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} e^{s (x- x')}\,\d s = {1\over2\pi i} \int_{\sigma-i\infty}^{\sigma+i\infty} e^{s (x- x')}\,\d s = {1\over2\pi i} \int_{-\infty}^{\infty} e^{(\sigma+i\omega) (x- x')}\,i\d \omega = = {e^{\sigma (x- x')}\over2\pi} \int_{-\infty}^{\infty} e^{i\omega (x- x')}\,\d \omega = e^{\sigma (x- x')}\delta(x - x') =\delta(x - x')$

and it can be derived from the Fourier transform by transforming a function $U(x)$ :

$U(x) = \begin{cases} f(x)e^{-\sigma x} &\text{for $x\ge0$}\cr 0 &\text{for $x<0$}\cr \end{cases}$

and making a substitution $s = \sigma + i\omega$ :

$L[f(x)] \equiv \bar f(s) = F[U(x)] \equiv \tilde U(\omega) = \int_{-\infty}^{\infty} U(x) e^{-i\omega x}\,\d x = \int_0^{\infty} f(x) e^{-\sigma x} e^{-i\omega x}\,\d x = \int_0^{\infty} f(x) e^{-s x}\,\d x L^{-1}[\bar f(s)] \equiv f(x) = U(x) e^{\sigma x} = F^{-1}[\tilde U(\omega)]e^{\sigma x} = F^{-1}[\bar f(s)]e^{\sigma x} = F^{-1}[\bar f(\sigma+i\omega)e^{\sigma x}] = {1\over2\pi}\int_{-\infty}^{\infty} \bar f(\sigma + i\omega)e^{\sigma x} e^{i\omega x}\,\d \omega = {1\over2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty} \bar f(s) e^{s x}\,\d s = \sum_{s_0} \res_{s=s_0} (\bar f(s) e^{s x})$

Where the bar ( $\bar f$ ) means the Laplace transform and tilde ( $\tilde U$ ) means the Fourier transform.

3.9. Hilbert Transform¶

The Hilbert transform is:

$H[f(x)](t) \equiv \bar f(t) = \mathrm{p.v.} {1\over\pi t} * f(t) = {1\over\pi} \mathrm{p.v.} \int_{-\infty}^{\infty} {f(x)\over t-x}\,\d x$

By applying the Fourier transform to both sides of the equation, we get:

$F[H[f(x)](t)](\omega) = F\left[\mathrm{p.v.} {1\over\pi t} * f(t)\right](\omega) = F\left[\mathrm{p.v.} {1\over\pi t}\right](\omega) F[f(t)](\omega) = -i\sign(\omega) F[f(t)](\omega)$

So the Hilbert transform can be calculated using a Fourier transform as:

$H[f](t) = F^{-1}\left[-i\sign(\omega) F[f](\omega)\right](t)$

The inverse Hilbert transform can then be calculated by inverting:

$F^{-1}\left[-i\sign(\omega) F[f(x)](\omega)\right](t) = H[f](t) = \bar f(t) -i\sign(\omega) F[f(x)](\omega) = F[\bar f(t)](\omega) F[f(x)](\omega) = i \sign(\omega) F[\bar f(t)](\omega) f(x) = F^{-1}\left[i \sign(\omega) F[\bar f(t)](\omega)\right](x)$

so we get:

$H^{-1}[\bar f(t)](x) = f(x) = F^{-1}\left[i \sign(\omega) F[\bar f(t)](\omega)\right](x)= = - F^{-1}\left[-i \sign(\omega) F[\bar f(t)](\omega)\right](x) = = - H[\bar f(t)](x) = = - {1\over\pi} \mathrm{p.v.} \int_{-\infty}^{\infty} {\bar f(t)\over x-t}\,\d t$

From this it also follows:

$\bar f(t) = - H[H[\bar f(t)](x)](t)$

or

$H[H[f(x)](t)](x) = -f(x)$

In other words, by applying the Hilbert transform twice, the result is the negative of a function.