The Fast Fourier Transform

Introduction

Fourier analysis has many applications in signal processing; in this note we will discuss an algorithm that helps us more quickly determine understand a signal. Our motivating example comes from music theory, and breaking a sound wave into components.

Let f be a continuous 1-periodic function on $\mathbb{R}$ , i.e. assume f(x) = f(x+1) for all x. We know a family of continuous functions f with this property, namely

$\sin(2\pi n x), \cos(2\pi n x)$

and it turns out that these are the building block for all 1-periodic functions. The idea is that we can decompose f as

where for all real t,

We won't discuss issues of convergence here, but the idea is that any function of the form a sin(mt) + b cos(mt) is just a simple wave with period $2\pi/m.$ and if we are willing to use a large enough N, then

However, we could be dealing with a whole range of frequencies, not just integer frequencies, and in this case our combination of simple sound waves is more "smeared out", but f(x) is still a combination of the simple periodic functions $e^{2\pi i x \xi}$ for $\xi \in \mathbb{R}:$

where $\hat{f}$ is the Fourier transform of f.)

In the case of a discrete spectrum of frequencies, the coefficients be computed as

In the case of a continuous spectrum, we have

There are many physical applications of this, but we'll stick to one example. Assuming a measured sound wave is approximately a periodic function we can decompose it into simple waves with period n (frequency 1/n) and amplitude

Perhaps the most familiar idealized example of this is breaking a musical chord into individual notes, each of which has its own frequency.

Discrete Fourier Transform

In computer science we want to find a discrete approximation, and we do this by sampling a finite number of values of a function and approximating sums with integrals.

If we want to assume we are sampling over a long period of time and also sampling frequently, we can the frequencies we care about are from -M to M and we sample every value of x that is a multiple of 1/M and get a good enough approximation of our data on the whole real line.

(This gives us 2M^2 + 1 data points).

However, due to the nature of the exponential function, our analysis is going to be exactly the same for any collection of equally spaced points, and we can play around with the number of data points a little bit to make it more convenient.

For simplicity we'll say the value of f is sampled at N equally spaced points, and define the discrete Fourier transform as follows:

Then we can in fact say that for each of the values of n from 0 through N - 1,

Again, recall that f is the data that is detected by our microphone and the values of a_n are what we actually care about to determine the actual signal.

Now on the surface, it would appear that computing the whole list of values

requires $O(N^2)$ arithmetic operations in addition to calculating N values of the (complex) exponential function - after all, for each of N Fourier coefficients a_n we compute N products.

The Fast Fourier Transform

We can greatly reduce the number of computations needed by thinking of the Fourier transform as a whole instead of as a N individual values to compute, and in fact reduce to O(N log N) arithmetic operations for suitable N. We will describe this in the case where N is a power of 2. This method is known as the radix-2 Cooley-Tukey algorithm.

Let's define our FFT to be a function that takes as input an integer N = 2^m and a list L of values of length N. Here's python-like pseudo-code that computes the FFT much faster. (There are packages that handle complex numbers in python so this is actually pretty close to what we would write.)

def FFT(L, N):
    if N == 1:
        return L[0]
    evens = L[::2]
    odds = L[1::2]
    A = FFT(evens, N // 2)
    B = FFT(odds, N // 2)
    u = e^(-2pi i/N)
    U = [1]
    for n in range(N // 2):
        if n > 0:    
            U.append(U[-1]*u)
    C = [A[n] + U[n] * B[n] for n in range(N//2)]
    D = [A[n] - U[n] * B[n] for n in range(N//2)]
    return C + D

We remark that A + B is list concatenation, not vector addition. For the purposes of a human doing the verification, U[n] is just u^n, but we want have to compute all of these powers in O(n) time.

It really only takes a couple of lines of algebra to verify that the recursive formula for the Fourier transform is correct at each position: That is, we take the Fourier transforms for our function sampled at even-numbered positions only and at odd-numbered positions only, and simply check that for each k, our recursive formula has the correct coefficient of L[k]. Then we add them together, making one phase shift for the odd-numbered values to get the front half of the full sequence of coefficients, do the same with a different phase shift for the back half.

We can see there is a serious speedup. If an algorithm required quadratic time, two times through the algorithm on data of size N would still take $\Theta(N^2)$ additional steps to reach the amount of time needed to perform the algorithm on data of size N. But we can clearly see that only O(N) additional steps are needed to make this leap, so the FFT is definitely much faster than computation directly from the definition.