Discrete Fourier Transform

Consider the following integral:

$$F(k) = \int_0^{2\pi}f(x) e^{-ikx}\,\textrm{d}x$$ (2.1)

Assume that function f(x) is sampled within the interval $[0, 2\pi]$on a regularly spaced set of points $x_l,\> l=0\ldots n-1$. In that case we can approximate integral (2.1) by evaluating a finite sum

 

$$F_k = \alpha \sum_{l=0}^{n-1} f(x_l) e^{-ik(2\pi/n)l}$$

$$= \alpha\bigg( f(x_0) e^{-ik(2\pi/n)0} + f(x_1) e^{-ik(2\pi/n)1} + f(x_2) e^{-ik(2\pi/n)2} + \cdots$$

$$\quad + f(x_{n-1}) e^{-ik(2\pi/n)(n - 1)} \bigg)$$ (2.2)

If every point x_l sits at the left hand side of its interval of length $\Delta x = 2\pi/n$, then the above evaluates to:

$$F_k = \alpha\bigg( f(x_0) e^{-ik x_0} + f(x_1) e^{-ik x_1} + f(x_2) e^{-ik x_2} + \cdots$$

$$\quad + f(x_{n-1}) e^{-ik x_{n - 1}} \bigg)$$ (2.3)

and for $\alpha = \Delta x = 2\pi/n$ we get that

$$F_k \approx F(k)$$ (2.4)

Let us substitute $\omega_n = e^{-2\pi i/n}$ in equation (2.2):

$$F_k = \alpha \sum_{l=0}^{n-1} \omega_n^{kl}f(x_l)$$ (2.5)

We can rewrite this expression in the form of a matrix multiplication:

$$\left( \begin{array}{l} F_0 \\ F_1 \\ F_2 \\ \vdots ... ...1) \\ f(x_2) \\ \vdots \\ f(x_{n-1}) \end{array}\right)$$ (2.6)

The discrete Fourier transform matrix shown in equation (2.6) is a complex Vandermonde matrix. Such matrices appear in polynomial interpolation. It is easy to write a simple Fortran program that builds a discrete Fourier transform matrix. But instead of doing it with Fortran, this time we'll use our interactive Fortran, Octave, to build a $4\times4$ discrete Fourier transform matrix.

Octave, version 2.0.13 (sparc-sun-solaris2.6).
Copyright (C) 1996, 1997, 1998 John W. Eaton.
This is free software with ABSOLUTELY NO WARRANTY.
For details, type `warranty'.

octave:1> w = exp(-2 * pi * i / 4)
w = -i
octave:2> F = ones(4,4)
F =

  1  1  1  1
  1  1  1  1
  1  1  1  1
  1  1  1  1

octave:3> (0:3)  
ans =

  0  1  2  3

octave:4> F(2,:) = w .^ (0:3)
F =

  1  1  1  1
  1 -i -1  i
  1  1  1  1
  1  1  1  1

octave:5> F(3,:) = w .^ (2 * (0:3))
F =

  1  1  1  1
  1 -i -1  i
  1 -1  1 -1
  1  1  1  1

octave:6> F(4,:) = w .^ (3 * (0:3))
F =

  1  1  1  1
  1 -i -1  i
  1 -1  1 -1
  1  i -1 -i

octave:7>

Another way to accomplish the same would be as follows:

octave:8> F = ones(4,4)
F =

  1  1  1  1
  1  1  1  1
  1  1  1  1
  1  1  1  1

octave:9> F(2,:) = w .^ (0:3)
F =

  1  1  1  1
  1 -i -1  i
  1  1  1  1
  1  1  1  1 

octave:10> F(3,:) = F(2,:) .* F(2,:)
F =

  1  1  1  1
  1 -i -1  i
  1 -1  1 -1
  1  1  1  1

octave:11> F(4,:) = F(2,:) .* F(3,:)
F =

  1  1  1  1
  1 -i -1  i
  1 -1  1 -1
  1  i -1 -i

octave:13>

This is a faster way because multiplication is cheaper than exponentiation.

Once we have built the discrete Fourier transform matrix we can attempt to find a numerical Fourier transform of, say, a sine function:

octave:39> x = [ 0, 2*pi/4, 2*2*pi/4, 3*2*pi/4 ]
x =

  0.00000  1.57080  3.14159  4.71239

octave:40> x = sin(x)
x =

   0.00000   1.00000   0.00000  -1.00000

octave:41> 2 * pi / 4 * F * x'
ans =

   0.00000 + 0.00000i
   0.00000 - 3.14159i
   0.00000 + 0.00000i
  -0.00000 + 3.14159i

octave:42>

Our numerical result, neglecting the last harmonic, is that

$$\int_0^{2\pi} \sin{x} \, e^{-ikx} \,\textrm{d}x = -\pi i \delta(k - 1)$$

which is just right, because

$$\begin{eqnarray*}&& \int_0^{2\pi} \sin{x} \, e^{-ikx} \,\textrm{d}x = \int_0^{2\... ...nt_0^{2\pi} \sin{x} \sin{x} \,\textrm{d}x = -i \delta(k - 1) \pi \end{eqnarray*}$$

Recall that you can check that integration easily with Emacs calc, e.g.,

alg' integ(sin(x)^2, x, 0, 2 * pi)
     pi

but remember to do it in the radians mode.

What about that last harmonic? Repeating the same computation for n = 8 returns:

octave:43> w = exp(-2 * pi * i / 8)
w = 0.70711 - 0.70711i
octave:44> F = ones(8,8);
octave:45> for k = 2:8
> F(k,:) = w .^ ((k - 1) * (0:7));
> end
octave:46> x = 2 * pi / 8 * (0:7);
octave:47> x = sin(x);
octave:48> 2 * pi / 8 * F * x'
ans =

  -0.00000 + 0.00000i
   0.00000 - 3.14159i
   0.00000 + 0.00000i
   0.00000 + 0.00000i
   0.00000 + 0.00000i
   0.00000 + 0.00000i
   0.00000 + 0.00000i
  -0.00000 + 3.14159i

octave:49>

So we are still getting the same result, $-i \pi \delta(k - 1)$, and the only difference is that the last harmonic has been moved further away.

Let us try a signal with higher harmonics:

octave:45> x = 2 * pi / 8 * (0:7);
octave:46> y = sin(x) + 0.5 * sin (2 * x);
octave:47> 2 * pi / 8 * F * y'
ans =

   0.00000 + 0.00000i
   0.00000 - 3.14159i
   0.00000 - 1.57080i
  -0.00000 + 0.00000i
   0.00000 + 0.00000i
   0.00000 + 0.00000i
  -0.00000 + 1.57080i
  -0.00000 + 3.14159i

octave:48> y = cos(x) + sin(2 * x) - cos(3 * x);
octave:49> 2 * pi / 8 * F * y'
ans =

   0.00000 + 0.00000i
   3.14159 + 0.00000i
   0.00000 - 3.14159i
  -3.14159 - 0.00000i
  -0.00000 + 0.00000i
  -3.14159 - 0.00000i
  -0.00000 + 3.14159i
   3.14159 + 0.00000i

octave:50>

What is happening is that harmonics are inserted from both ends. The reason for this is that the discrete Fourier transform returns amplitudes that correspond to both positive and negative wave numbers. Amplitudes for positive wave numbers are returned in slots labeled with $l = 1 \ldots n/2 - 1$. Amplitudes for negative wave numbers are returned in reverse order in slots $l = n/2 + 1 \dots n - 1$ and the Nyquist freqency amplitude

$$f_c = \frac{1}{2\Delta}$$ (2.7)

or

$$k_c = 2\pi f_c = \frac{2\pi}{2\Delta},$$ (2.8)

where $\Delta$ is the sampling interval, is returned in slot l = n/2. The amplitude that corresponds to the zero frequency, i.e., a constant shift, is returned in slot l=0. Because in Fortran arrays are normally numbered by default from 1 through n (but this can be changed with an appropriate declaration) you've got to add 1 to l in order to refer to an appropriate slot in a Fortran array.

In our case $\Delta = 2\pi/8$, k_c = 4, and the Nyquist frequency Amplitude is -0.0 + 0.0 i.

How did those negative frequencies creep into our sum:

$$F_k = \alpha \sum_{l=0}^{n-1} f(x_l) e^{-ik(2\pi/n)l}$$

Observe that F_k is periodic in k with period equal to n. Therefore:

F_{n - k} = F_-k (2.9)

In other words only half of the vector returned by a discrete Fourier transform corresponds to positive wave numbers. So, if you want, say, m amplitudes associated with positive wave numbers, you must sample your signal on at least 2m + 2 points^2.1. Conversely, if you sample on 2m + 2 points the maximum frequency amplitude will correpond not to $\omega = 2m + 1$ but to k_c = m+1.

What is going to happen if your signal contains harmonics that are higher than the Nyquist frequency k_c? It turns out that all of the power spectral density that lies outside of the range [-k_c, +k_c] will be sneakily moved into that range, i.e., any amplitude that corresponds to a frequency outside that range is going to be aliased into that range.

There is no way to overcome aliasing other than to know the bandwidth limit of a signal in advance and then sample sufficiently densely, so that the Nyquist frequency is higher than the bandwidth limit.