fft

COT 4401 Top 10 Algorithms
Chris Lacher
Fast Fourier Transform [FFT]

Resources

Chapter 30 of [Cormen] provides an excellent discussion with background on polynomials, complex roots of unity, and a mathematical (application neutral) approach.
Maxim Raginsky's lecture notes present a more direct discussion from the EE Signal Analysis perspective.
FFT played a role in finding bug in Intel Skylake processor.
Here's an "English" intro: Smoothies Unmixed
There are many other sources for both code and applications available.

Background (Polynomials)

We'll loosely follow [Cormen] and use their terminology. A neat mathematical fact is their Theorem 30.1:

Theorem 1 (Uniqueness of interpolating polynomial).
For any set {(x₀,y₀), (x₁,y₁), ..., (x_n,y_n)} of number pairs, with x₀, x₁, ... , x_n all distinct, there is a unique polynomial A(x) of degree <= n such that y_i = A(x_i) for i = 0, ..., n.

To get a feel for this result, look at some special cases:

[n = 1]: we have two points (x₀,y₀) and (x₁,y₁), with x₀ != x₁. There is a unique straight line passing through these two points, and because x₀ != x₁ this line is not vertical. Any such line is given by an equation of the form y = mx + b, where m is the slope of the line and b is the y-intercept. A(x) = mx + b is the interpolating polynomial. (Exercise: find formulas expressing m and b in terms of the point coordinates.)

[n = 2]: we have three points with distinct x coordinates. There is a unique parabola (a straight line in the degenerate case where the points are colinear) passing through these points, with an equation of the form y = ax² + bx + c. The interpolating polynomial is A(x) = ax² + bx + c.

This theorem is true even when the pairs are of complex numbers! We might then have complex number coefficients for the interpolating polynomial, but when the pairs of complex numbers have certain symmetry properties the interpolating polynomial will have real coefficients.

A companion result is

Theorem 2 (Factorization of polynomials).
Any polynomial A(x) with real coefficients and degree n can be factored as a product of degree 1 (linear) and degree 2 (quadratic) polynomials with real coefficients:

A(x) = c(x - r₁) ... (x - r_k)Q₁(x) ... Q_m(x)

where c is the leading coefficient of A(x), r₁, ..., r_k are the real roots of A(x), and each Q(x) has two non-real roots that are conjugate pairs of complex numbers. (Note that n = k + 2m.) We can find the roots of the quadratic factors using the "quadratic formula". The factorization shown above is unique: any two such factorizations differ only in the order in which we write the factors.

Example 1:

x⁴ - 5x³ + 10x² -10x + 4 = (x - 1)(x - 2)(x² - 2x + 2)

= (x - 1)(x - 2)(x - 1 - i)(x - 1 + i)

roots: 1, 2, 1+i, 1-i

Example 2:

x⁵ - 6x⁴ + 12x³ - 12x² +11x - 6 = (x - 1)(x - 2)(x - 3)(x² + 1)

= (x - 1)(x - 2)(x - 3)(x - i)(x + i)

roots: 1, 2, 3, i, -i

where i² = -1. (Recall that the conjugate of the complex number a + bi is a - bi and that (a + bi)(a - bi) = a² + b².)

These two facts together tell us that a polynomial of degree n is completely determined by its values on (n + 1) distinct points. Thus to specify a polynomial of degree n or smaller we need only specify its values on n+1 points.

About proofs

Theorem 1 has a straightforward verification that is constructive and thus yields an algorithm for calculating A(x) from a set of point values {(x₀,y₀), (x₁,y₁), ..., (x_n,y_n)}. As shown in [Cormen, p. 903], the algorithm reduces to calculating the inverse of the so-named Vandermonde Matrix

1 x₀ x₀² x₀³ .  .  .   x₀ⁿ

1 x₁ x₁² x₁³ .  .  .   x₁ⁿ

. . . . .  .  .   .

1 x_n x_n² x_n³ .  .  .   x_nⁿ

Note that inverting an n x n matrix is a Θ(n³) operation. A faster algorithm for finding A(x) from its point values uses Lagrange's Formula (see [Cormen]) but is still Θ(n²). The FFT solves this problem in time O(n log n), given a judicious choice for the points x_i.

The proof of Theorem 2 is significantly more involved and has a long history. (The result for complex polynomials is often called the Fundamental Theorem of Algebra.) There were conjectures, claims of proof, and claims of counterexample throughout the period 1608 - 1920. Gauss is generally credited with the first correct proof in 1799, but even that has a serious "gap" that was finally filled by Ostrowski in 1920. No "elementary" proof exists, nor is there a practical constructive proof.

Framing FFT as a multiplication of polynomials problem.

First let's agree on some terminology. A polynomial A(x) is normally given in coefficient representation:

A(x) = a₀ + a₁x + a₂x² + . . . + a_nxⁿ

The theory above makes it clear that there is an equivalent point-value representation

{(x₀,y₀), (x₁,y₁), ... , (x_n,y_n)}, y_i = A(x_i), or
{(x₀,A(x₀)), (x₁,A(x₁)), ... , (x_n,A(x_n))}

for any set of n+1 distinct numbers x₀, x₁, ... , x_n.

Now suppose we have two polynomials A(x) and B(x) of degree n and we need to multiply them. Multiplying using the standard coefficient representation is straightforward and learned in basic algebra. The degree k term in the product is the sum-of-products

Σ_i+j=k a_ib_j

from which we can see that calculating the product polyniomial A(x)B(x) using the coefficient representations is a Θ(n²) process. Our goal is to use three steps to reduce this to O(n log n):

Translate A(x) and B(x) to point-value representations. [FFT - O(n log n)]
Multiply the point-value representations to obtain a point-value representation of the product. [straightforward evaluation - O(n)]
Translate the product point-value representation to coefficient form. [FFT - O(n log n)]

Roots of Unity

You are all familier with the "imaginary" number i that is typically (and not quite correctly) described as "the square root of -1". The incorrectness of that statement is due to the fact there are two such complex numbers - two solutions to the equation

x² = -1 | equivalently x² + 1 = 0

We generally assign i to one of these, and the other is then -i. Note that if ω = i or -i then ω² = -1 and that ω⁴ = 1: i and -i are 4th roots of 1. There are 4 4th roots of 1 - the other two are 1 and -1.

More generally, a complex nth root of unity is a complex number such that

ωⁿ = 1

There are exactly n such roots, and we know a formula to calculate them:

e^2πik/n , for k = 0, 1, 2, ..., n - 1.

(We are using the greek letter π = "pi" in this formula. Unfortunately it seems to render very much like n.) Continuing with the complex number revelations, there is a trigonometric representation for complex numbers, from which we can derive

e^iu = cos(u) + i sin(u) for any real number u

This last formula allows us to interpret the nth roots of unity as n points on the unit circle with angles that are integer multiples of 2π/n: n points on the unit circle that are mutually equidistant from each other and include 1 = 1 + 0i (the case with angle 0 or 2π).

Following [Cormen] we define the principal nth root of unity to be

ω_n = e^2πi/n

(the case k = 1).

Exercise. Show that the powers ω_n⁰ , ω_n¹ , ... , ω_n^n-1 of ω_n are all of the nth roots of unity.

The complex unit circle is the set of all complex numbers with modulus 1, that is:

{z = x + yi | x and y are real and x² + y² = 1}

A very handy property of the complex unit circle is that it is closed under multiplication: the product of two numbers z and w on the unit circle is also on the unit circle. The trigonometric representation of unit circle numbers gives deeper insight into how the unit circle behaves under multiplication. Note that any number z on the unit circle is determined by its angle θ in the representation

z = e^{i θ} = cos(θ) + i sin(θ)

where θ is the angle z makes with the x axis.

Exercise. Show that if z and w are on the unit circle then the angle of the product zw is the sum of the angles of z and w.

Exercise. Show that

ω_dn^dk = ω_n^k

for any non-negative integers d, k, n. That is, the dnth principal root of unity raised to the power dk is equal to the principal nth root of unity raised to the power k

How FFT works

We will rephrase FFT below, but it is worth looking at the essential cleverness used to find the point-value representation, given:

A polynomial A(x) of degree less than n
n is a power of 2
The points are the nth roots of unity

Note that we can write A(x) as a sum of even power terms followed by a sum of odd power terms:

A(x) = (a₀ + a₂x² + ... + a_n-2x^n-2) + (a₁x + a₃x³ + ... + a_n-1x^n-1)

Then we have

A(x) = B(x²) + x C(x²)

where

B(x) = a₀ + a₂x + ... + a_n-2x^{n/2 - 2}

and

C(x) = a₁ + a₃x + ... + a_n-1x^{n/2 - 2}

Evaluating at the nth roots of unity, we have

A(ω_n^k) = B(ω_n^2k) + ω_n^k C(ω_n^2k) = B(ω_n/2^k) + ω_n^k C(ω_n/2^k)

(the last equality using the Exercise on powers of roots of unity and the assumption that n is even.) We have expressed the values of A(x) on the nth roots of unity in terms of values of B(x) and C(x) on the n/2 th roots of unity. Note also that B(x) and C(x) have degree less than n/2, so we have reduced the original calculation into two such calculations 1/2 the size of the original. This leads to a recursive calculation with each stage requiring O(n) steps, and the depth of the recursion is log₂ n. Hence we have the ingredients of an O(n log n) algorithm to obtain the point-value representation, given assumptions 1-3.

DFT

Fourier Analysis refers to the broad and vital topic of expressing functions in terms of their periodic components. These functions may operate on one or several variables that may represent time, space, space-time, or other more exotic domains. The periodic components would represent such things as pure tones (when the function values represent sound) or pure colors of light (when the function values represent a picture). A Fourier Transform is a specific mapping from functions of one value type to another, expressing a function as a sum (or integral) of periodic components.

Consider as an example a function f(t) representing a sound over some interval of time [0, t₁). A Fourier transform of f would be a function g(θ) representing the periodic components of the sound as a weighted sum (or integral) of pure tones (the "sound spectrum").

Another example is a function f(x,y) representing representing a color value over some rectangle in 2-D space. A Fourier transform of f would be a function g(θ) representing the periodic components of the color as a weighted sum (or integral) of pure colors (the "RGB" representation).

The study of Fourier analysis is a huge and still very active field, with applications as far-reaching as power grid stability, digital sound representation and filtering, and compression of digital images. (The Joint Photographic Experts Group, or JPEG, uses Fourier analysis in its image compression algorithms.) Obviously, detailed delving into any of these areas is beyond the scope of this class. Just as obvious, any algorithm that improves the computational complexity of Fourier analysis, especially of digital data, has vital impact throughout our modern world. So with the backdrop of our specific polynomial multiplication problem and the far-reaching digital signal processing field, we now lay out the digital Fourier transform and its improvement, the fast Fourier transform.

We will specialize from here on to the case of discrete-time functions of one variable. These values might be digital samplings from a continuous function, or they might be values in an inherently discrete system. (Or, they may be the coefficients of a polynomial of degree less than n.) Whatever their origin, assume we have a vector

a = (a₀, a₁, ... , a_n-1)

of real or complex components a_k. The discrete Fourier transform [DFT] of a is the vector

y = (y₀, y₁, ... , y_n-1)

where

y_k = Σ_j=0^n-1 a_jω_n^kj = a₀ + a₁ω_n^k + a₂ω_n^2k + ... + a_n-1ω_n^k(n-1)

(where ω_n is the principal nth root of unity). We use the notation

y = DFT_n(a)

to denote our digital Fourier transform.

Theorem 3. DFT_n(a) is Θ(n²).

FFT

The formula defining y_k above requires n multiplications of complex numbers, and hence DFT_n(a) requires Θ(n²) multiplications. Considering that, for example in the case of digital representation of sound or color, n may be very large, this quadratic complexity can be a barrier to useful computation of DFT_n, particularly in real time. The fast Fourier transform [FFT] is a method for computing DFT_n in time Θ(n log n) under the assumption that n is a power of 2. (This last assumption is not as constraining as it might appear. In applications where the assumption is inconvenient, the vector a can be "padded" with zero values out to the next power of 2 before applying FFT.)

The FFT follows the calculation we introduced above for converting polynomials to point-value representation. This specific version is from [Cormen]:

FFT (a, n)  // a is an array of size n (or smaller)
{
  if (n <= 1) return a;             // base case for recursion
  wn = exp(2*pi*i/n);               // wn = principal nth root of unity
  w  = 1;                           // w  is maintained as wn^k, updated in post-processing loop
  b  = (a[0], a[2], ... , a[n-2]);  // even indexed values
  c  = (a[1], a[3], ... , a[n-1]);  // odd indexed values
  yb = FFT(b, n/2);                 // recursive call
  yc = FFT(c, n/2);                 // recursive call
  for (k = 0; k < n/2; ++k)
  {                               // loop invariant: w = wn^k at beginning of loop body
    ya[k]     = yb[k] + w*yc[k];  // see (1) below
    ya[k+n/2] = yb[k] - w*yc[k];  // see (2) below
    w = w*wn;                     // w = wn^(k+1) for next iteration
  }
  return ya;
}

Proof (1): show that yb[k] + w*yc[k] is equal to ya[k] for 0 ≤ k < n/2

Looking separately at the two terms of the sum, we have:

yb_k = Σ_j=0^n/2-1 b_j (ω_n/2)^kj (inductive hypothesis/recursive call)

= Σ a_2j (ω_n/2)^kj

= Σ a_2j (ω_n)^2kj

(ω_n)^k yc_k = (ω_n)^k (Σ_j=0^n/2-1 c_j (ω_n/2)^kj)

= (ω_n)^k (Σ a_2j+1 (ω_n/2)^kj)

= (ω_n)^k(Σ a_2j+1 (ω_n)^2kj)

= Σ a_2j+1 (ω_n)^k+2kj

= Σ a_2j+1 (ω_n)^k(2j+1)

Therefore yb_k + (ω_n)^k yc_k = Σ_{j even} a_j (ω_n)^kj + Σ_{j odd} a_j (ω_n)^kj = ya_k

Proof (2): show that yb[k] - w*yc[k] is equal to ya[k+n/2] for 0 ≤ k < n/2

First observe that (ω_n)^k+(n/2) = - (ω_n)^k: since n is even, (ω_n)^n/2 = -1 (half way around the unit circle), so (ω_n)^k+(n/2) = (ω_n)^k(ω_n)^n/2 = -(ω_n)^k. After making this substitution, the remainder of the proof follows like that for (1).

Theorem 4. FFT_n(a) is Θ(n log n).

Conclusion

Exercise. Write a short obituary of Fourier. Make his family proud in a few paragraphs.

Exercise. Suppose you are on a job interview at Google, and you are asked to describe the FFT and why it is important. What is your answer?

Exercise. Give an informal argument that FFT_n has runtime Θ(n log n).

Exercise. Who is generally credited with discovery of FFT? Whose work does that discovery build on (other than Fourier, of course).

index↑

x⁴ - 5x³ + 10x² -10x + 4	= (x - 1)(x - 2)(x² - 2x + 2)
	= (x - 1)(x - 2)(x - 1 - i)(x - 1 + i)
	roots: 1, 2, 1+i, 1-i

x⁵ - 6x⁴ + 12x³ - 12x² +11x - 6	= (x - 1)(x - 2)(x - 3)(x² + 1)
	= (x - 1)(x - 2)(x - 3)(x - i)(x + i)
	roots: 1, 2, 3, i, -i

1	x₀	x₀²	x₀³	. . .	x₀ⁿ
1	x₁	x₁²	x₁³	. . .	x₁ⁿ
.	.	.	.	. . .	.
1	x_n	x_n²	x_n³	. . .	x_nⁿ

yb_k	= Σ_j=0^n/2-1 b_j (ω_n/2)^kj (inductive hypothesis/recursive call)
	= Σ a_2j (ω_n/2)^kj
	= Σ a_2j (ω_n)^2kj
(ω_n)^k yc_k	= (ω_n)^k (Σ_j=0^n/2-1 c_j (ω_n/2)^kj)
	= (ω_n)^k (Σ a_2j+1 (ω_n/2)^kj)
	= (ω_n)^k(Σ a_2j+1 (ω_n)^2kj)
	= Σ a_2j+1 (ω_n)^k+2kj
	= Σ a_2j+1 (ω_n)^k(2j+1)