-
Notifications
You must be signed in to change notification settings - Fork 20
/
fft-recursive.lisp
78 lines (73 loc) · 3.16 KB
/
fft-recursive.lisp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
(defpackage :cp/fft-recursive
(:use :cl)
(:export #:fft-float #:dft! #:inverse-dft! #:convolve!)
(:documentation "Provides FFT by simple recursion. You should use cp/fft
instead. I keep it just for my reference."))
(in-package :cp/fft-recursive)
(deftype fft-float () 'double-float)
(defun dft! (f)
(declare ((simple-array (complex fft-float) (*)) f))
(prog1 f
(let ((n (length f))) ; must be power of 2
(unless (= n 1)
(let ((f0 (make-array (floor n 2) :element-type '(complex fft-float)))
(f1 (make-array (floor n 2) :element-type '(complex fft-float))))
(dotimes (i (floor n 2))
(setf (aref f0 i) (aref f (* i 2))
(aref f1 i) (aref f (+ (* i 2) 1))))
(dft! f0)
(dft! f1)
(let ((zeta (cis (/ #.(coerce (* 2 pi) 'fft-float) n)))
(power-zeta #.(coerce #c(1d0 0d0) '(complex fft-float))))
(declare ((complex fft-float) power-zeta))
(dotimes (i n)
(let ((subindex (mod i (ash n -1))))
(setf (aref f i)
(+ (aref f0 subindex)
(* power-zeta (aref f1 subindex))))
(setf power-zeta (* power-zeta zeta))))))))))
(defun inverse-dft! (f)
(declare ((simple-array (complex fft-float) (*)) f))
(labels ((%idft! (f)
(declare ((simple-array (complex fft-float) (*)) f))
(prog1 f
(let ((n (length f))) ; must be power of 2
(unless (= n 1)
(let ((f0 (make-array (floor n 2) :element-type '(complex fft-float)))
(f1 (make-array (floor n 2) :element-type '(complex fft-float))))
(dotimes (i (floor n 2))
(setf (aref f0 i) (aref f (* i 2))
(aref f1 i) (aref f (+ (* i 2) 1))))
(%idft! f0)
(%idft! f1)
(let ((zeta (cis (/ #.(coerce (* -2 pi) 'fft-float) n)))
(power-zeta #.(coerce #c(1d0 0d0) '(complex fft-float))))
(declare ((complex fft-float) power-zeta))
(dotimes (i n)
(let ((subindex (mod i (ash n -1))))
(setf (aref f i)
(+ (aref f0 subindex)
(* power-zeta (aref f1 subindex))))
(setf power-zeta (* power-zeta zeta)))))))))))
(let* ((n (length f))
(/n (/ (coerce n 'fft-float))))
(%idft! f)
(dotimes (i n f)
(setf (aref f i) (* (aref f i) /n))))))
(declaim (inline power2-p))
(defun power2-p (x)
"Checks if X is a power of 2."
(zerop (logand x (- x 1))))
(declaim (inline convolve!))
(defun convolve! (g h)
(declare ((simple-array (complex fft-float) (*)) g h))
(assert (and (power2-p (length g))
(power2-p (length h))
(= (length g) (length h))))
(let ((n (length g)))
(dft! g)
(dft! h)
(let ((f (make-array n :element-type '(complex fft-float))))
(dotimes (i n)
(setf (aref f i) (* (aref g i) (aref h i))))
(inverse-dft! f))))