FFT algorithm and program

Fast Fourier transform: idea

The DFT of $N$ samples is $X_k=\sum_{n=0}^{N-1}x_n W_N^{kn}$ with $W_N=e^{-i2\pi/N}$, needing $N^2$ operations. The Cooley-Tukey radix-2 FFT (for $N=2^m$) splits the series into even and odd samples:

$$X_k=E_k+W_N^kO_k,\qquad X_{k+N/2}=E_k-W_N^kO_k,\quad k=0,\dots,\tfrac N2-1$$

and repeats the split $\log_2N$ times. The cost is $\tfrac N2\log_2N$ complex multiplications.

Algorithm

1. Read $N$ (a power of 2) and the samples $x_n$ at interval $\Delta t$.
2. Reorder the samples in bit-reversed order.
3. For stage length $len=2,4,8,\dots,N$: compute $W=e^{-i2\pi/len}$; for each block, for each $k=0\dots len/2-1$ form $u=x_{i+k}$, $v=x_{i+k+len/2}\,W^k$, then set $x_{i+k}=u+v$, $x_{i+k+len/2}=u-v$.
4. After the last stage $x$ holds $X_k$. Amplitude $=2|X_k|/N$, frequency $f_k=k/(N\Delta t)$.

 x0 --o----o----o-- X0
 x4 --o----|    |
 x2 --o----o    |   (butterflies:
 x6 --o----|    |    a+Wb , a-Wb)
 ...

Program (C)

#include <stdio.h>
#include <math.h>
#include <complex.h>
#define PI 3.14159265358979323846

/* in-place radix-2 FFT, n must be a power of 2 */
void fft(double complex x[], int n)
{
    for (int i = 1, j = 0; i < n; i++) {      /* bit reversal */
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) { double complex t = x[i]; x[i] = x[j]; x[j] = t; }
    }
    for (int len = 2; len <= n; len <<= 1) {  /* butterfly stages */
        double complex wlen = cexp(-2.0 * I * PI / len);
        for (int i = 0; i < n; i += len) {
            double complex w = 1.0;
            for (int k = 0; k < len / 2; k++) {
                double complex u = x[i + k];
                double complex v = x[i + k + len / 2] * w;
                x[i + k]           = u + v;
                x[i + k + len / 2] = u - v;
                w *= wlen;
            }
        }
    }
}

int main(void)
{
    int n = 8; double fs = 8.0; double complex x[8];
    for (int i = 0; i < n; i++) {             /* 1 Hz (amp 2) + 3 Hz (amp 1) */
        double t = i / fs;
        x[i] = 2*sin(2*PI*1*t) + 1*sin(2*PI*3*t);
    }
    fft(x, n);
    for (int k = 0; k < n; k++)
        printf("f=%.1f Hz  amp=%.3f\n", k*fs/n, 2*cabs(x[k])/n);
    return 0;
}

Example: time domain and frequency domain

Signal $x(t)=2\sin(2\pi t)+\sin(6\pi t)$ sampled at 8 Hz for 1 s ($N=8$). The program output is:

(components at $k=5,6,7$ are the mirror image of $k=3,2,1$.)

• Time domain shows the amplitude of vibration against time, peak value, duration, arrival time and decay of the signal, but the individual components are mixed together.
• Frequency domain shows the dominant frequencies (1 Hz and 3 Hz), their amplitudes (2 and 1) and phase. For a structure, the peak frequency is the natural frequency, so resonance with earthquake or machine frequency can be checked; noise is separated by filtering.
• Frequency resolution is $\Delta f=1/(N\Delta t)=1$ Hz; the highest frequency identified is the Nyquist frequency $f_s/2=4$ Hz.