Skip to content

Commit 91d5625

Browse files
committed
Added Cooley-Tukey in Lisp
1 parent 7bc0e22 commit 91d5625

File tree

3 files changed

+110
-0
lines changed

3 files changed

+110
-0
lines changed

CONTRIBUTORS.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -54,3 +54,4 @@ This file lists everyone, who contributed to this repo and wanted to show up her
5454
- Amaras
5555
- Jonathan Dönszelmann
5656
- Ishaan Verma
57+
- Delphi1024
Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
2+
(defun coefficient (time-index freq-index dft-len)
3+
"Calculates a single twiddle factor for the Fourier Transform."
4+
(exp (- (/ (* #c( 0 1) 2.0 pi time-index freq-index)
5+
dft-len))))
6+
7+
(defun dft (data)
8+
"Performs the Discrete Fourier Transform"
9+
(let ((dft-len (length data)))
10+
(loop for freq-index from 0 below dft-len collect
11+
(loop for time-index from 0 below dft-len sum
12+
(* (coefficient time-index freq-index dft-len) (elt data time-index))))))
13+
14+
(defun merge-sub-ffts (evens odds)
15+
"Combines the FFTs of the even and odd indices"
16+
(let* ((fft-length (+ (length evens) (length odds)))
17+
(indices (loop for i from 0 below (length odds) collect i))
18+
(odd-terms (mapcar #'(lambda (v i) (* v (coefficient 1.0 i fft-length)))
19+
odds
20+
indices)))
21+
(concatenate 'list
22+
(mapcar #'+ evens odd-terms)
23+
(mapcar #'- evens odd-terms))))
24+
25+
(defun cooley-tukey-rec (data)
26+
"Performs the DFT using the recursive Cooley-Tukey method."
27+
(if (<= (length data) 1)
28+
data
29+
(loop for i from 0 below (length data)
30+
if (evenp i)
31+
collect (elt data i) into evens
32+
else
33+
collect (elt data i) into odds
34+
finally
35+
(return (merge-sub-ffts (cooley-tukey-rec evens)
36+
(cooley-tukey-rec odds))))))
37+
38+
(defun reverse-bits (value num-bits)
39+
"Reverses the bits of a value"
40+
(if (= num-bits 1)
41+
value
42+
(let* ((num-low-bits (floor (/ num-bits 2)))
43+
(num-high-bits (- num-bits num-low-bits))
44+
(bit-mask (- (expt 2 num-low-bits) 1))
45+
(lower-half (logand value bit-mask))
46+
(upper-half (ash value (- num-low-bits))))
47+
(logior (ash (reverse-bits lower-half num-low-bits) num-high-bits)
48+
(reverse-bits upper-half num-high-bits)))))
49+
50+
(defun bit-shuffle-indices (data)
51+
"Rearanges the elements in a list according to their bit-reversed indices."
52+
(loop
53+
with num-bits = (floor (log (length data) 2))
54+
for i from 0 below (length data)
55+
collect (elt data (reverse-bits i num-bits))))
56+
57+
(defun butterfly (a b coeff)
58+
"Calculates a single butterfly."
59+
(values (+ a (* coeff b)) (- a (* coeff b))))
60+
61+
(defun butterfly-group (data start stride)
62+
"Calculates a single group of butterflies."
63+
(dotimes (i stride)
64+
(let* ((a-index (+ start i))
65+
(b-index (+ start i stride))
66+
(a (elt data a-index))
67+
(b (elt data b-index))
68+
(coeff (coefficient 1.0 i (* 2 stride))))
69+
(multiple-value-bind (sum difference) (butterfly a b coeff)
70+
(setf (elt data a-index) sum)
71+
(setf (elt data b-index) difference)))))
72+
73+
(defun cooley-tukey-tailrec (data stride)
74+
"Actual iterative implementation of the Cooley-Tukey method."
75+
(if (>= stride (length data))
76+
data
77+
(progn
78+
(loop for i from 0 below (length data) by (* 2 stride) do
79+
(butterfly-group data i stride))
80+
(cooley-tukey-tailrec data (* 2 stride)))))
81+
82+
(defun cooley-tukey-iter (data)
83+
"Performs the DFT using the iterative Cooley-Tukey method."
84+
(cooley-tukey-tailrec (bit-shuffle-indices data) 1))
85+
86+
(defun approx-eql (list1 list2)
87+
(let ((diffs (mapcar #'(lambda (e1 e2) (abs (- e1 e2)))
88+
list1
89+
list2)))
90+
(loop for d in diffs always (< d 1e-9))))
91+
92+
(defun test-fft (data)
93+
(let ((dft-result (dft data))
94+
(rec-result (cooley-tukey-rec data))
95+
(iter-result (cooley-tukey-iter data)))
96+
(format T "~&DFT and recursive Cooley-Tukey approx. equal: ~a"
97+
(approx-eql dft-result rec-result))
98+
(format T "~&DFT and iterative Cooley-Tukey approx. equal: ~a"
99+
(approx-eql dft-result iter-result))
100+
(format T "~&Recursive Cooley-Tukey and iterative Cooley-Tukey approx. equal: ~a"
101+
(approx-eql rec-result iter-result))))
102+
103+
(test-fft '(0.0 0.25 0.5 0.75 0.0 -0.25 -0.5 -0.75))

contents/cooley_tukey/cooley_tukey.md

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -89,6 +89,8 @@ For some reason, though, putting code to this transformation really helped me fi
8989
[import:3-15, lang:"javascript"](code/javascript/fft.js)
9090
{% sample lang="rs" %}
9191
[import:24-37, lang:"rust"](code/rust/fft.rs)
92+
{% sample lang="lisp" %}
93+
[import:2-12, lang:"lisp"](code/clisp/fft.lisp)
9294
{% endmethod %}
9395

9496
In this function, we define `n` to be a set of integers from $$0 \rightarrow N-1$$ and arrange them to be a column.
@@ -142,6 +144,8 @@ In the end, the code looks like:
142144
[import:17-39, lang="javascript"](code/javascript/fft.js)
143145
{% sample lang="rs" %}
144146
[import:39-55, lang:"rust"](code/rust/fft.rs)
147+
{% sample lang="lisp" %}
148+
[import:14-36, lang:"lisp"](code/clisp/fft.lisp)
145149
{% endmethod %}
146150

147151
As a side note, we are enforcing that the array must be a power of 2 for the operation to work.
@@ -255,6 +259,8 @@ Some rather impressive scratch code was submitted by Jie and can be found here:
255259
[import, lang:"javascript"](code/javascript/fft.js)
256260
{% sample lang="rs" %}
257261
[import, lang:"rust"](code/rust/fft.rs)
262+
{% sample lang="lisp" %}
263+
[import, lang:"lisp"](code/clisp/fft.lisp)
258264
{% endmethod %}
259265

260266
<script>

0 commit comments

Comments
 (0)