The Fast Fourier Transform Algorithm

December 26, 2015

The last section on the divide-and-conquer design paradigm on the book Algorithms by S. Dasgupta, C. H. Papadimitriou, and U. V. Vazirani focuses on a very interesting and well known algorithm: The Fast Fourier Transform algorithm. The Fast Fourier Transform is an algorithm which takes a coefficient representation of a polynomial and changes it to its equivalent point-wise representation. It is widely used in a varieaty of tasks in computer science such as: getting rid of noise in data, pattern matching, digital filtering, and others.

The big idea in the Fast Fourier Transform algorithm is that a change in the mathematical representation in a problem can simplify and have efficiency benefits over others. To better undesrtand this idea, let’s examine the following polynomial: f(x) = 5 + 2x + x^2. This poynomial has degree 2, and 3 numbers represent it (its coefficients). In general, a polynomial of n - 1 degree has n coefficients. Now, let’s investigate what values does this polynomial has at x equals to 0, 1, and 2.

  x = 0, then f(x) = 5
  x = 1, then f(x) = 8
  x = 2, then f(x) = 13

The natural question is then, if n points are given, do these numbers define a polynomial of n - 1 degree? Yes, we can solve a set of linear equations, and the answer will be equal to f(x) = 5 + 2x + x^2. So what does this mean? We have been able to show that a polynomial can be represented in two different ways:

A coefficient representation
A point-wise representation

The algorithmic steps of the Fast Fourier Transform algorithm takes (assume n is a power of 2):

  A(x) = a0 + a1x + a2x^2 + ... + a(n-1)x(n-1)

Re-writes it as:

  a0      + a2x^2       + a4x^4        + a(n-2)x^(n-2)
    + a1x        + a3x^3       + a4x^4 + a(n-1)x^(n-1)

Then, define new polynomials (the divide step):

  Ae(y) = a0 + a2y + a4y^2 + a6^y3 + ... + a(n-2)y^(n-2/2)
  Ao(y) = a1 + a3y + a5y^2 + a6^y3 + ... + a(n-1)y^(n-2/2)

And finally, A(w) can be expressed as (the combine step):

  A(w) = Ae(w^2) + w * Ao(w^2)

The only issue though, is that the points which the algorithm needs to use to compute A(w) for an n (power of 2) need to be carefully chosen. The method the author suggests is using the complex nth-roots of unity. By using the complex nth-roots of unity for an n power of 2, then at successive levels of recursion tree we will have (n/2^k)th roots of unity for k=0, 1, 2, 3, ....

The Fast Fourier Transform Algorithm

Pseudocode: (Page 71)

function FFT(A, ω)
  Input: Coefficient representation of a polynomial A(x) of degree ≤ n − 1, where n is a power of 2
  Output: Value representation A(ω^0), . . . , A(ω^n−1)

  if ω = 1: return A(1)
  express A(x) in the form Ae(x^2) + xAo(x^2)
  call FFT(Ae, ω^2) to evaluate Ae at even powers of ω
  call FFT(Ao, ω^2) to evaluate Ao at even powers of ω
  for j = 0 to n − 1:
    compute A(ω^j) = Ae(ω^2j) + ω^jAo(ω^2j)
  return A(ω^0), . . . , A(ω^n−1)

Implementation:

class FFT
  # Input: n coefficients
  # Output: Point-wise representation of the n coefficients
  # Vec size must be a power of 2
  def fft(vec)
    return vec if vec.size <= 1
    # Split A(x) into its odd and even powers
    a_even = vec.each_with_index.select { |_, i| i.even? }.map { |i, _| i }
    a_odd = vec.each_with_index.select { |_, i| i.odd? }.map { |i, _| i }

    # Express A(x) in the form Ae(x^2) + xAo(x^2)
    fft_even = fft(a_even) * 2
    fft_odd  = fft(a_odd) * 2

    # Compute Ae(x^2) + xAo(x^2)
    fft_even.zip(fft_odd).each_with_index.map { |(even, odd), i|
      even + odd * omega(i, vec.size)
    }
  end

  private

  # Calculates (e ^ (2πik/n))
  def omega(i, n)
    Math::E ** Complex(0, -2 * Math::PI * i / n)
  end
end

Tests:

require 'minitest/autorun'
require './fft'

describe FFT do

  before do
    @fft = FFT.new
  end

  describe "#fft" do
    it "should return a point-wise representation of a polynomial of degree 1" do
      assert_equal(@fft.fft([1, 1]), [
        Complex(2, +0),
        Complex(0.0, -1.2246467991473532e-16)
      ])
    end

    it "should return a point-wise representation of a polynomial of degree 3" do
      assert_equal(@fft.fft([1, 1, 0, 0]), [
        Complex(2, +0),
        Complex(1.0, -1.0),
        Complex(0.0, -1.2246467991473532e-16),
        Complex(0.9999999999999998, +1.0)
      ])
    end

    it "should return a point-wise representation of a polynomial of degree 7" do
      assert_equal(@fft.fft([1, 1, 1, 1, 0, 0, 0, 0]), [
        Complex(4, +0),
        Complex(1.0, -2.414213562373095),
        Complex(-1.2246467991473532e-16, -1.2246467991473532e-16),
        Complex(1.0, -0.4142135623730949),
        Complex(0.0, -2.4492935982947064e-16),
        Complex(0.9999999999999998, +0.41421356237309515),
        Complex(1.2246467991473532e-16, -1.224646799147353e-16),
        Complex(0.9999999999999994, +2.414213562373095)
      ])
    end
  end
end

Conclusion

The Fast Fourier Transform algorithm is a really nice way to demonstrate how mathematical ingenuity plays a big role in the design and analysis of algorithms. Changing the representation of the data given to us for a specific task can bring great benefits to both the understanding of the overall process and the efficiency in the computations needed to produce the desired output.