diff options
Diffstat (limited to 'src/math/fft.org')
| -rw-r--r-- | src/math/fft.org | 141 |
1 files changed, 141 insertions, 0 deletions
diff --git a/src/math/fft.org b/src/math/fft.org new file mode 100644 index 0000000..45e9595 --- /dev/null +++ b/src/math/fft.org @@ -0,0 +1,141 @@ +#+title: Fast Fourier Transform +#+subtitle: Implementation in Common Lisp +#+author: Thomas Albers Raviola +#+date: <2023-03-06> +#+setupfile: ../../math_options.org + +In this article I attept to explain the principles of the fast fourier transform +(FFT) and develop in incremental steps an implementation in common lisp. + +*NOTE*: This is part of my personal notes and as such may contain mistakes or +not be correct at all. If you find any issues please contact me at my email +address. + +* Introduction + +The FFT is an algorithm for computing the discrete fourier transform (DFT) of an +array of values that is meant to be faster than the usual naive implementation. +Instead of having the usual \[ O(n^2) \] behaviour it is an \[ O(n\log(n)) \] +algorithm. + +* Derivation + +To derive the FFT we first consider the definition of the DFT + +Let \[ \xi_k \] be an array of N, possible complex, numbers. The DFT of \[ \xi_k +\] is an array of the same length \[ x_k \] defined as: + +{{{beg-eqn}}} +x_k = \sum_{n=0}^{N-1} \xi_n e^{- \frac{2 \pi i}{N} n k} +{{{end-eqn}}} + +A possible implementation of this in code would be something like the following: + +#+begin_src common-lisp + (defun dft (x) + "Naive discrete fourier transform, not fft" + (flet ((transform (k) + (loop :with i = (complex 0 1) + :for n :below (length x) + :sum (* (aref x n) (exp (- (/ (* 2 pi i k n) (length x)))))))) + (map-seq-iota #'transform (length x)))) +#+end_src + +But we can do a lot better in terms of performance. To this end we first split +the sum {{{ref}}} into even and odd members + +{{{beg-align}}} +x_k &= \sum_{n=0}^{N/2-1} \xi_{2n} e^{- \frac{2 \pi i}{N} (2n) k} + +\sum_{n=0}^{N/2-1} \xi_{2n+1} e^{- \frac{2 \pi i}{N} (2n+1) k}\\ +&= \sum_{n=0}^{N/2-1} \xi_{2n} e^{- \frac{2 \pi i}{N/2} n k} + +e^{-\frac{2 \pi i}{N}k} \sum_{n=0}^{N/2-1} \xi_{2n+1} e^{- \frac{2 \pi i}{N/2} n +k}\\ +&=E_k + e^{-\frac{2 \pi i}{N}k} O_k +{{{end-align}}} + +Where \[ E_k \] and \[ O_k \] themselves are FFTs of a smaller array each. +Furthermore notice how thanks to the cyclic nature of the complex exponential +function it is possible to half the number of computations: + +{{{beg-align}}} +x_{k+\frac{N}{2}} &= E_{k+\frac{N}{2}} + e^{-\frac{2 \pi i}{N} (k + +\frac{N}{2})} O_{k+\frac{N}{2}}\\ +&= E_k - e^{-\frac{2 \pi i}{N} k} O_k +{{{end-align}}} + +Thus + +{{{beg-eqn}}} +\begin{aligned} +x_k &= E_k + e^{-\frac{2 \pi i}{N} k} O_k\\ +x_{k+\frac{N}{2}} &= E_k - e^{-\frac{2 \pi i}{N} k} O_k +\end{aligned}\,\text{,}\quad k=0,\dots,\frac{N}{2}-1 +{{{end-eqn}}} + +* Implementing the FFT + +With the reduction identities we derived for the FFT a preliminary +implementation could be: + +#+begin_src common-lisp + (defun fft (sequence) + (labels + ((fft-step (n step offset) + (let ((ret (make-array n))) + (if (= n 1) + (setf (aref ret 0) (aref sequence offset)) + (loop + :with n2 = (/ n 2) + :with ek = (fft-step n2 (* 2 step) offset) + :with ok = (fft-step n2 (* 2 step) (+ offset step)) + :for k :below n2 + :for c = (exp (- (/ (* 2 pi (complex 0 1) k) n))) + :do + (setf (aref ret k) (+ (aref ek k) (* c (aref ok k))) + (aref ret (+ k n2)) (- (aref ek k) (* c (aref ok k)))))) + ret))) + (fft-step (length sequence) 1 0))) +#+end_src + +The function ~fft-step~ needs to be able to access both even and odd terms, for +this we use the ~step~ and ~offset~ arguments. ~step~ indicate the distance +between consecutive elements of the input to the fft. If we wanted to compute +the FFT of the even terms, this distance would be 2. For the fft inside this +outer fft of the even entries we also want to split further into even an odd +members, this ~step~ needs to be 4 and so on. ~offset~ represents how many +elements must be skipped until the first element of the subarray. + +Notice that we are allocating in each step a new array to return. The FFT can +however be computed with a single array allocation. For that we store the values +of \[ E_k \] in the first half of the return array and the values of \[ O_k \] +in the second half. Finally we use the so called "butterfly" operations to +compute the final result. + +#+caption: Diagram of the abstract data flow inside the FFT. Observe how the FFT itself is define recursively halving the number of items in each step +[[/svg/fft.svg]] + +#+begin_src common-lisp + (defun fft* (sequence) + (let ((result (make-array (length sequence)))) + (labels ((fft-step (n step offset roffset) + (if (= n 1) + (setf (aref result roffset) (aref sequence offset)) + (let ((n2 (/ n 2))) + (fft-step n2 (* 2 step) offset roffset) ; e_k + (fft-step n2 (* 2 step) (+ offset step) (+ roffset n2)) ; o_k + (loop :for k :below n2 + :for j = (+ roffset k) + :for c = (exp (- (/ (* 2 pi (complex 0d0 1d0) k) n))) + :for ek = (aref result j) + :for ok = (aref result (+ j n2)) + :do + (setf (aref result j) (+ ek (* c ok)) + (aref result (+ j n2)) (- ek (* c ok)))))))) + (fft-step (length sequence) 1 0 0) + result))) +#+end_src + +There are still many optimizations that could be made to speed up our code, like +twiddle factors, cache optimizations and more. But this goes beyond the scope I +planned for this article. Many would probably be also beyond the scope of my +current abilities. |