{"id":479,"date":"2021-05-10T11:30:00","date_gmt":"2021-05-10T11:30:00","guid":{"rendered":"http:\/\/www.ethanepperly.com\/?p=479"},"modified":"2021-05-09T20:52:56","modified_gmt":"2021-05-09T20:52:56","slug":"big-ideas-in-applied-math-the-fast-fourier-transform","status":"publish","type":"post","link":"https:\/\/www.ethanepperly.com\/index.php\/2021\/05\/10\/big-ideas-in-applied-math-the-fast-fourier-transform\/","title":{"rendered":"Big Ideas in Applied Math: The Fast Fourier Transform"},"content":{"rendered":"\n

<\/p>\n\n\n\n

The famous law of the instrument<\/a> states that “when all you have is a hammer, every problem looks like a nail.” In general, this tendency is undesirable: most problems in life are not nails and could better be addressed by a more appropriate tool. However, one can also review the law of the instrument in a more positive framing: when presented with a powerful new tool, it is worth checking how many problems it can solve. The fast Fourier transform (FFT) is one of the most important hammers in an applied mathematician’s toolkit. And it has made many seemingly unrelated problems look like nails.<\/p>\n\n\n\n

In this article, I want to consider three related questions:<\/p>\n\n\n\n

  1. What is the FFT\u2014what problem is it solving and how does it solve it fast?<\/li>
  2. How can the idea<\/em>s behind the FFT be used to solve other problems?<\/li>
  3. How can the FFT be used as a building block in solving a seemingly unrelated problem?<\/li><\/ol>\n\n\n\n

    The FFT is widely considered<\/a> one of the most important numerical algorithms, and as such every sub-community of applied mathematics is inclined to see the most interesting applications of the FFT as those in their particular area. I am unapologetically victim to this tendency myself, and thus will discuss an application of the FFT that I find particularly beautiful and surprising. In particular, this article won’t focus on the manifold applications of the FFT in signal processing, which I think has been far better covered by authors more familiar with that field.<\/p>\n\n\n\n

    The Discrete Fourier Transform<\/h2>\n\n\n\n

    At its core, the FFT is a fast algorithm to compute \"n\" complex numbers \"\hat{f}_0,\ldots,\hat{f}_{n-1}\" given \"n\" real or complex numbers \"f_0,\ldots,f_{n-1}\" defined by the formula1<\/a><\/sup>The factor of \"\frac{1}{\sqrt{n}}\" is not universal. It is common to omit the factor in (1) and replace the \"\tfrac{1}{\sqrt{n}}\" in Eq. (2) with a \"\tfrac{1}{n}\". We prefer this convention as it makes the DFT a unitary<\/a> transformation. When working with Fourier analysis, it is important to choose formulas for the (discrete) Fourier transform and the inverse (discrete) Fourier transform which form a pair in the sense they are inverses of each other.<\/span>\n\n\n

    (1) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    The outputs \"\hat{f} = (\hat{f}_0,\ldots,\hat{f}_{n-1})\" is called the discrete Fourier transform<\/strong> (DFT) of \"f = (f_0,\ldots,f_{n-1})\". The FFT is just one possible algorithm to evaluate the DFT. <\/p>\n\n\n\n

    The DFT has the following interpretation. Suppose that \"f\" is a periodic function defined on the integers with period \"n\"\u2014that is, \"f(j + n) = f(j)\" for every integer \"j\". Choose \"f_0,\ldots,f_{n-1}\" to be the values of \"f\" given by \"f_j = f(j)\" for \"j=0,1,2,\ldots,n-1\". Then, in fact, \"\hat{f}_0,\ldots,\hat{f}_{n-1}\" gives an expression for \"f\" as a so-called trigonometric polynomial2<\/a><\/sup>The name “trigonometric polynomial” is motivated by Euler’s formula which shows that \"e^{(2\pi i/n)jk} = (\cos (\tfrac{2\pi i}{n} j) + i\sin(\tfrac{2\pi i}{n} j))^k\", so indeed the right-hand side of Eq. (2) is indeed a “polynomial” in the “variables” \"\sin(\tfrac{2\pi i}{n} j)\" and \"\cos(\tfrac{2\pi i}{n} j)\"<\/span>:<\/p>\n\n\n

    (2) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    This shows that (1) converts function values \"f_0,f_1,\ldots,f_{n-1}\" of a periodic function \"f\" to coefficients \"\hat{f}_0,\hat{f}_1,\ldots,\hat{f}_{n-1}\" of a trigonometric polynomial representation of \"f\", which can be called the Fourier series<\/strong> of \"f\". Eq. (2), referred to as the inverse discrete Fourier transform<\/strong>, inverts this, converting coefficients \"\hat{f}_0,\hat{f}_1,\ldots,\hat{f}_{n-1}\" to function values \"f_0,f_1,\ldots,f_{n-1}\". <\/p>\n\n\n\n

    Fourier series are an immensely powerful tool in applied mathematics. For example again, if \"f\" represents a sound wave produced by a chord on a piano, its Fourier coefficients \"\hat{f}\" represents the intensity of each pitch comprising the chord. An audio engineer could, for example, compute a Fourier series for a piece of music and zero out Fourier coefficients, thus reducing the amount of data needed to store a piece of music. This idea is indeed part of the way audio compression standards like MP3<\/a> work. In addition to many more related applications in signal processing, the Fourier series is also a natural way to solve differential equations, either by pencil and paper<\/a> or by computer via so-called Fourier spectral methods<\/a>. As these applications (and more to follow) show, the DFT is a very useful computation to perform. The FFT allows us to perform this calculation fast.<\/p>\n\n\n\n

    The Fast Fourier Transform<\/h2>\n\n\n\n

    The first observation to make is that Eq. (1) is a linear transformation<\/a>: if we think of Eq. (1) as describing a transformation \"f \mapsto \hat{f}\", then we have that \"\widehat{\alpha f + \beta g} = \alpha \hat{f} + \beta \hat{g}\". Recall the crucial fact from linear algebra that every linear transformation can be represented by a matrix-vector muliplication<\/a>.3<\/a><\/sup>At least in finite dimensions; the story for infinite-dime\u000ensional vector spaces<\/a> is more complicated.<\/span> In my experience, one of the most effective algorithm design strategies in applied mathematics is, when presented with a linear transformation, to write its matrix down and poke and prod it to see if there are any patterns in the numbers which can be exploited to give a fast algorithm. Let’s try to do this with the DFT. <\/p>\n\n\n\n

    We have that \"\hat{f} = F_n f\" for some \"n\times n\" matrix \"F_n\". (We will omit the subscript \"n\" when its value isn’t particularly important to the discussion.) Let us make the somewhat non-standard choice of describing rows and columns of \"F\" by zero-indexing<\/a>, so that the first row of \"F\" is row \"0\" and the last is row \"n-1\". Then we have that \"\hat{f}_k = \sum_{j=0}^{n-1} F_{kj} f_j\". Comparing with Eq. (1), we see that \"F_{kj} = \tfrac{1}{\sqrt{n}} e^{(-2\pi i/n) jk}\". Let us define \"\omega_n = e^{-2\pi i/n}\". Thus, we can write the matrix \"F\" out as<\/p>\n\n\n

    (3) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    This is a highly structured matrix. The patterns in this matrix are more easily seen for a particular value of \"n\". We shall focus on \"n = 8\" in this discussion, but what will follow will generalize in a straightforward way to \"n\" any power of two (and in less straightforward ways to arbitrary \"n\"\u2014we will return to this point at the end).<\/p>\n\n\n\n

    Instantiating Eq. (3) with \"n = 8\" (and writing \"\omega = \omega_8\"), we have<\/p>\n\n\n

    (4) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    To fully exploit the patterns in this matrix, we note that \"\omega\" represents a clockwise rotation<\/a> of the complex plane by an eighth of the way around the circle. So, for example \"\omega^{21}\" is twenty-one eighths of a turn or simply just \"21-16 = 5\" turns. Thus \"\omega^{21} = \omega^5\" and more generally \"\omega^m = \omega^{m \operatorname{mod} 8}\". This allows us to simplify as follows:<\/p>\n\n\n

    (5) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    Now notice that, since \"\omega\" represents a clockwise rotation of an eighth of the way around the circle, \"\omega^2\" represents a quarter turn of the circle. This fact leads to the surprising observation we can actually find the DFT matrix \"F_4\" for \"n = 4\" hidden inside the DFT matrix \"F_8\" for \"n = 8\"!<\/p>\n\n\n\n

    To see this, rearrange the columns of \"F_8\" to interleave every other column. In matrix language this is represented by right-multiplying with an appropriate4<\/a><\/sup>In fact, this permutation has a special name: the perfect shuffle<\/a>.<\/span> permutation matrix \"\Pi\":<\/p>\n\n\n

    (6) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    The top-left \"4\times 4\" sub-block is precisely \"F_4\" (up to scaling). In fact, defining the diagonal matrix \"\Omega_4 = \operatorname{diag}(1,\omega,\omega^2,\omega^3)\" (called the twiddle factor<\/a>) and noting that \"\omega^4 = -1\", we have<\/p>\n\n\n

    (7) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    The matrix \"F_8\" is entirely built up of simple scalings of the smaller DFT matrix \"F_4\"! This suggests the following decomposition to compute \"F_8x\":<\/p>\n\n\n

    (8) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    Here \"x_1\" represent the even-indexed entries of \"x\" and \"x_2\" the odd-indexed entries. Thus, we see that we can evaluate \"F_8 x\" by evaluating the two expressions \"F_4x_1\" and \"F_4x_2\". We have broken our problem into two smaller problems, which we then recombine into a solution of our original problem.<\/p>\n\n\n\n

    How then, do we compute the smaller DFTs \"F_4x_1\" and \"F_4x_2\"? We just use the same trick again, breaking, for example, the product \"F_4x_1\" into further subcomputations \"F_2x_{11}\" and \"F_2x_{12}\". Performing this process one more time, we need to evaluate expressions of the form \"F_1x_{111}\", which are simply given by \"F_1 x_{111}= x_{111}\" since the matrix \"F_1\" is just a \"1\times 1\" matrix whose single entry is \"1\".<\/p>\n\n\n\n

    This procedure is an example of a recursive algorithm<\/a>: we designed an algorithm which solves a problem by breaking it down into one or more smaller problems, solve each of the smaller problems by using this same algorithm, and then reassemble the solutions of the smaller problems to solve our original problem. Eventually, we will break our problems into such small pieces that they can be solved directly, which is referred to as the base case<\/a> of our recursion. (In our case, the base case is multiplication by \"F_1\"). Algorithms using this recursion in this way are referred to as divide-and-conquer algorithms<\/a>.<\/p>\n\n\n\n

    Let us summarize this recursive procedure we’ve developed. We want to compute the DFT \"y = F_n x\" where \"n\" is a power of two. First, we use the DFT to recursively compute \"F_{n/2}x_1\" and \"F_{n/2}x_2\". Next, we combine these computations to evaluate \"y = F_n x\" by the formula<\/p>\n\n\n

    (9) <\/span>   <\/span>\"\begin{equation*}<\/p><\/p>\n\n\n

    This procedure is the famous fast Fourier transform (FFT), whose modern incarnation was presented by Cooley and Tukey in 1965 with lineage that can be traced back to work by Gauss in the early 1800s. There are many variants of the FFT using similar ideas.<\/p>\n\n\n\n

    Let us see why the FFT is considered “fast” by analyzing its operation count. As is common for divide-and-conquer algorithms, the number of operations for computing \"F_n x\" using the FFT can be determined by solving a certain recurrence relation<\/a>. Let \"T(n)\" be the number of operations required by the FFT. Then the cost of computing \"F_n x\" consists of<\/p>\n\n\n\n