Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Newer
Older
100644 170 lines (159 sloc) 6.226 kB
910f46a @pkhuong Real transforms
authored
1 (in-package "NAPA-FFT.IMPL")
2
3 (declaim (inline swap))
4 (defun swap (x)
5 (declare (type complex-sample x))
6 #+(and sbcl complex-float-vops)
7 (sb-vm::swap-complex x)
8 (complex (imagpart x) (realpart x)))
9
10 (declaim (inline fft-swizzled-reals))
11 (defun fft-swizzled-reals (vec scale)
12 (declare (type complex-sample-array vec))
13 (let* ((size (length vec))
14 (n/2 (truncate size 2)))
15 (fft vec :dst vec :size n/2 :scale scale)
16 (let ((x (aref vec 0)))
17 (setf (aref vec n/2) (complex (imagpart x))
18 (aref vec 0) (complex (realpart x))))
19 (loop for i of-type index from 1 upto (truncate n/2 2)
20 for j of-type index from (1+ n/2)
21 do (let* ((x (aref vec i))
22 (xs (swap x))
23 (y (aref vec (- n/2 i)))
24 (ys (swap y))
25 (xj (* .5d0 (+ (conjugate xs) ys)))
26 (-xj (* .5d0 (+ (conjugate ys) xs))))
27 (setf (aref vec j) xj
28 (aref vec (- size i)) -xj
29 (aref vec i) (- x (napa-fft.gen::mul+i xj))
30 (aref vec (- n/2 i)) (- y (napa-fft.gen::mul+i -xj)))))
31 vec))
32
33 (defvar *rfft-twiddles* nil)
34 (defvar *rifft-twiddles* nil)
35
36 (defun %get-radix-2-twiddle (n direction)
37 (let* ((n/2 (truncate n 2))
38 (twiddle (make-array n/2 :element-type 'complex-sample)))
39 (let ((root (if (plusp direction)
40 (complex 1d0 0d0)
41 (complex 0d0 1d0)))
42 (mul (exp (* -2 pi (complex 0d0 direction) (/ 1d0 n)))))
43 (declare (type complex-sample root mul))
44 (loop for k below n/2 do
45 (setf (aref twiddle k) root
46 root (* root mul)))
47 twiddle)))
48
49 (defun get-radix-2-twiddle (n direction)
50 (assert (power-of-two-p n))
51 (assert (member direction '(-1d0 1d0)))
52 (let ((twiddles (ecase direction
53 (1d0
54 (or *rfft-twiddles*
55 (setf *rfft-twiddles*
56 (make-array 32 :initial-element nil))))
57 (-1d0
58 (or *rifft-twiddles*
59 (setf *rifft-twiddles*
60 (make-array 32 :initial-element nil))))))
61 (len (lb n)))
62 (or (aref twiddles len)
63 (setf (aref twiddles len)
64 (%get-radix-2-twiddle n direction)))))
65
66 (defmacro with-scale ((scale) &body body)
67 `(ecase ,scale
68 ((nil 1)
69 (flet ((scale (x)
70 x))
71 (declare (inline scale))
72 ,@body))
73 ((t :inv)
74 (flet ((scale (x)
75 (* x .5d0)))
76 (declare (inline scale))
77 ,@body))
78 ((:sqrt sqrt)
79 (flet ((scale (x)
80 (* x ,(sqrt .5d0))))
81 (declare (inline scale))
82 ,@body))))
83
edea1a6 @pkhuong DST in R[I]FFT is keyword argument
authored
84 (defun rfft (vec &key dst
d85007e @pkhuong Use new real-sample[-array] types
authored
85 size
86 (scale nil))
87 (declare (type scaling scale)
910f46a @pkhuong Real transforms
authored
88 (optimize speed))
d85007e @pkhuong Use new real-sample[-array] types
authored
89 (let* ((vec (real-samplify vec))
90 (size (or size (length vec))))
91 (declare (type real-sample-array vec)
92 (type size size))
93 (assert (>= (length vec) size))
94 (assert (power-of-two-p size))
95 (let* ((n size)
96 (n/2 (truncate n 2))
97 (dst (or dst
98 (make-array n :element-type 'complex-sample)))
99 (twiddle (get-radix-2-twiddle n 1d0)))
100 (declare (type complex-sample-array twiddle)
101 (type complex-sample-array dst))
102 (assert (>= (length dst) size))
103 (locally (declare (optimize (safety 0)))
104 (with-scale (scale)
105 (loop for i of-type index below n by 2
106 for j of-type index from 0
107 do (setf (aref dst j)
108 (scale (complex (aref vec i)
109 (aref vec (+ i 1)))))))
110 (fft-swizzled-reals dst scale)
111 (loop for i of-type index below n/2
112 for j of-type index from n/2
113 do
114 (let* ((x (aref dst i))
115 (y (* (aref dst j) (aref twiddle i))))
116 (setf (aref dst i) (+ x y)
117 (aref dst j) (- x y)))))
118 dst)))
910f46a @pkhuong Real transforms
authored
119
120 (defun windowed-rfft (signal-vector center length
121 &key (window-fn 'hann)
122 dst
123 (scale nil))
d85007e @pkhuong Use new real-sample[-array] types
authored
124 (declare (type index length)
910f46a @pkhuong Real transforms
authored
125 (optimize speed))
126 (assert (power-of-two-p length))
d85007e @pkhuong Use new real-sample[-array] types
authored
127 (let* ((signal-vector (real-samplify signal-vector))
128 (input-window (extract-centered-window signal-vector center length))
910f46a @pkhuong Real transforms
authored
129 (window (window-vector window-fn length)))
d85007e @pkhuong Use new real-sample[-array] types
authored
130 (declare (type real-sample-array signal-vector input-window window))
910f46a @pkhuong Real transforms
authored
131 (map-into input-window #'* input-window window)
edea1a6 @pkhuong DST in R[I]FFT is keyword argument
authored
132 (rfft input-window
133 :dst dst
910f46a @pkhuong Real transforms
authored
134 :scale scale)))
135
d85007e @pkhuong Use new real-sample[-array] types
authored
136 (defun rifft (vec &key dst size (scale t))
137 (declare (type scaling scale)
910f46a @pkhuong Real transforms
authored
138 (optimize speed))
d85007e @pkhuong Use new real-sample[-array] types
authored
139 (let* ((vec (complex-samplify vec))
140 (size (or size (length vec))))
141 (declare (type size size))
142 (assert (>= (length vec) size))
143 (assert (power-of-two-p size))
144 (let* ((n size)
145 (n/2 (truncate n 2))
146 (twiddle (get-radix-2-twiddle n -1d0))
147 (dst (or dst
148 (make-array n :element-type 'double-float))))
149 (declare (type complex-sample-array twiddle)
150 (type real-sample-array dst))
151 (assert (>= (length dst) size))
152 (locally (declare (optimize (safety 0)))
153 (loop for i of-type index below n/2
154 for j of-type index from n/2
155 do (let* ((x (aref vec i))
156 (y (aref vec j))
157 (sum (+ x y))
158 (sub (* (- x y)
159 (aref twiddle i))))
160 (setf (aref vec i) (+ sum sub))))
161 (ifft vec :dst vec :size n/2 :scale scale)
162 (with-scale (scale)
163 (loop for i of-type index below n/2
164 for j of-type index by 2
165 do (let ((x (scale (aref vec i))))
166 (setf (aref dst j) (realpart x)
167 (aref dst (1+ j)) (imagpart x))))))
168 dst)))
910f46a @pkhuong Real transforms
authored
169
Something went wrong with that request. Please try again.