jsEval "curl_module('../haskell/Diagram.ob')"
jsEval "curl_module('Matrix.ob')"
DFT Speedrun
Texts on Fourier analysis typically begin by explaining why you’d want to decompose a given signal into sinusoid waves. However, this is no gentle introduction to the subject. Rather, this is a speedrun through the discrete Fourier Transform with linear algebra and a dash of functional programming. For a saner approach, try Julius O. Smith III, Mathematics of the Discrete Fourier Transform (DFT).
We view the discrete Fourier transform as a change of basis in \(\mathbb{C}^N ,\) one that rewrites standard coordinates in terms of an orthogonal basis built using the \(N\)th roots of unity:
\[ \newcommand{\u}{\mathbf{u}} \newcommand{\v}{\mathbf{v}} \newcommand{\x}{\mathbf{x}} \newcommand{\y}{\mathbf{y}} \newcommand{\X}{\mathbf{X}} \DeclareMathOperator{\DFT}{dft} \DeclareMathOperator{\BAR}{bar} \DeclareMathOperator{\TAILWHIP}{tailwhip} \DeclareMathOperator{\ROTATE}{rotate} \DeclareMathOperator{\STRETCH}{stretch} \DeclareMathOperator{\REPEAT}{repeat} \DeclareMathOperator{\SELECT}{select} \DeclareMathOperator{\ALIAS}{alias} z^0, z^1, …, z^{N-1} \]
where \( z = e^{2\pi i / N} ,\) the (primitive) \(N\)th root of unity with the smallest positive angle.
We draw the cases \(N = 2,…,8\) on the complex plane, with lines connecting one root of unity to the next so they appear as regular polygons. Observe each the \(N\)-gon looks the same if rotated by \(1/N\) of a turn, or if it is flipped upside-down. We exploit these symmetries later.
-- The Matrix module is unaware of complex numbers, so its dot product -- fails to conjugate the second argument. import Matrix hiding (scale, dot) -- The Diagram module treats complex numbers as 2D vectors, so its -- dot product means something else. import Diagram hiding (dot) dot u v = sum $ zipWith (*) u $ conjugate <$> v ngons = hcat $ intersperse (strutX 1) $ scale 1.5 . cyclogon <$> [2..8] jsEval_ $ "ngonsDiv.innerHTML = `" ++ svg 24 ngons ++ "`;"
For an integer \(m\), define \( \v_m \) to be the vector consisting of each of the roots raised to the power of \(m:\)
\[ \v_m = (z^0, z^m, ..., z^{(N-1)m}) \]
Since \(z^N = 1,\) we can liken this to a card game with \(N\) people in a circle where the turn order follows every \(m\)th person (a far more humane version of the Josephus problem). Some examples in a 10-player game:
skip k xs = take n $ (xs!!) . (`mod` n) <$> [0, k..] where n = length xs skip 1 [0..9] skip -2 [0..9] skip (10-2) [0..9] skip 3 [0..9] skip 5 [0..9] skip 10 [0..9]
Consider the dot product of two such vectors:
\[ \v_a \cdot \v_b = \sum_{k=0}^{N-1} z^{(a - b)k} \]
When \(a = b \bmod N\) this simplifies to \(N\), implying \(|\v_a|^2 = N.\) Otherwise the geometric series sums to:
\[ \frac{1 - z^{N(a - b)}}{1 - z^{a - b}} = 0 \]
For example:
cast = fromIntegral zs n = take n $ iterate (cis(tau/cast n) *) 1 -- Roots of unity. trunc (a :+ b) = go a :+ go b where go x = if abs x < 0.000001 then 0 else x trunc $ skip 3 (zs 7) `dot` skip 3 (zs 7) -- N = 7, v_3 . v_3 trunc $ skip 3 (zs 7) `dot` skip 4 (zs 7) -- N = 7, v_3 . v_4 trunc $ skip 1 (zs 8) `dot` skip 2 (zs 8) -- N = 8, v_1 . v_2 trunc $ skip 4 (zs 8) `dot` skip 4 (zs 8) -- N = 8, v_4 . v_4
Thus \( \{\v_0, …, \v_{N-1}\} \) is an orthogonal set, and in fact must be an orthogonal basis of \(\mathbb{C}^N\) because it has \(N\) members. We could normalize each vector to produce an orthonormal basis \( \{\hat{\v}_0, …, \hat{\v}_{N-1}\} \) and then define a nice clean DFT by the map taking a given vector \(\x\) to the vector whose \(k\)th component is simply \(\x\cdot\hat{\v}_k.\)
However, in practice, we prefer to avoid division by \(\sqrt{N},\) so we leave the basis unnormalized. Define the discrete Fourier transform (DFT) of \(\x\) to be the vector \(\DFT(\x)\) whose \(k\)th component is:
\[ \begin{align} \DFT_k(\x) &= \x \cdot \v_k \\ &= \sum_{n=0}^{N-1} x_n \overline{z^{nk}} \\ &= \sum_{n=0}^{N-1} x_n e^{-2\pi i k n / N} \end{align} \]
Let \(\X = (X_0, …, X_{N-1}) = \DFT(\x). \) The inverse discrete Fourier transform (IDFT) is:
\[ \x = \frac{1}{N} \sum_{k=0}^{N-1} X_k \v_k \]
In particular, the \(n\)th component is:
\[ \begin{align} x_n &= \frac{1}{N} \sum_{k=0}^{N-1} X_k z_{nk} \\ &= \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{2\pi i k n / N} \end{align} \]
Since the DFT is essentially a change of basis, and an orthogonal one at that, it is a linear operator that can be represented by a matrix \(W\) so that \(\X = W \x:\)
\[ W = \begin{bmatrix} z^0 & z^0 & ... & z^0 \\ z^0 & z^{-1} & ... & z^{-(N-1)} \\ \vdots & \vdots & \ddots & \vdots \\ z^0 & z^{-(N-1)} & ... & z^{-(N-1)^2} \\ \end{bmatrix} \]
The inverse is:
\[ W^{-1} = \frac{1}{N} \begin{bmatrix} z^0 & z^0 & ... & z^0 \\ z^0 & z^{1} & ... & z^{N-1} \\ \vdots & \vdots & \ddots & \vdots \\ z^0 & z^{N-1} & ... & z^{(N-1)^2} \\ \end{bmatrix} = N^{-1} \overline{W} \]
The following code shows the DFT and IDFT matrices for \(N = 3.\) We can check by eye that \(W^{-1} = N^{-1} \overline{W}.\)
dftMatrix n = [skip k $ conjugate <$> zs n | k <- [0..n-1]] undftMatrix n = [skip k $ (/cast n) <$> zs n | k <- [0..n-1]] dft v = matrixFun (dftMatrix $ length v) v undft v = matrixFun (undftMatrix $ length v) v dftMatrix 3 undftMatrix 3 map (map ((/3) . conjugate)) $ dftMatrix 3
We’re well over the finish line, but let’s go for a higher completion percentage. We’ve already unlocked two achievements: the above imply the DFT editions of the power theorem:
\[ \u \cdot \v = \frac{1}{N} \DFT(\u) \cdot \DFT(\v) \]
\[ |\v|^2 = \frac{1}{N} |\DFT(\v)|^2 \]
Reverse gear
Some famous properties of the DFT arise from symmetries of the roots of unity, which we express as permutations.
The DFT is a bunch of dot products, so we manufacture identities by noticing that the dot product of two vectors remains unchanged if we permute the coordinates of each vector in the same way. This is because we wind up multiplying the same pairs of coordinates together before taking their sum. That is, if \(P\) is a permutation then:
\[ \x \cdot \v = P\x \cdot P\v \]
Replacing \(\v\) with \(P^{-1} \v\) gives:
\[ \x \cdot P^{-1} \v = P\x \cdot \v \]
For example, instead of iterating the roots of unity in counterclockwise order, we can go clockwise:
\[ \v_{-1} = (z^0, z^{-1}, …, z^{-(N-1)}) = (z^0, z^{N-1}, …, z^1) \]
More generally:
\[ \v_{-k} = (z^0, z^{(N-1)k}, …, z^k) \]
We see the roots of unity are permuted: the first coordinate stays in place, while the others are reversed. Define \(\TAILWHIP : \mathbb{C}^N → \mathbb{C}^N\) to be this permutation:
\[ \TAILWHIP : (x_0, x_1, …, x_{N-1}) \mapsto (x_0, x_{N-1}, …, x_1) \]
tailwhip (h:t) = h:reverse t tailwhip "tailwhip" tailwhip <$> ["banana", "flatcar", "hydra", "penal", "retinue", "tabu"]
From \(z^N = 1,\) we have \(\v_{-k} = \v_{N-k},\) and hence:
\[ \v_{N-k} = \TAILWHIP(\v_k) \]
skip -2 $ zs 5 tailwhip $ skip 2 $ zs 5
Since \(\TAILWHIP = \TAILWHIP^{-1}:\)
\[ \begin{align} \DFT_{N-k} (\x) &= \x \cdot \v_{N-k} \\ &= \x \cdot \TAILWHIP(\v) \\ &= \TAILWHIP(\x) \cdot \v_k \\ &= \DFT_k(\TAILWHIP(\x)) \end{align} \]
which implies:
\[ \TAILWHIP \circ \DFT = \DFT \circ \TAILWHIP \]
tailwhip $ dft [1, 2, 3 :+ 4] dft $ tailwhip [1, 2, 3 :+ 4]
DFT Squared
I prefer the following proof, because functional programmers naturally wonder what happens if you repeatedly apply the same function. For linear operators in particular, this line of inquiry often bears fruit.
Consider the entry of the matrix \(W^2\) at row \(r\) and column \(c:\)
\[ \sum_{k=0}^{N-1} z^{-kr} z^{-kc} = \sum_{k=0}^{N-1} z^{-k(r+c)} \]
Earlier we showed this is \(N\) when \(r + c = 0 \bmod N\) and zero otherwise. Thus \(W^2 = N T\) where:
\[ T = \begin{bmatrix} 1 & 0 & … & 0 \\ 0 & 0 & & 1 \\ \vdots & & \unicode{x22F0} & \\ 0 & 1 & & 0 \\ \end{bmatrix} \]
which describes the \(\TAILWHIP\) permutation. Then \( W^2 W = W W^2 \) implies:
\[ T W = W T \]
which is a succinct form of \( \TAILWHIP \circ \DFT = \DFT \circ \TAILWHIP.\) QED.
We also have \(W^4 = N^2, \) or alternatively \( W^{-1} = N^{-2} W^3.\)
Above, we found \( W^{-1} = N^{-1}\overline{W} \) thus \(\overline{W} = N^{-1}W^3.\)
Let’s try out the \(N = 5\) case:
m5 = Matrix $ dftMatrix 5 m5 trunc <$> m5^2 trunc <$> m5^4
Upside-down
The vertical symmetry of the roots of unity can be expressed via complex conjugation, which we denote with an overbar. However, we also want a first-class name for complex conjugation, so we define \(\BAR : \mathbb{C}^N → \mathbb{C}^N \) by \(\x \mapsto \overline{\x}: \)
\[ (x_0, …, x_{N-1}) \mapsto (\overline{x_0}, …, \overline{x_{N-1}}) \]
Flipping the roots of unity vertically leads to a familiar permutation:
\[ \begin{align} \BAR(\v_k) &= (\overline{z^0}, \overline{z^k}, …, \overline{z^{(N-1)k}} ) \\ &= (z^0, z^{-k}, …, z^{-(N-1)k} ) \\ &= \v_{N-k} \end{align} \]
The fact that complex conjugation distributes over the dot product looks clunky with our notation: \(\BAR(\u \cdot \v) = \BAR(\u) \cdot \BAR(\v).\) We use this identity below:
\[ \begin{align} \DFT_k(\BAR(\x)) &= \BAR(\x) \cdot \v_k \\ &= \BAR(\x) \cdot \BAR(\BAR(\v_k)) \\ &= \BAR( \x \cdot \BAR(\v_k) ) \\ &= \BAR (\x \cdot \v_{N-k} ) \\ &= \BAR (\DFT_{N-k}(\x)) \end{align} \]
Therefore:
\[ \DFT \circ \BAR = \BAR \circ \TAILWHIP \circ \DFT = \TAILWHIP \circ \BAR \circ \DFT \]
where we have used the fact that \(\BAR\) commutes with any permutation.
As \(\TAILWHIP = \TAILWHIP^{-1}, \) left-composing both sides by this function yields:
\[ \TAILWHIP \circ \DFT \circ \BAR = \BAR \circ \DFT \]
When each coordinate of \(\x\) is real:
\[ \begin{align} & \BAR(\x) = \x \\ & \implies \DFT \circ \TAILWHIP = \TAILWHIP \circ \DFT = \BAR \circ \DFT \end{align} \]
bar = map conjugate tailwhip $ dft [3,1,4,1,5] bar $ dft [3,1,4,1,5]
Let \(T\) be the matrix representing the tailwhip permutation. If \(T \x = \x\) then we say \(\x\) is even and if \(T \x = -\x\) then we say \(\x\) is odd.
Thus the previous identity means that if \(\x\) is real, then \(\Re(\X)\) is even and \(\Im(\X)\) is odd where \(X = \DFT(\x).\) Similarly, replacing each component \(X_k\) of \(\X\) with its magnitude \(|X_k|\) results in an even vector, and replacing each of them with with its angle \(\angle X_k\) results in an odd vector.
The identity \(T W = W T\) implies if \(\x\) is even then \(\DFT(\x)\) is even.
If \(\x\) is real and even, that is \(T\x = \x = \overline{\x},\) then:
\[ \overline{W\x} = \overline{W} \overline{\x} = \overline{W} \x = N^{-1} W^3 \x = W T \x = W \x \]
or in other words, \(\DFT(\x)\) is also real and even.
Going around in circles
We now exploit the rotational symmetry of the roots of unity. For any integer \(d,\) define the permutation \(\ROTATE_d\) by:
\[ \ROTATE_d : (x_0, …, x_{N-1}) \mapsto (x_{-d}, …, x_{N-1-d}) \]
where subscripts are reduced modulo \(N\) to one of \([0..N-1].\)
rotate d xs = drop k xs ++ take k xs where k = -d `mod` length xs rotate 3 <$> ["alloy", "kingpin", "shotgun"] rotate -2 <$> ["ads", "elbow", "enlist", "stripe"] rotate 3 <$> [[0..9], [0..4]] rotate -2 <$> [[0..9], [0..4]]
Then:
\[ \begin{align} \ROTATE_d(\v_k) &= \ROTATE_d(z^0, z^k, …,z^{(N-1)k}) \\ &= (z^{-kd}, …, z^{-k}, z^0, z^k, …,z^{(N-1-d)k}) \\ &= z^{-kd} \v_k \end{align} \]
Taking the dot product with \(\x\) on both sides proves the shift theorem:
\[ \DFT_k \circ \ROTATE_d = (z^{-kd} \cdot) \circ \DFT_k \]
where we’ve borrowed the operator section syntax of Haskell so that \((a \cdot)\) denotes the function that multiplies by \(a\).
This is a good time to bust out the Hadamard product (\(\odot\)), which multiplies together corresponding elements of its inputs:
\[ \DFT \circ \ROTATE_d = (\v_{-d} \odot) \circ \DFT \]
hadamard = zipWith (*) trunc <$> dft (rotate 2 [2, 7, 1 :+ 8, 2 :+ 8]) trunc <$> skip (-2) (zs 4) `hadamard` dft [2, 7, 1 :+ 8, 2 :+ 8]
Productive Results
We turn our attention to operations that involve products of coordinates.
The convolution of \(\x\) and \(\y\) is given by:
\[ (\x * \y)_n = \sum_{m=0}^{N-1} x_m y_{n - m} = \x \cdot (\ROTATE_n \circ \TAILWHIP \circ \BAR)(\y) \]
where subscripts are reduced modulo \(N\) to one of \([0..N-1].\)
The Convolution Theorem:
\[ \DFT_k(\x * \y) = \DFT_k(\x) \DFT_k(\y) \]
Again we can use the Hadamard product to remove the \(k\):
\[ \DFT(\x * \y) = \DFT(\x) \odot \DFT(\y) \]
Proof: We rely on interchanging the order of summations, so we expand the definitions and slog through:
\[ \begin{align} \DFT_k(\x * \y) &= \sum_{n=0}^{N-1} (\x * \y)_n z^{-kn} \\ &= \sum_{n=0}^{N-1} \sum_{m=0}^{N-1} x_m y_{n - m} z^{-kn} \\ &= \sum_{m=0}^{N-1} x_m \sum_{n=0}^{N-1} y_{n - m} z^{-kn} \\ &= \sum_{m=0}^{N-1} x_m \DFT_k (\ROTATE_m (\y)) \\ &= \sum_{m=0}^{N-1} x_m z^{-km} \DFT_k (\y) \\ &= \DFT_k(\x) \DFT_k (\y) \end{align} \]
where we apply the shift theorem in one of the steps. QED.
convolve x = map (dot x) . take (length x) . iterate (rotate 1) . tailwhip . bar dft $ [1,2,3] `convolve` [4,5,6] dft [1,2,3] `hadamard` dft [4,5,6]
Dual of the convolution theorem:
\[ \DFT(\x \odot \y) = \frac{1}{N} \DFT(\x) * \DFT(\y) \]
Proof: With matrices, the convolution theorem reads:
\[ W(\x * \y) = W(\x) \odot W(\y) \]
Replace \(\x, \y\) by \(W^{-1}\x, W^{-1}\y :\)
\[ W(W^{-1}\x * W^{-1}\y) = \x \odot \y \]
Apply \(W\) to both sides, and write the equation the other way:
\[ W(\x \odot \y) = W^2(W^{-1}\x * W^{-1}\y) \]
Above we showed \(W^2 = N T,\) where \(T\) is the matrix representing the \(\TAILWHIP\) permutation, hence:
\[ W(\x \odot \y) = NT(W^{-1}\x * W^{-1}\y) \]
Recall we argued that a dot product remains unchanged if we permute the coordinates of both its arguments in the same way. In a similar fashion, since \(T\) swaps \(k\) with \(N-k,\) from the definition of convolution, we have:
\[ T(\x * \y) = T\x * T\y \]
Therefore:
\[ \begin{align} W(\x \odot \y) &= NT(W^{-1}\x * W^{-1}\y) \\ &= N(W^{-1}T\x * W^{-1}T\y) \\ &= N(W^{-1}N^{-1}W^2 \x * W^{-1}N^{-1}W^2 \y) \\ &= N^{-1}(W\x * W\y) \end{align} \]
QED. It’s times like these we’d prefer a tidier DFT defined with an orthonormal basis, for then there would be no pesky \(N\) factors polluting our expressions.
trunc <$> dft ([1,2,3] `hadamard` [4,5,6]) trunc <$> (/3) <$> dft [1,2,3] `convolve` dft [4,5,6]
The correlation of two vectors \(\x\) and \(\y\) is:
\[ (\x \star \y)_n = \sum_{m=0}^{N-1} \overline{x_m} y_{n + m} = \ROTATE_{-n}(\y) \cdot \x \]
The correlation theorem:
\[ \DFT (\x\star\y) = \overline{\DFT(\x)} \odot \DFT(\y) \]
As usual, subscripts are reduced modulo \(N\) to one of \([0..N-1].\)
Proof:
\[ \begin{align} (\x \star \y)_n &= \sum_{m=0}^{N-1} \overline{x_m} y_{n + m} \\ &= \sum_{m'=0}^{N-1} \overline{x_{N-m'}} y_{n-m'} \\ &= ((\TAILWHIP \circ \BAR)(\x) * \y)_n \\ \end{align} \]
Hence:
\[ \DFT(\x \star \y) = \DFT((\TAILWHIP \circ \BAR)(\x) * \y) \]
The theorem follows from the convolution theorem and the identity involving \(\BAR\) we found above. QED.
correlate x y = map (`dot` x) $ take (length x) $ iterate (rotate -1) y trunc <$> dft ( [1, 2, 3] `correlate` [7:+ -6, 8:+ -5, 9:+ -4]) trunc <$> bar (dft [1, 2, 3]) `hadamard` dft [7:+ -6, 8:+ -5, 9:+ -4]
Stretch Goals
Other results relate to functions that change the length of inputs.
The \(\STRETCH_L\) function increases the length of a vector by factor of \(L\) by inserting \(L - 1\) zeroes after each element.
stretch ell = concatMap (: replicate (ell-1) 0)
The \(\REPEAT_L\) function repeats the elements of the input vector \(L\) times to produce a vector that is \(L\) times as long.
repeat' ell = concat . replicate ell
The stretch theorem, or repeat theorem:
\[ \DFT \circ \STRETCH_L = \REPEAT_L \circ \DFT \]
By abuse of notation \(\DFT\) refers to two different maps. They are both discrete Fourier transforms, but act on different vector spaces. The left one is \(\mathbb{C}^{NL} → \mathbb{C}^{NL}\) and the right one is \(\mathbb{C}^N → \mathbb{C}^N.\)
Proof: Let \(z' = z^{1/L}.\) To disambiguate, write \(\DFT'\) for the DFT on \(\mathbb{C}^{NL}\) and \(\DFT\) for the DFT on \(\mathbb{C}^N.\) For \(k \in [0..N-1]\) we have:
\[ \begin{align} \DFT'_k(\STRETCH_L(\x)) &= \STRETCH_L(\x) \cdot (z'^0, z'^k, ..., z'^{(N-1)k}) \\ &= (x_0, 0..0, x_1, 0..0, ..., x_{N-1}, 0..0) \cdot (z'^0, z'^k, ..., z'^{(N-1)k}) \end{align} \]
where \(0..0\) denotes a run of \(L - 1\) zeroes. Then only every \(L\)th coordinate can contribute nontrivially to the dot product, hence:
\[ \begin{align} \DFT'_k(\STRETCH_L(\x)) &= (x_0, x_1, ..., x_{N-1}) \cdot (z'^0, z'^{Lk}, ..., z'^{L(N-1)k}) \\ &= (x_0, x_1, ..., x_{N-1}) \cdot (z^0, z^k, ..., z^{(N-1)k}) \\ &= \DFT_k(\x) \end{align} \]
where we have used \(z'^L = z.\)
Since \(z^N = 1\), if \(k = qN + r\) for some integer \(q\) and \(r \in [0..N-1],\) retracing the above steps leads to:
\[ \DFT'_k(\STRETCH_L(\x)) = \DFT_r(\x) \]
thus taking \(k\) over the entire range of \([0..NL-1]\) shows:
\[ \DFT'(\STRETCH_L(\x)) = \REPEAT_L(\DFT(\x)) \]
QED.
The \(\SELECT_L\) function picks out every \(L\)th coordinate of a vector, whose length which must be a multiple of \(L:\)
chunksOf i ls = map (take i) (go ls) where -- See Data.List.Split. go [] = [] go l = l : go (drop i l) select l = map head . chunksOf l select 3 <$> ["benevolently", "speedways", "theatrically"]
The \(\ALIAS_L\) function takes a vector whose length must be a multiple of \(L,\) splits it into \(L\) vectors of the same length and adds them all together:
alias l v = foldl1 (zipWith (+)) $ chunksOf (length v `div` l) v alias 2 [1..10] alias 5 [1..10]
The downsampling theorem, or alias theorem:
\[ \DFT \circ \SELECT_L = (/L) \circ \ALIAS_L \circ \DFT \]
Proof:
Let \(\x\) be a vector of length \(NL\) and as usual let \(z\) be the \(N\)th root of unity with the smallest positive angle, and let \(\DFT\) be the DFT for vectors of size \(N.\) (Even though \(\x\) has length \(LN\) we still want \(\DFT\) to refer to the length-\(N\) version of the DFT.)
Let \(z' = z^{1/L}\) and define:
\[\v'_k = (z'^0, z'^k, ..., z'^{(LN-1)k})\]
Again, to disambiguate, let \(\DFT'\) be the DFT for vectors of size \(LN,\) that is:
\[ \DFT'_k (\x) = \x \cdot \v'_k \]
As with earlier proofs, we rely on various geometric series summing to something simple, often zero. However, there are many of them, and we just have to wade through the details. We start with:
\[ \begin{align} \ALIAS_L (\DFT' (\x)) &= \ALIAS_L((\x \cdot \v'_0, \x \cdot \v'_1, ..., \x \cdot \v'_{N-1})) \\ \end{align} \]
Since the dot product is bilinear, the \(k\)th component of this vector is:
\[ \x \cdot \v'_k + \x \cdot \v'_{N+k} + ... + \x \cdot \v'_{(L-1)N + k} = \x \cdot (\v'_k + \v'_{N+k} + ... + \v'_{(L-1)N + k}) \]
We write down the coefficients of the vectors in the sum:
\[ \begin{matrix} \v'_k &= & z'^0 & z'^k & z'^{2k} & ... & z'^{(LN-1)k} \\ \v'_{N+k} &= &z'^0 & z'^{N+k} & z'^{2(N+k)} & ... & z'^{(LN-1)(N+k)} \\ &\vdots \\ \v'_{(L-1)N + k} &= & z'^0 & z'^{(L-1)N+k} & z'^{2((L-1)N+k)} & ... & z'^{(LN-1)((L-1)N+k)} \\ \end{matrix} \]
Then the \(c\)th coefficients sum to:
\[ z'^{ck} + z'^{c(N+k)} + ... + z'^{c((L-1)N+k)} = z'^{ck} (z'^0 + z'^{cN} + ... + z'^{(L-1)cN}) \]
If \(c = mL\) for some integer \(m\), each term of the sum is a power of \(z'^{LN} = 1\) so the entire expression reduces to:
\[ z'^{mLk}(1 + ... + 1) = L (z'^L)^{mk} = L z^{mk} \]
Otherwise \(z'^{cN} \ne 1\) so the geometric series sums to:
\[ \frac{1 - z'^{LcN}}{1- z'^{cN}} = 0 \]
Thus:
\[ \v'_k + \v'_{N+k} + ... + \v'_{(L-1)N + k} = L(z^0, 0..0, z^k, 0..0, ..., z^{(L-1)k}, 0..0) \]
where each run of zeroes is of length \(L-1.\) Then the above dot product becomes:
\[ \begin{align} \x \cdot L(z^0, 0..0, z^k, 0..0, ..., z^{(L-1)k}, 0..0) &= \SELECT_L(\x) \cdot L(z^0, z^k, ..., z^{(L-1)k}) \\ &= L \SELECT_L(\x) \cdot \v_k \\ &= L \DFT_k (\SELECT_L(\x)) \end{align} \]
QED.
map (3*) $ dft $ select 3 $ zipWith (:+) [11..19] [21..29] alias 3 $ dft $ zipWith (:+) [11..19] [21..29]
A permutation by any other name…
We’ve been stubbornly viewing the DFT as a function taking a vector of \(\mathbb{C}^N\) when the input really should be seen as an infinite periodic signal \( …, x[-1], x[0], x[1], … \) that repeats every \(N\) samples, namely \(x[n] = x[m]\) whenever \(n = m \bmod N.\)
Thus instead of "tailwhip", normal people say "flip", because the tailwhip permutation we defined can simply be written \(x[-n].\) This also explains the terms "even" and "odd".
Similarly, instead of "rotate", they say "shift", because the rotation by \(d\) we defined can simply be written \(x[n - d].\)
Fast Fourier Transform
Let \(N = N_1 N_2.\) Write \(n = N_1 n_2 + n_1\) and \(k = N_2 k_1 + k_2.\)
Then:
\[\begin{align}
X_k &= \sum_{n_1=0}^{N_1-1} \sum_{n_2=0}^{N_2-1} x_{N_1 n_2 + n_1}
e^{-\frac{2\pi i}{N_1 N_2} (N_1 n_2 + n_1)(N_2 k_1 + k_2) } \\
&= \sum_{n_1=0}^{N_1-1}
\left[ e^{-\frac{2\pi i}{N}n_1 k_2} \right]
\left(
\sum_{n_2=0}^{N_2-1} x_{N_1 n_2 + n_1}
e^{-\frac{2\pi i}{N_2} n_2 k_2}
\right)
e^{-\frac{2\pi i}{N_1} n_1 k_1}
\end{align}\]
and we have reduced the original DFT to \(N_1 + N_2\) smaller DFTs; doing so recursively yields the fast Fourier transform (FFT). The expression in square brackets is known as a "twiddle factor", and is an \(N\)th root of unity that can be precomputed.