4.9b - FFT: history and algorithms

See "Chapter 4 Fourier Analysis" in "Signal Processing for Communications" by Paolo Prandoni and Martin Vetterli, available from http://www.sp4comm.org/ , and the slides for the entire module 4 of the Digital Signal Processing Coursera course, in https://spark-public.s3.amazonaws.com/dsp/slides/module4-0.pdf .

And if you are registered for this Coursera course, see https://class.coursera.org/dsp-001/wiki/view?page=week3

In this second half of module 4.9, we are going to look specifically

at one example of a Divide and Conquer algorithm, it's the Cooley-Turkey Fast Fourier Transform,

we will see a matrix view of this algorithm, as well as flow-graphs,

and then conclude by looking at more general cases of Fast Fourier Transforms.

So far, we have only had conceptual discussions

of what the Divide and Conquer algorithm would be.

We now need to be concrete and actually develop an algorithm for the Discrete Fourier Transform,

based on the idea of Divide and Conquer.

We'll do this by analyzing an algorithm called the decimation-in-time algorithm.

It'll take us a number of slides until we get to all the details of it.

But this is very important,

because you're going to understand the mechanics of one of the most famous algorithms of all times.

Start with the definition of the DFT, so X[k] is the sum from 0 to N-1

of (?) x[n], multiplied by the root of unity to the power nk.

k goes from 0 to N -1, W is the usual Nth root of unity.

Now, as very often, a good guess is half of the answer to a fast algorithm.

The good guess is the following: break the input into even and odd index terms,

that's why it's called decimation-in-time,

because it's like taking the sequence and throwing out half of the entries,

that's decimation by a factor of 2.

And so, x[n] is mapped into x[2n], the even index samples, and x[2n] + 1, the odd index samples,

and now n goes from 0 to N/2 - 1.

Recall that N is assumed to be a power of 2, so N/2 is again an integer.

Then, we need to do something similar to the output.

But there it's not even and odd, it's first half and second half.

So, X[k], k going from 0 to N minus 1, becomes the first half.

X[k] plus the second half X[k+N/2], where now k runs from 0 to N/2 -1.

Let us now look at the equations for the DFT,

and use this separation of the inputs into even and odd inputs.

So X[k] now is the sum of 2 sums of length N/2,

over the even and odd ones, with respective powers of W.

And we separate these two sums, we simplify the exponent,

and then we need to use a little trick here.

The trick is the following:

If you look at Nth root of unity, you take its power of 2,

then this is equal to e to the -j4π/N, which is equal to e to the -j2π divided by N/2.

And this of course is the root of unity of order N/2.

Using this, we can write W to the power 2nk

as W of root N/2, to the power nk.

We use this and we see now that we have 2 sums of lengths N/2 with respect to the N/2th root of unity

and we write this as the sum of X' k and X'' k,

where X'' k is actually multiplied by W to the k.

Here W is still the nth (?) root of unity.

We can now say this explicitly in the following way:

we have 2 half-size DFTs, which we call X', X''.

We have to multiply the second DFT by W k.

That will require N/2 complex multiplications and we add the results together.

With this, we have found the first half of our result, namely X[k],

for k going from zero to N/2 -1.

Now we need to find out how to compute the second half of the DFT, and then we'll be done.

We have our first half of the DFT: let us consider the second half,

So that's X[k+N/2].

This is made up of two sums over the even and the odd inputs.

And the sum of the odd inputs is pre-multiplied by W to the power k + N/2.

We need again a little trick and this is to simplify W to the power N/2,

this is equal to e to the - j Nπ/N, this is e to the -jπ, which is equal to -1.

With this, we see that

the second sum is simply premultiplied by W to the k and a change of sign,

so we get as a net result that the second half of the DFT is actually equal to X' - Wk X''

So it's the same two parts we had before.

We simply -- subtraction inside of a sum.

In words: we have now computed both parts.

So we divide the input into even and odd samples,

compute 2 DFTs of size N/2.

We multiply this output of the second DFT by the roots of unity of order n (?) raised to the kth power.

That uses capital N/2 multiplications.

Then we recover the outputs by taking the sum and the difference of these 2 parts,

and that will give us X[k] and X[k+N/2].

It is easiest to see this in a diagram, and we do the diagram for a DFT of size 8.

So we have x[0] to x[7].

They are divided into the even parts, the upper half, the odd part, the lower half.

And we compute 2 DFTs of size N/2, in this case,

that would be DFT's of size 4.

Then we simply recombine the respective outputs by sum and differences,

and we get X[0], which is a sum of the first output of both DFTs.

The difference which will give you the X[4] output, which is the difference of these first outputs, etc.

And with this, we have first half, second half.

In this very simple diagram, we see very explicitly this Divide and Conquer algorithm at work.

So what is the complexity now?

So to split the DFT into two pieces of size N/2, that was free.

You simply had to select even and odd index inputs.

Then we compute two DFT's of size N/2.

Recall that the complexity's quadratic, so we have this N square over 2 complexity here.

We merge the two results with multiplications: there are N/2 of them.

And the total complexity, at this stage, is N/2 + N/2 (?),

which, indeed, is smaller than the original complexity, which is quadratic in N.

So, in general, we have about half the complexity of the initial problem.

Now, this would not be a big deal: that's a factor of 2 of savings, nice to have.

But we are aiming for much more than this.

What about repeating this process?

We know that, since N is a power of 2, we can divide in 2, 4, 8, etcetera.

And we can go until we end up with trivial DFTs of size 2.

These simply amount to a sum and a difference.

This will require log N minus 1 (?) steps, because we end up with a DFT of size 2.

Each step requires a merger and uses an order N/2 multiplications and order N additions,

so the total will be N/2 (log-in-base-2 N-1) multiplications and N log-in-base-2 N additions.

So the saving is of order log-in-base-2 N divided by N.

The key result, the "take home" message here, is that a DFT of size N

requires order N log-in-base-2 N operations.

Let us look at a few examples to see how big the savings is.

So when N is equal to 1024, that's 2 to the power of 10,

so the saving is of the order of factor 100.

If N is 2 to the power of 20, which is roughly a million, a little bit more,

then the saving is of the order of 50,000: so that's a big deal, indeed.

We have seen a recursive formulation using sums on even and odd indexed inputs.

We have seen a blog diagram (?) view of the DFT.

Let's look now at a matrix factorization view.

This is important because many implementations, for example in Matlab or FreeMat

will actually use numerical linear algebra,

and then the matrix factorization is actually important.

So in this case we look at W4, the size 4 DFT, so recall you have to separate even and odd samples,

that's the 1st matrix on the right.

Then we compute two DFT's of size 2.

These are sum and differences, it's a block diagonal matrix.

Then we have to multiply the odd outputs by W k:

here, these are simply j's and minus j, and sum and difference, which is included in here,

and we get the overall result, with these three small matrices.

This uses eight additions, no multiplications.

Let's now move to the next bigger size.

This is N=8.

Now this is going to be big because each matrix is of size 8 by 8,

and we have a factorization to take care of,

so the first matrix W8 in its initial form we have seen before.

The first step is to separate even from odd indexed samples.

Call this D 8 for decimation of size 8.

You see the matrix here, you can sort of recognize it's structure,

so the 1st half takes out the even inputs, the 2nd half takes out the odd inputs.

No arithmetic operations so far, just re-indexing of inputs.

The step 2 is to compute two DFTs of size N/2, on these even and odd index samples.

So, each submatrix is a DFT of size 4 and the result is a block diagonal matrix.

On the diagonal, you have the W4 matrices, off the diagonal, you have the 0 matrix.

This requires 2 DFTs of size 4, according to what we have just seen,

that's a total of 16 additions.

Step 3 is to multiply the output of the 2nd DFT by W to the k.

This will be a diagonal matrix with an identity matrix in the 1st half of the diagonal

and the 2nd half is made up of the successive powers of W.

This will require 2 multiplications because W square is actually -j, and that's free.

So we can summarize this as a block matrix

with I4, Λ4 on the diagonal, 04 which are the zero matrices, off the diagonal.

And the Λ4 is the one that contains the complex multiplications by the roots of unity

Step 4: we have to recombine the final output

to obtain X[k] and X[k+N/2] by a sum and difference.

So this matrix is also very structured. The first half computes the sum.

So we can put this in block matrix form by having two identity matrices of size 4.

And the second half computes the differences,

so that's the identity matrix of size 4 minus itself: that requires 8 additions.

So in total, we have the product of 4 matrices, D8, the selector matrix,

then the one that is diagonal with the DFTs of size 4,

then the one that is diagonal with the I4, Λ4, and finally the sum and differences.

So the total will be 24 additions and 2 multiplications.

Now this is much smaller than what we would expect from an N square complexity.

N square would be in this case 64 order, 64 additions and multiplications.

So you see that the big saving here is multiplications.

That is the case here because the DFT is very small

and so it quickly comes down to trivial DFT's.

We can write out what we just did in an explicit flow-graph.

Such a flow-graph simply indicates how inputs -- x[0] to x[7| -- go through the flow-graph,

get multiplied by the roots of unity when needed, then take sum and differences etc.

And as you go from left to right, you end up with the output X[0] to X[7].

What you can recognize visually are the butterflies, which are simply the sum and differences

that appear again and again in the recursive computation of the DFT

using the Fast Fourier Transform algorithm.

Now, is this a big deal?

In image processing, so, for example, on your digital camera or on your mobile phone,

the image is taken in small blocks of 8 by 8.

One computes a transform that is called a Discrete Cosine Transform,

we shall see this later in this class

when we study the JPEG compression algorithm,

and the DCT is a little bit similar to a DFT, it's like a real DFT.

It has a fast algorithm that is very much inspired by what we just saw.

Let's compare direct multiplication: so 8 by 8 is 64 inputs;

squared, that's 4000 multiplications; it's a dominant cost because it's, it's fixed point multiplications.

The transform can be computed in rows and columns separately,

so you have 16 DFT's, essentially.

We have seen that such an algorithm would take 2 multiplications,

so we would have a total of 32 multiplications.

So we have a saving of the order of 2 orders of magnitude.

And that's a big deal.

So, instead of compressing your image in 100 seconds, it will take 1 second, for example.

So far, we have only considered lengths N which were powers of 2.

But of course, other lengths are also important,

and the message here is that for every possible length N, there is a fast algorithm.

So the first case is N is equal to a power of a prime

We looked (?) at the case p is equal to 2, and we will derive the Cooley-Tukey algorithm.

Very famous, very simple and very important.

When N is a product of 2 coprime numbers,

then there is a way to reorder this into a two-dimensional DFT of size N1 by N2,

this is known as the Good-Thomas algorithm, it was invented in the 50s

and reinvented in the 60s again.

When N is a prime, we can do a mapping of the Transform

that results in a convolution of size N-1.

This might not be good news,

except that we'll see that the convolution of size N-1 can also be solved with a DFT,

and so if N is, for example, 17, then N-1 is 16,

and this we can solve with a usual Cooley-Tukey fast algorithm.

By combining all the tricks involved,

minimizing the number of multiplications, reordering the terms,

we have the so-called Winograd Fast Fourier Transform,

which was invented by Winograd in the late 1970's,

and which really caps the work on fast algorithm development from a theoretical point of view.

Let's us look at this case N=5 that we had seen before: the matrix is familiar.

We take out one piece of it and we remove the columns of 1 and the row of 1,

and we are left with a matrix that has, of course, all powers of W from 1 to 5,

but each power appears only once

-- that's the important thing here -- and it can be shown that if you reorder rows and columns properly,

you end up with this matrix that we see on the right side which is actually a circulant matrix.

So you can see that the 2nd row is simply the 1st row,

circularly shifted, and so on, to the bottom.

And such a circulant matrix, we will see, can be diagonalized by a DFT of size 4,

so we are back to computing a DFT but now of a size 1 less.

When N is a product of 2 coprime numbers.

Here, we take N=6.

It can be written as a 2-dimensional DFT of size 3 by 2.

And such a DFT can be computed by having 3 DFTs of size 2, and 2 DFTs of size 3.

We write this here in a factorized form,

so we have permutation matrices R, Q and P and we have block diagonal matrices.

The first one on the left is block diagonal with DF-- three DFTs of size 2.

The block diagonal matrix on the right has two DFTs of size three.

Except for permutations, there was no cost in transforming a DFT of size 6

into 3 DFTs of size 2, and 2 DFTs of size 3, so this is very efficient.

The conclusion is don't worry be happy:

any length of DFT can be mapped into a Fast Fourier Transform algorithm.

Therefore this is an important message because you don't have to take your signal,

zero-pad so as to get a length that is a power of 2, for example.

This is often done but it's a mistake, if you really want to compute a DFT of size 17,

then there is an algorithm that will compute it and it will compute it fast.

I also want to say that there are extremely good packages out there,

something called, for example, the Fastest Fourier Transform in the West, or the SPIRAL program.

And these packages allow you to not worry,

and simply find the best possible Fourier Transform for any length.

And as we had seen before, this makes a big difference.

We told-- we are talking about orders of magnitude in speed up.

Let me conclude this module on the Fourier Transform.

The Fourier transform is a magical tool.

It's a beautiful mathematical tool, it is also a very practical tool.

We have seen the DTFD, which is good for math.

We have seen the DFT, which is good for computation.

And we have seen the Fast Fourier Transform,

which allows to compute Fourier transforms very rapidly.

This fascination with Fourier can be seen here in a tribute on a license plate.