-
Notifications
You must be signed in to change notification settings - Fork 6
Expand file tree
/
Copy pathforward.lisp
More file actions
156 lines (151 loc) · 6.89 KB
/
forward.lisp
File metadata and controls
156 lines (151 loc) · 6.89 KB
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
(in-package "NAPA-FFT.GEN")
(defvar *fwd-base-case* 32)
(defun %gen-flat-dif (n scale window)
(with-vector (n)
(let ((last n))
(labels ((rec (start n)
(cond
((= n 1)
(when (and (= n last)
(or (not (eql scale 1d0))
window))
(scale start scale window)))
((= n 2)
(when (and (= n last)
(or (not (eql scale 1d0))
window))
(scale start scale window)
(scale (1+ start) scale window))
(butterfly start (1+ start)))
(t
(let* ((n/2 (truncate n 2))
(start2 (+ start n/2))
(n/4 (truncate n/2 2))
(start3 (+ start2 n/4)))
(dotimes (i n/2)
(when (= last n)
(scale (+ i start) scale window)
(scale (+ i start2) scale window))
(butterfly (+ i start)
(+ i start2)))
(rec start n/2)
(dotimes (count n/4)
(let ((i (+ count start2))
(j (+ count start3))
(k (+ n/2 +twiddle-offset+
(* 2 count))))
(rotate j nil 3/4)
(butterfly i j)
(rotate i k (/ (* -1 count) n))
(rotate j (1+ k) (/ (* -3 count) n))))
(rec start2 n/4)
(rec start3 n/4))))))
(rec 0 n)))))
(defun gen-flat-dif (n &key (scale 1d0) window)
(let ((table (load-time-value (make-hash-table :test #'equal)))
(key (list n scale window)))
(or (gethash key table)
(setf (gethash key table)
(%gen-flat-dif n scale window)))))
(defun gen-dif (n &key (scale 1d0) window)
(let ((defs '())
(last n))
(labels ((name (n)
(intern (format nil "~A/~A" 'dif n)
"NAPA-FFT.GEN"))
(gen (n &aux (name (name n)))
(when (member name defs :key #'first)
(return-from gen))
(cond
((<= n *fwd-base-case*)
(push
`(,(name n) (start)
(declare (type index start)
(ignorable start))
,(if (= n last)
(gen-flat-dif n :scale scale :window window)
(gen-flat-dif n)))
defs))
(t
(gen (truncate n 4))
(gen (truncate n 2))
(let* ((n/2 (truncate n 2))
(n/4 (truncate n 4))
(name/2 (name n/2))
(name/4 (name n/4))
(body
`(,(name n) (start)
(declare (type index start))
(for (,n/2 (i start)
,@(and (= n last)
window
`((k window-start))))
(let ((x ,(if (= n last)
`(%window (%scale (aref vec i) ,scale)
,window
,(if window 'k 0))
`(aref vec i)))
(y ,(if (= n last)
`(%window (%scale (aref vec (+ i ,n/2))
,scale)
,window
,(if window `(+ k ,n/2) 0))
`(aref vec (+ i ,n/2)))))
(setf (aref vec i) (+ x y)
(aref vec (+ i ,n/2)) (- x y))))
(,name/2 start)
(for (,n/4 (i start)
(k ,(+ n/2 +twiddle-offset+) 2))
(let ((x (aref vec (+ i ,n/2)))
(y (mul-i (aref vec (+ i ,(+ n/2 n/4)))))
(t1 (aref twiddle k))
(t2 (aref twiddle (1+ k))))
(setf (aref vec (+ i ,n/2))
(* t1 (+ x y))
(aref vec (+ i ,(+ n/2 n/4)))
(* t2 (- x y)))))
(,name/4 (+ start ,n/2))
(,name/4 (+ start ,(+ n/2 n/4))))))
(push body defs))))))
(gen n)
`(labels (,@(nreverse defs))
(declare (inline ,(name n)))
(,(name n) start)))))
(defun %dif (vec start n twiddle)
(declare (type complex-sample-array vec twiddle)
(type index start)
(type size n))
(labels ((rec (start n)
(declare (type index start)
(type size n))
(cond ((>= n 4)
(let* ((n/2 (truncate n 2))
(start2 (+ start n/2))
(n/4 (truncate n/2 2))
(start3 (+ start2 n/4)))
(for (n/2 (i start)
(j start2))
(let ((x (aref vec i))
(y (aref vec j)))
(setf (aref vec i) (+ x y)
(aref vec j) (- x y))))
(rec start n/2)
(for (n/4 (i start2)
(j start3)
(k (+ n/2 +twiddle-offset+) 2))
(let ((x (aref vec i))
(y (mul-i (aref vec j)))
(t1 (aref twiddle k))
(t2 (aref twiddle (1+ k))))
(setf (aref vec i) (* t1 (+ x y))
(aref vec j) (* t2 (- x y)))))
(rec start2 n/4)
(rec start3 n/4)))
((= n 2)
(let ((s0 (aref vec start))
(s1 (aref vec (1+ start))))
(setf (aref vec start) (+ s0 s1)
(aref vec (1+ start)) (- s0 s1)))
nil))))
(rec start n)
vec))