The Fast Fourier Transform ========================== The :index:`fast Fourier transform`, abbreviated by :index:`FFT`, computes the discrete Fourier transform for :math:`n` points in time :math:`O(n \log(n))`. Time and Frequency Domains -------------------------- Data is often periodic, e.g.: seasonal temperature fluctuations. To fit periodic data, a basis of cosines and sines works well. The problem we solve is an :index:`interpolation` problem, stated as follows. Given :math:`(t_i, f_i)`, :math:`i=0,1,\ldots,n`, compute the coefficients of .. math:: y(t) = c_0 + c_1 \cos(2 \pi t) + c_2 \sin(2 \pi t), so :math:`y(t)` fits the data :math:`(t_i, f_i)` best. A Fourier series approximation is written as .. math:: y(t) = a_0 + \sum_{k=0}^n \left( a_k \cos(2 \pi k t) + b_k \sin(2 \pi k t) \right). The Fast Fourier Transform solves the interpolation problem in time :math:`O(n \log(n))` for :math:`n` data points. There is excellent software available for computing the FFT. Below is a session with Julia which illustrates the transition from the time to the frequency domain. :: using Plots using FFTW dt = 0.01 t = 0:dt:4 # amplitudes 3 and 5 at frequencies 4 and 2 y = 3*sin.(4*2*pi*t) + 5*sin.(2*2*pi*t) plot(t, y, yticks=-7:1:7, label="Amplitude in seconds") F = fft(y) n = length(y)/2 amps = abs.(F)/n freq = [0:79]/(2*n*dt) plot(freq, amps[1:80], xticks=0:1:20, label="Frequency (Hz) versus Amplitude") Observe that the test signal is a linear combination of two sine functions, running at frequencies 4 and 2, and corresponding amplitudes 3 and 5, as can be seen from the definition of :math:`y(t)` below: .. math:: y(t) = 3 \sin(4 \cdot 2 \pi t) + 5 \sin(2 \cdot 2 \pi t). The Julia code makes two plots, shown in :numref:`figShowFFT` (the plot of :math:`y(t)` in the :index:`time domain`), and in :numref:`figShowFFTspectrum`, which shows the essential characterics of the signal in the :index:`frequency domain`. .. _figShowFFT: .. figure:: ./figShowFFT.png :align: center The plot of :math:`y(t) = 3 \sin(4 \cdot 2 \pi t) + 5 \sin(2 \cdot 2 \pi t)`. .. _figShowFFTspectrum: .. figure:: ./figShowFFTspectrum.png :align: center The plot of the Fourier transform of :math:`y(t) = 3 \sin(4 \cdot 2 \pi t) + 5 \sin(2 \cdot 2 \pi t)`. Observe the spikes at 2 and 4 (the frequencies) with corresponding heights 3 and 5 (the amplitudes). The Fast Fourier Transform -------------------------- For :math:`n=4`, consider :math:`\widehat{\bf x} = \frac{1}{\sqrt{4}} M_4 {\bf x}`, with :math:`\omega = e^{-i 2 \pi/4}` and .. math:: M_4 = \left[ \begin{array}{cccc} \omega^0 & \omega^0 & \omega^0 & \omega^0 \\ \omega^0 & \omega^1 & \omega^2 & \omega^3 \\ \omega^0 & \omega^2 & \omega^4 & \omega^6 \\ \omega^0 & \omega^3 & \omega^6 & \omega^9 \end{array} \right] = \left[ \begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & \omega & \omega^2 & \omega^3 \\ 1 & \omega^2 & 1 & \omega^2 \\ 1 & \omega^3 & \omega^2 & \omega \end{array} \right] In :numref:`figComplex4RootsUnity` is the corresponding plot of the four roots of unity. .. _figComplex4RootsUnity: .. figure:: ./figComplex4RootsUnity.png :align: center The four complex roots of unity. Now we introduce the permute and partition operations on the matrix. .. math:: M_4 = \left[ \begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & \omega & \omega^2 & \omega^3 \\ 1 & \omega^2 & 1 & \omega^2 \\ 1 & \omega^3 & \omega^2 & \omega \end{array} \right] \quad P = \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right] Multiplying :math:`M_4` with :math:`P` swaps columns 2 and 3. .. math:: M_4 P = \left[ \begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega & \omega^3 \\ 1 & 1 & \omega^2 & \omega^2 \\ 1 & \omega^2 & \omega^3 & \omega \end{array} \right] = \left[ \begin{array}{cc|cc} 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega & \omega^3 \\ \hline 1 & 1 & \omega^2 & \omega^2 \\ 1 & \omega^2 & \omega^3 & \omega \end{array} \right] Observe the repeated occurrence of :math:`M_2`: .. math:: M_2 = \left[ \begin{array}{cc} 1 & 1 \\ 1 & \omega^2 \end{array} \right], \quad F_2 = \frac{1}{\sqrt{2}} M_2. The matrix :math:`D_2` is obtained via a multiplication with :math:`\omega`. .. math:: M_4 P = \left[ \begin{array}{cc|cc} 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega & \omega^3 \\ \hline 1 & 1 & \omega^2 & \omega^2 \\ 1 & \omega^2 & \omega^3 & \omega \end{array} \right], \quad M_2 = \left[ \begin{array}{cc} 1 & 1 \\ 1 & \omega^2 \end{array} \right], \quad F_2 = \frac{1}{\sqrt{2}} M_2. Multiplying the second row by :math:`\omega` of :math:`M_2`, via a matrix multiplication: .. math:: \left[ \begin{array}{cc} 1 & 0 \\ 0 & \omega \end{array} \right] \left[ \begin{array}{cc} 1 & 1 \\ 1 & \omega^2 \end{array} \right] = \left[ \begin{array}{cc} 1 & 1 \\ \omega & \omega^3 \end{array} \right] = D_2 M_2, where :math:`D_2` is the multiplier matrix: .. math:: D_2 = \left[ \begin{array}{cc} 1 & 0 \\ 0 & \omega \end{array} \right]. Consider flipping signs: .. math:: M_4 P = \left[ \begin{array}{cc|cc} 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega & \omega^3 \\ \hline 1 & 1 & \omega^2 & \omega^2 \\ 1 & \omega^2 & \omega^3 & \omega \end{array} \right], \quad M_2 = \left[ \begin{array}{cc} 1 & 1 \\ 1 & \omega^2 \end{array} \right], \quad D_2 = \left[ \begin{array}{cc} 1 & 0 \\ 0 & \omega \end{array} \right]. and we have .. math:: \left[ \begin{array}{cc} \omega^2 & \omega^2 \\ \omega^3 & \omega \end{array} \right] = \left[ \begin{array}{cc} -1 & -1 \\ -\omega & -\omega^3 \end{array} \right] = - D_2 M_2. The *butterfly matrix* is illustrated in :numref:`figButterfly`. .. _figButterfly: .. figure:: ./figButterfly.png :align: center The 2-dimensional butterfly matrix. .. math:: M_2 = \left[ \begin{array}{cc} 1 & 1 \\ 1 & \omega^2 \end{array} \right] \quad \mbox{and} \quad D_2 = \left[ \begin{array}{cc} 1 & 0 \\ 0 & \omega \end{array} \right] we rewrite :math:`M_4 P` as .. math:: \left[ \begin{array}{cc|cc} 1 & 1 & 1 & 1 \\ 1 & \omega^2 & \omega & \omega^3 \\ \hline 1 & 1 & \omega^2 & \omega^2 \\ 1 & \omega^2 & \omega^3 & \omega \end{array} \right] = \left[ \begin{array}{rr} M_2 & D_2 M_2 \\ M_2 & -D_2 M_2 \end{array} \right] = \left[ \begin{array}{rr} I_2 & D_2 \\ I_2 & -D_2 \end{array} \right] \left[ \begin{array}{cc} M_2 & 0 \\ 0 & M_2 \end{array} \right]. The matrix .. math:: B_4 = \left[ \begin{array}{rr} I_2 & D_2 \\ I_2 & -D_2 \end{array} \right], \quad I_2 = \left[ \begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array} \right], is called the :index:`butterfly matrix`. For :math:`n=4`, the effect of the *butterfly matrix* :math:`B_4` is shown below: .. math:: B_4 \left[ \begin{array}{c} \widehat{x_0} \\ \widehat{x_2} \\ \widehat{x_1} \\ \widehat{x_3} \end{array} \right] = \left[ \begin{array}{rrrr} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & \omega \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -\omega \end{array} \right] \left[ \begin{array}{c} \widehat{x_0} \\ \widehat{x_2} \\ \widehat{x_1} \\ \widehat{x_3} \end{array} \right] = \left[ \begin{array}{l} \widehat{x_0} + \widehat{x_1} \\ \widehat{x_2} + \widehat{x_3} \omega \\ \widehat{x_0} - \widehat{x_1} \\ \widehat{x_2} - \widehat{x_3} \omega \end{array} \right]. For any :math:`n`, we apply the recursive factorization .. math:: M_n = B_n \left[ \begin{array}{cc} M_{n/2} & 0 \\ 0 & M_{n/2} \end{array} \right] P_n where the partition :math:`P_n` reorders the even and odd indexed numbers, with the butterfly matrix .. math:: B_n = \left[ \begin{array}{rr} I_{n/2} & D_{n/2} \\ I_{n/2} & - D_{n/2} \end{array} \right], where :math:`I_{n/2}` is the identity matrix of dimension :math:`n/2`, and :math:`D_{n/2}` is the diagonal matrix with on its diagonal .. math:: 1, \omega, \omega^{2}, \ldots, \omega^{n/2-1}. To the recurrence :math:`T(n) \leq 2 T(n/2) + O(n)`, we apply the following theorem: .. topic:: Theorem on resolving an n*log(n) recurrence. If :math:`T(n) \leq 2 T(n/2) + c_1 n` and :math:`T(2) \leq c_2`, for positive constants :math:`c_1` and :math:`c_2`, then :math:`T(n)` is :math:`O(n \log(n))`. .. topic:: Theorem on the cost of FFT. On vectors of size :math:`n`, the FFT runs in :math:`O(n \log(n))` time. For the inverse FFT, :math:`F^{-1}_n = \overline{F}_n`, so .. math:: F^{-1}_n y = \overline{F}_n y = \overline{F_n \overline{y}}, and no extra work is needed. Filter Design ------------- To design filters, we apply the convolution property of the discrete Fourier transform, shown as a diagram in :numref:`figDFTconvolutionDiagram`. The inverse discrete Fourier transform (iDFT) applied to the componentwise product :math:`\widehat{x} \bullet \widehat{y}` of the discrete Fourier transforms (DFTs) :math:`\widehat{x}$ and $\widehat{y}`, respectively of :math:`x` and :math:`y`, equals the convolution :math:`x \star y`. .. _figDFTconvolutionDiagram: .. figure:: ./figDFTconvolutionDiagram.png :align: center The convolution property of the discrete Fourier transform. As the fast Fourier fourier transform of vectors of size :math:`n` can be computed in time :math:`O(n \log(n))`, the :math:`O(n^2)` cost of a convolution is reduced to time :math:`O(n \log(n))`, as illustrated in :numref:`figFFTconvolutionDiagram`. .. _figFFTconvolutionDiagram: .. figure:: ./figFFTconvolutionDiagram.png :align: center The FFT reduces the cost of a convolution to :math:`O(n \log(n))`. Exercises --------- 1. Use your ``FourierMatrix`` of Exercise 2 of the previous lecture, section 4.1, to compare the outcome of :math:`F_n {\bf x}` computed with the :math:`F_n` computed with your ``FourierMatrix`` and the output of ``FFTW.fft``. * Do the comparison on a random vector :math:`\bf x` of size :math:`n = 8`. * Does the ``FFTW.fft`` use the normalization factor :math:`1/\sqrt{n}`? 2. Write a Julia function ``ButterflyMatrix`` which takes on input the dimension :math:`n` and which returns the butterfly matrix of dimension :math:`n`. To verify the decomposition, define a function ``mFourierMatrix`` that takes on input the dimension :math:`n` and which returns the Fourier matrix times :math:`\sqrt{n}`. For :math:`n=8`, verify the relationship between ``mFourierMatrix`` (*n*) and ``ButterflyMatrix`` (*n*). 3. Generate random vectors for size :math:`2^k`, for :math:`k=10, 11, \ldots, 20`. * For each random vector, record the time for ``FFTW.fft``. * Do you observe the :math:`O(n \log(n))` time? Bibliography ------------ 1. Matteo Frigo and Steven G. Johnson: **The Design and Implementation of FFTW3**. *Proc. IEEE* 93(2): 216--231, 2005. 2. Charles R. MacCluer, Chapter 4 of *Industrial Mathematics. Modeling in Industry, Science, and Government*, Prentice Hall, 2000. 3. Timothy Sauer, chapter 10 of *Numerical Analysis*, second edition, Pearson, 2012.