1 ;;;; -*- Mode: Lisp; Syntax: ANSI-Common-Lisp; Base: 10 -*-
2 ;;;; *************************************************************************
3 ;;;; FILE IDENTIFICATION
5 ;;;; Name: transform.lisp
6 ;;;; Purpose: Transformation functions
7 ;;;; Programmer: Kevin M. Rosenberg
8 ;;;; Date Started: March 2009
10 ;;;; This file and CL-FFTW3 are Copyright (c) 2009-2011 by Kevin M. Rosenberg
12 ;;;; FFTW3 users are granted the rights to distribute and use this software
13 ;;;; as governed by the terms of the Lisp Lesser GNU Public License
14 ;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL.
15 ;;;; *************************************************************************
19 (defconstant +size-double+ 8)
21 (defun fftw-r2hc-1d (in &key (start 0) (count nil) (normalize t))
22 (declare #.*standard-optimize-settings*)
23 (multiple-value-bind (start count) (normalize-input-range in start count)
24 (declare (fixnum start count))
27 (out-c (make-array count :element-type 'double-float))
28 (in-cf (fftw-malloc (* count +size-double+)))
29 (out-cf (fftw-malloc (* count +size-double+))))
30 (declare (fixnum in-pos))
33 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
36 (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-r2hc+
37 (logior +fftw-estimate+ +fftw-destroy-input+))))
39 (fftw-destroy-plan p))
42 (let ((factor (coerce (/ 1 count) 'double-float)))
43 (declare (double-float factor))
46 (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i)))))
49 (setf (aref out-c i) (cffi:mem-aref out-cf :double i))))
56 (defun make-fftw-r2r-1d-multi (count type)
57 (declare #.*standard-optimize-settings*)
58 (declare (fixnum count type))
59 (let* ((out-c (make-array count :element-type 'double-float))
60 (in-cf (fftw-malloc (* count +size-double+)))
61 (out-cf (fftw-malloc (* count +size-double+)))
62 (plan (fftw-plan-r2r-1d count in-cf out-cf type
63 (logior +fftw-estimate+ +fftw-destroy-input+))))
64 (make-instance 'fftw-r2r-1d-multi
65 :in-count count :out-c out-c :plan plan
66 :in-cf in-cf :out-cf out-cf)))
68 (defun fftw-r2r-1d-multi (multi in &key (start 0) (count nil) (normalize t))
69 (declare #.*standard-optimize-settings*)
70 (multiple-value-bind (start count) (normalize-input-range in start count)
71 (declare (fixnum start count))
72 (unless (equal (in-count multi) count)
73 (error "Different plan and vector lengths."))
75 (let ((in-cf (in-cf multi))
76 (out-cf (out-cf multi))
79 (declare (fixnum in-pos))
83 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
86 (fftw-execute (plan multi))
88 (let ((factor (coerce (/ 1 count) 'double-float)))
89 (declare (double-float factor))
92 (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i)))))
95 (setf (aref out-c i) (cffi:mem-aref out-cf :double i))))
100 (defun fftw-hc2r-1d (in)
101 (declare #.*standard-optimize-settings*)
102 (let* ((count (length in))
103 (out-c (make-array count :element-type 'double-float))
104 (in-cf (fftw-malloc (* count +size-double+)))
105 (out-cf (fftw-malloc (* count +size-double+))))
106 (declare (fixnum count))
110 (setf (cffi:mem-aref in-cf :double i) (aref in i)))
112 (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-hc2r+
113 (logior +fftw-estimate+ +fftw-destroy-input+))))
115 (fftw-destroy-plan p))
119 (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))
124 (defun fftw-r2c-1d (in &key (start 0) (count nil) (normalize t))
125 (declare #.*standard-optimize-settings*)
126 (multiple-value-bind (start count) (normalize-input-range in start count)
127 (declare (fixnum start count))
129 (let* ((out-n (1+ (floor count 2)))
130 (in-cf (fftw-malloc (* count +size-double+)))
131 (out-cf (fftw-malloc (* 2 out-n +size-double+)))
133 (declare (fixnum in-pos out-n))
137 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
140 (let ((p (fftw-plan-dft-r2c-1d count in-cf out-cf
141 (logior +fftw-estimate+ +fftw-destroy-input+))))
143 (fftw-destroy-plan p))
145 (let ((out (make-array out-n :element-type '(or rational complex)))
147 (declare (fixnum opos))
149 (let ((factor (coerce (/ 1 count) 'double-float)))
150 (declare (double-float factor))
153 (let ((impos (1+ opos)))
154 (declare (fixnum impos))
155 (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor)
156 (* (cffi:mem-aref out-cf :double impos) factor)))
160 (let ((impos (1+ opos)))
161 (declare (fixnum impos))
162 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
163 (cffi:mem-aref out-cf :double impos)))
169 (defun make-fftw-r2c-1d-multi (count &key (complex-output t))
170 (declare #.*standard-optimize-settings*)
171 (declare (fixnum count))
172 (let* ((out-n (1+ (floor count 2)))
173 (in-cf (fftw-malloc (* count +size-double+)))
174 (out-cf (fftw-malloc (* 2 out-n +size-double+)))
175 (plan (fftw-plan-dft-r2c-1d count in-cf out-cf
176 (logior +fftw-estimate+ +fftw-destroy-input+))))
177 (declare (fixnum out-n))
179 (unless complex-output
182 (make-instance 'fftw-r2c-1d-multi
183 :in-count count :out-count out-n :plan plan
186 :complex-output complex-output
189 (make-array out-n :element-type '(or double-float complex) :initial-element 0d0)
190 (make-array out-n :element-type 'double-float)))))
192 (defgeneric destroy-fftw-multi (multi))
193 (defmethod destroy-fftw-multi ((multi fftw-multi))
194 (declare #.*standard-optimize-settings*)
196 (fftw-destroy-plan (plan multi))
197 (setf (plan multi) nil))
198 (dolist (slot '(in-cf out-cf))
199 (when (slot-value multi slot)
200 (fftw-free (slot-value multi slot))
201 (setf (slot-value multi slot) nil)))
202 (dolist (slot '(plan-even plan-odd))
203 (when (and (slot-exists-p multi slot)
204 (slot-value multi slot))
205 (fftw-destroy-plan (slot-value multi slot))
206 (setf (slot-value multi slot) nil))))
209 (defun fftw-r2c-1d-multi (multi in &key (start 0) (normalize t))
210 (declare #.*standard-optimize-settings*)
211 (multiple-value-bind (start count) (normalize-input-range in start (in-count multi))
212 (declare (fixnum start count))
214 (unless (equal (length in) (in-count multi))
215 (error "In count of multi plan doesn't equal length of in vector."))
217 (let ((in-cf (in-cf multi))
218 (out-cf (out-cf multi))
222 (in-count (in-count multi))
223 (out-count (out-count multi)))
224 (declare (fixnum in-pos opos in-count out-count))
226 (dotimes (i in-count)
228 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
231 (fftw-execute (plan multi))
235 ((complex-output multi)
237 (let ((factor (coerce (/ 1 count) 'double-float)))
238 (declare (double-float factor))
239 (dotimes (i out-count)
241 (let ((impos (1+ opos)))
242 (declare (fixnum impos))
243 (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor)
244 (* (cffi:mem-aref out-cf :double impos) factor)))
246 (dotimes (i out-count)
248 (let ((impos (1+ opos)))
249 (declare (fixnum impos))
250 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
251 (cffi:mem-aref out-cf :double impos)))
255 (let ((o2 (* out-count 2))
256 (factor (coerce (/ 1 count) 'double-float)))
258 (double-float factor))
261 (setf (aref out i) (* factor (aref out i))))))))
264 (defun fftw-c2r-1d (in)
265 (declare #.*standard-optimize-settings*)
266 (let* ((count (length in))
267 (out-n (if (complex-is-even-p in)
269 (1+ (* 2 (1- count)))))
270 (out-c (make-array out-n :element-type 'double-float))
271 (in-cf (fftw-malloc (* 2 count +size-double+)))
272 (out-cf (fftw-malloc (* 2 count +size-double+)))
274 (declare (fixnum count out-n pos))
279 (let ((c (aref in i)))
280 (declare (complex c))
281 (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float))
283 (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float))
286 (let ((p (fftw-plan-dft-c2r-1d out-n in-cf out-cf
287 (logior +fftw-estimate+ +fftw-destroy-input+))))
289 (fftw-destroy-plan p))
292 (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))
298 (defun make-fftw-c2r-1d-multi (count)
299 (declare #.*standard-optimize-settings*)
300 (declare (fixnum count))
301 (let* ((out-count-even (* 2 (1- count)))
302 (out-count-odd (1+ (* 2 (1- count))))
303 (in-cf (fftw-malloc (* 2 count +size-double+)))
304 (out-cf (fftw-malloc (* 2 (max out-count-odd out-count-even) +size-double+))) ;; use larger size
305 (plan-odd (fftw-plan-dft-c2r-1d out-count-odd in-cf out-cf
306 (logior +fftw-estimate+ +fftw-destroy-input+)))
307 (plan-even (fftw-plan-dft-c2r-1d out-count-even in-cf out-cf +fftw-measure+)))
308 (declare (fixnum out-count-even out-count-odd))
310 (make-instance 'fftw-c2r-1d-multi
316 :out-count-even out-count-even
317 :out-count-odd out-count-odd
318 :out-even (make-array out-count-even :element-type 'double-float)
319 :out-odd (make-array out-count-odd :element-type 'double-float))))
321 (defun fftw-c2r-1d-multi (multi in)
322 (declare #.*standard-optimize-settings*)
323 (let* ((is-even (complex-is-even-p in))
324 (out (if is-even (out-even multi) (out-odd multi)))
325 (out-count (if is-even (out-count-even multi) (out-count-odd multi)))
326 (plan (if is-even (plan-even multi) (plan-odd multi)))
327 (count (in-count multi))
328 (in-cf (in-cf multi))
329 (out-cf (out-cf multi))
331 (declare (fixnum count out-count pos))
336 (let ((c (aref in i)))
337 (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float))
339 (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float))
344 (dotimes (i out-count)
346 (setf (aref out i) (cffi:mem-aref out-cf :double i)))
350 (defun fftw-c-1d (in &key (direction +fftw-forward+)
351 (start 0) (count nil) (normalize t))
352 (declare #.*standard-optimize-settings*)
353 (multiple-value-bind (start count) (normalize-input-range in start count)
354 (declare (fixnum start count))
357 (out (make-array out-n :element-type '(or rational complex)))
358 (in-cf (fftw-malloc (* 2 count +size-double+)))
359 (out-cf (fftw-malloc (* 2 out-n +size-double+)))
362 (declare (fixnum out-n ipos opos))
368 (let ((c (aref in ipos)))
370 (setf (cffi:mem-aref in-cf :double opos) (coerce (realpart c) 'double-float))
372 (setf (cffi:mem-aref in-cf :double opos) (coerce (imagpart c) 'double-float))
375 (let ((p (fftw-plan-dft-1d count in-cf out-cf direction
376 (logior +fftw-estimate+ +fftw-destroy-input+))))
378 (fftw-destroy-plan p))
382 (let ((factor (coerce (/ 1 count) 'double-float)))
383 (declare (double-float factor))
386 (let ((impos (1+ opos)))
387 (declare (fixnum impos))
388 (setf (aref out i) (complex (* factor (cffi:mem-aref out-cf :double opos))
389 (* factor (cffi:mem-aref out-cf :double impos)))))
393 (let ((impos (1+ opos)))
394 (declare (fixnum impos))
395 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
396 (cffi:mem-aref out-cf :double impos))))