by Connor Boyle
tags: mathematicssoftware
This post is part 2 in my ongoing series on fast Fourier transform algorithms. You may wish to read part 1 before reading this article.
One of the most useful applications of the fast Fourier transform is in efficiently calculating cyclic convolutions. Paradoxically, cyclic convolutions can also be used to help efficiently calculate discrete Fourier transforms. This post will show the former; I will explain the latter point in future blog posts.
If \(y\) and \(z\) are two sequences of equal length \(N = \lvert y \rvert = \lvert z \rvert\), then the cyclic convolution of those two sequences, \(y \circledast z\), is defined as follows:
\[\left( y \circledast z \right) [j] \triangleq \sum_{m=0}^{N-1} y[m] \cdot z[(j-m) \% N]\] \[| y \circledast z | \triangleq N = |y| = |z|\]Where \(\%\) is the modulo operator.1
The visualization below shows two input sequences \(y\) (vertically, along the left) and \(z\) (horizontally, along the top), the product of each of their elements in the big matrix in the middle. The output convolution, \(y \circledast z\) is shown on the bottom. You can discover the value and origin of an element by hovering the cursor over it. Those with a mouse can click-and-drag the input sequences to change their values.
\(N = \lvert y \rvert = \lvert z \rvert =\)
A cyclic convolution between two sequences \(y\) and \(z\) (length \(\lvert y \rvert = \lvert z \rvert = N\))—calculated naively (see above)—would require \(O(N^2)\) basic operations. However, using the Fourier transform, we can reduce this requirement to \(O(N \log(N))\) operations (assuming \(N\) is highly composite and we calculate the Fourier transforms using the Cooley-Tukey algorithm).
Let’s prove this:
\[\mathcal{F} \{ y \circledast z \}[k] = \sum_{j=0}^{|y \circledast z| - 1} (y \circledast z)[j] \cdot W_{|y \circledast z|}^{jk}\] \[= \sum_{j=0}^{N - 1} \sum_{m=0}^{N-1} y[m] \cdot z[(j-m) \% N] \cdot W_{N}^{jk}\] \[= \sum_{m=0}^{N-1} y[m] \sum_{j=0}^{N - 1} z[(j-m) \% N] \cdot W_{N}^{jk}\]The indices of the inner summation can be split2 into two ranges, from \(0\) to \(m-1\) and from \(m\) to \(N-1\):
\[= \sum_{m=0}^{N-1} y[m] \left( \left( \sum_{j=0}^{m - 1} z[(j-m) \% N] \cdot W_{N}^{jk} \right) + \left( \sum_{j=m}^{N - 1} z[(j-m) \% N] \cdot W_{N}^{jk} \right) \right)\]For the entire outer summation, note that:
\[0 \leq m \leq N-1\] \[1-N \leq -m \leq 0\] \[-N \leq -m \leq 0\]For the first inner summation:
\[\sum_{j=0}^{m - 1} z[(j-m) \% N] \cdot W_{N}^{jk}\] \[0 \leq j \leq m-1\]Since \(-N \leq -m\):3
\[-N \leq j-m \leq -1\] \[(j - m) \% N = j - m + N\]Therefore:
\[\sum_{j=0}^{m - 1} z[(j-m) \% N] \cdot W_{N}^{jk}\] \[= \sum_{j=0}^{m - 1} z[j-m + N] \cdot W_{N}^{jk}\] \[= \sum_{j=N}^{N + m - 1} z[j-m] \cdot W_{N}^{(j-N)k}\] \[= \sum_{j=N}^{N + m - 1} z[j-m] \cdot W_{N}^{jk}\]For the second inner summation:
\[\sum_{j=m}^{N - 1} z[(j-m) \% N] \cdot W_{N}^{jk}\] \[m \leq j \leq N-1\]Since \(-m \leq 0\):3
\[0 \leq j-m \leq N-1\] \[(j - m) \% N = j - m\]Therefore:
\[\sum_{j=m}^{N - 1} z[(j-m) \% N] \cdot W_{N}^{jk}\] \[= \sum_{j=m}^{N - 1} z[j-m] \cdot W_{N}^{jk}\]or, if we define \(\mathcal{F} \{ y \} \odot \mathcal{F} \{ z \}\) as the elementwise product of \(\mathcal{F} \{ y \}\) and \(\mathcal{F} \{ z \}\), i.e.
\[(\mathcal{F} \{ y \} \odot \mathcal{F} \{ z \})[k] \triangleq \mathcal{F} \{ y \}[k] \cdot \mathcal{F} \{ z \}[k]\]then, applying the inverse discrete Fourier transform:
\[(y \circledast z)[k] = \mathcal{F}^{-1} \{ \mathcal{F} \{ y \} \odot \mathcal{F} \{ z \} \}[k]\]Most other proofs of the convolution theorem found on the public-facing web seem to forgo an explicit modulo operator in favor of (sometimes implicitly) treating the input sequences as “\(N\)-periodic” infinitely-repeating signals. For example, Professor Julius Orion Smith III’s free educational book, Mathematics of the Discrete Fourier Transform (DFT) with Audio Applications, does exactly this in its passage on the convolution theorem.
I suspect this convention is intuitive to anyone with sufficient background in signal processing. However, for newcomers such as myself, the unexplained elision of the modulo operator can be rather confusing. It’s also not immediately obvious how to represent and interact with an \(N\)-periodic infinitely-repeating signal in a computer program; the modulo operator, on the other hand, should be quite intuitive to anyone with a computer science background, and much more easily translatable into computer code.
Footnotes:
For our purposes, let the modulo of an integer \(a\) by an integer \(b\) be defined using the remainder from floored division, i.e.
\[a \, \% \, b \triangleq a - b \lfloor \frac{a}{b} \rfloor\]Note that using this definition, the range of the modulo operator is exclusively nonnegative, unlike the modulo operator found in many programming languages such as Javascript and Python. ↩
It is possible for the upper bound of the first sum to be lesser than its lower bound (i.e. \(m = 0 \implies m - 1 < 0\)). Let summations with invalid ranges be defined as 0, i.e.:
\[m = 0 \implies \sum_{j=0}^{m - 1} z[(j-m) \% N] \cdot W_{N}^{jk} = 0\]This proof uses the following property of modular arithmetic; for any \(c \in \mathbb{Z}\):
\[cb \leq a < (c+1)b \implies a \% b = a - cb\]Find me on:
Posts:
Hosted on GitHub Pages — Theme by orderedlist