37  A5: Solutions

Assignment 5: Fast Fourier Transform

using LinearAlgebra, Plots, LaTeXStrings

Recall that we are interested in solving the periodic problem

\begin{align} &-u'' + p u' + q u = f \qquad \text{on } (0,1) \nonumber\\ % &u(0) = u(1) \nonumber\\ &u'(0) = u'(1). \end{align}

We will take p = 0 and q > 0. We approximate this problem with \mathcal L_h U = F where U_0 = U_n,

\begin{align} (\mathcal L_h U)_j := -\frac{U_{j-1} - 2 U_j + U_{j+1}}{h^2} % + q U_j, \end{align}

and F_j := f(x_j). By solving the linear system of equations \mathcal L_h U = F, we find an approximation to the continuous problem.

Recall that we may write U_j = \sum_{k=0}^{n-1} \hat{U}_k (\bm e_k)_j where (\bm e_k)_j := \frac{1}{\sqrt{n}} e^{\frac{2\pi i j k}{n}} is the Fourier basis (note that you will see different normalisations in different texts) and \hat{U}_k := \sum_{j=0}^{n-1} \overline{ (e_{\bm k})_j } U_j.

Exercise 1. Write \mathcal L_h U = F in terms of the Fourier coefficients (\hat{U}_k, \hat{F}_k) and hence explain how one may solve the linear system \mathcal L_h U = F by implementing the discrete Fourier transform and the inverse discrete Fourier transform.

Notice that (\bm e_k)_{j\pm1} = \frac{1}{\sqrt{n}} e^{\frac{2\pi i (j\pm1)k}{n}} = e^{\pm\frac{2\pi i k}{n}} (\bm e_k)_j and so

\begin{align} (\mathcal L_h \bm e_k)_{j} % &= \left( - \frac{e^{-\frac{2\pi i k}{n}} - 2 + e^{\frac{2\pi i k}{n}}}{h^2} + q \right) (\bm e_k)_{j} \nonumber\\ % &= \left(\frac{2}{h^2} (1 - \cos \tfrac{2\pi k}{n}) + q \right) (\bm e_k)_{j} \nonumber\\ % &=: \lambda_k (\bm e_k)_{j}. \end{align}

As a result, we have

\begin{align} \mathcal L_h U &= \sum_{k=0}^{n-1} \lambda_k \hat{U}_{k} \bm e_k \nonumber\\ % F &= \sum_{k=0}^{n-1} \hat{F}_k \bm e_k. \end{align}

Therefore, we can (i) compute the Fourier transform \hat{F}_k, (ii) define \hat{U}_k := \hat{F}_k / \lambda_k, and (iii) compute the inverse Fourier transform U_j.

Exercise 2. Write down a matrix \mathscr F such that \hat{U} = \mathscr F U and U = \overline{\mathscr F}^T \hat{U}.

(The conjugate transpose is normally written \mathscr F^\star := \overline{\mathscr F}^T)

Since \hat{U}_{k} = \sum_{j=0}^{n-1} \overline{(\bm e_k)_j} U_j, U_{j} = \sum_{k=0}^{n-1} (\bm e_k)_j \hat{U}_k, and (\bm e_k)_j = (\bm e_j)_k so we can define \mathscr F_{kj} := \overline{(\bm e_k)_j}.

Exercise 3. Take f(x) = x^2 and q = 1. Use your previous two answers to fill in the gaps in the following code sample to approximate the solution to -u''(x) + u(x) = x^2 with periodic boundary conditions.

q = 1
f(x) = x^2 
n = 2^12
h = 1/n
X = 0:h:1-h

t = @elapsed begin
    ℱ = [ (1/√n)*exp(-2π*1im*j*k/n) for j∈0:n-1, k∈0:n-1 ] 
    Fhat = ℱ*f.(X)
    Uhat = Fhat ./ [ (2/h^2)*(1 - cos(2π*k/n)) + q for k∈0:n-1 ] 
    U = ℱ' * Uhat 
end

U = [U..., U[1]]
plot( 0:h:1, real.(U); 
    m=1, legend=false,
    title=L"Numerical solution to $-u''(x) + u(x) = x^2$ with PBCs" )

Exercise 4. What is the computational cost involved in solving \mathcal L_h U = F by using the code above?

Here, we are multiplying dense matrices with vectors which is O(n^2) where n is the number of grid points.

A Fast Fourier Transform (FFT) refers to any algorithm that improves on this computational complexity. In order to get the main idea, it is useful to first consider the case where n = 2^m is a power of two. Let

\begin{align} E_k &:= \frac{1}{\sqrt{n}} \sum_{j=0}^{\frac{n}{2}-1} U_{2j} e^{-\frac{2\pi ijk}{n/2}}, \nonumber\\ % O_k &:= \frac{1}{\sqrt{n}} \sum_{j=0}^{\frac{n}{2}-1} U_{2j+1} e^{-\frac{2\pi ijk}{n/2}} \nonumber \end{align}

be (up to a constant) the Discrete Fourier Transform of the even and odd terms, respectively.

Exercise 5. Show that

\begin{align} \hat{U}_k &= E_k + e^{-\frac{2\pi i k}{n}} O_k \nonumber\\ \hat{U}_{k+\frac{n}2} &= E_k - e^{-\frac{2\pi i k}{n}} O_k \end{align}

for all k =0,\dots,\frac{n}{2}-1.

We have

\begin{align} \hat{U}_k &= \frac{1}{\sqrt{n}} \sum_{j=0}^{n-1} U_j e^{-\frac{2\pi i jk}{n}} \nonumber\\ % &= \frac{1}{\sqrt{n}} \sum_{j=0}^{\frac{n}{2}-1}\left[ U_{2j} e^{-\frac{2\pi i (2j) k}{n}} % + U_{2j+1} e^{-\frac{2\pi i (2j+1) k}{n}} \right] \nonumber\\ % &= \frac{1}{\sqrt{n}} \sum_{j=0}^{\frac{n}{2}-1}\left[ U_{2j} e^{-\frac{2\pi i j k}{n/2}} % + e^{-\frac{2\pi i k}{n}} U_{2j+1} e^{-\frac{2\pi i jk}{n/2}} \right] \nonumber\\ % &= E_k + e^{-\frac{2\pi i k}{n}} O_k. \end{align}

Moreover,

\begin{align} \hat{U}_{k+\frac{n}{2}} &= E_{k+\frac{n}{2}} + e^{-\frac{2\pi i (k+\frac{n}{2})}{n}} O_{k + \frac{n}{2}} \nonumber\\ % &= E_{k+\frac{n}{2}} - e^{-\frac{2\pi i k}{n}} O_{k + \frac{n}{2}} \end{align}

and E_{k+\frac{n}{2}} = E_{k} and O_{k + \frac{n}{2}} = O_k since e^{-\frac{2\pi ij(k+\frac{n}{2})}{n/2}} = e^{-2\pi ij}e^{-\frac{2\pi ijk}{n/2}} = e^{-\frac{2\pi ijk}{n/2}}.

This is the basis for the (radix-2) Cooley-Tukey algorithm where one computes the Fourier transform on \{0,1,\dots,n-1\} by iteratively subdividing the domain into two parts and computing (rescaled) Fourier transforms on each smaller piece (this is an example of a divide-and-conquer algorithm).

Exercise 6. When n=1, what is the discrete Fourier transform of \{U_0\}? Hence, complete the following code sample:

function FFT(u; normalise=true)
    n = length(u)
    if !ispow2(n)
        @warn "n = $n is not a power of two"
    else
        U = zeros(ComplexF64,n)
        if n == 1 
            U[1] = u[1]
        else 
            E = FFT(u[1:2:end]; normalise=false)
            O = FFT(u[2:2:end]; normalise=false)
            for k ∈ 1:Int(n/2)
                t = exp(-2π*1im*(k-1)/n) # Twiddle factor
                U[k]          = E[k] + t * O[k]
                U[Int(k+n/2)] = E[k] - t * O[k]
            end
        end
        return normalise ? U / √n : U
    end 
end

FFT (generic function with 1 method)

Exercise 7. By using the function FFT() and the complex conjugate conj(), implement a one-line function that computes the inverse discrete Fourier transform:

\begin{align} U_j = \frac{1}{\sqrt{n}} \sum_{k=0}^{n-1} \hat{U}_k e^{\frac{2\pi i jk}{n}}. \end{align}

iFFT(u) = conj( FFT(conj(u)) )

# this should be (close to) zero:
norm( iFFT( FFT( f.(X) ) ) - f.(X) )

7.516995118346703e-15

Here, we again solve -u'' + u = x^2 using our FFT method:

t_FFT = @elapsed begin 
    Fhat = FFT( f.(X) )
    Uhat = Fhat ./ [ (2/h^2)*(1 - cos(2π*k/n)) + q for k∈0:n-1 ] 
    U = iFFT( Uhat ) 
end
U = [U..., U[1]]
plot( 0:h:1, real.(U); 
    m=1, legend=false,
    title=L"Numerical solution to $-u''(x) + u(x) = x^2$ with PBCs" )

println("this is $(round(Int,t/t_FFT)) × faster than before!")

this is 174 × faster than before!

Exercise 8. Look up another use for discrete Fourier transform in the literature. Summarise the key ideas in a few sentences.

There’s loads of potential answers here! e.g. look at Brigham (1988).

e.g. 

• signal & audio processing: MP3, filtering (removing noise), Shazam!
• image compression (JPEG), deblurring, pattern recognition,
• solving PDEs - physics, chemistry, ….
• multiplying extremely large numbers/polynomials,
• medical imaging (e.g. MRI, CT),
• fingerprint recognition, ……

This is only possible because FFT are O(n \log n) rather than the naive O(n^2) (multiplying dense matrices): for example spotify uses a sampling rate of 44,100 Hz - one second of audio has n = 44,100 samples

n = 44100
println("n^2  = ", n^2) 
println("n log(n) = ", round(Int, n*log(n)))  

n^2  = 1944810000
n log(n) = 471615

470,000 vs 1.9 billion is a big saving!

Brigham, E. Oran. 1988. The Fast Fourier Transform and Its Applications. Englewood Cliffs, NJ: Prentice Hall.