31 A5: Fast Fourier Transform

Assignment 5: Fast Fourier Transform

Published

April 6, 2026

Note

Complete the following and submit to Canvas before April 15, 23:59
Late work will recieve 0%,
Each assignment is worth the same,
Please get in contact with the TA/instructor in plenty of time if you need help,
Before submitting your work, make sure to check everything runs as expected. Click Kernel > Restart Kernel and Run All Cells.
Feel free to add more cells to experiment or test your answers,
I encourage you to discuss the course material and assignment questions with your classmates. However, unless otherwise explicitly stated on the assignment, you must complete and write up your solutions on your own,
The use of GenAI is prohibited as outlined in the course syllabus. If I suspect you of cheating, you may be asked to complete a written or oral exam on the content of this assignment.

using LinearAlgebra, Plots, LaTeXStrings

Name:

Approx. time spent on this assignment:

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.

Your answer here

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)

Your answer here

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^8
h = 1/n
X = 0:h:1-h

# write down the matrix ℱ (exercise 1)
ℱ = # your code here
Fhat = ℱ*f.(X)

# write down the Fourier coeficients 
# of the numerical solution (exercise 2)
Uhat = Fhat ./ # your code here

# leave the following as it is 
# (in Julia A' is the conjugate transpose of A)
U = ℱ' * Uhat 
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?

Your answer here

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.

Your answer here

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 
            # your answer here:

            

        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) = # your code here

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

1.466075386898995e-15

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

Other Links