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
12 ;;;; This file, part of FFTW3, is Copyright (c) 2009 by Kevin M. Rosenberg
14 ;;;; FFTW3 users are granted the rights to distribute and use this software
15 ;;;; as governed by the terms of the Lisp Lesser GNU Public License
16 ;;;; (http://opensource.franz.com/preamble.html), also known as the LLGPL.
17 ;;;; *************************************************************************
21 (defconstant +size-double+ 8)
23 (defun fftw-r2hc-1d (in &key (start 0) (count nil) (normalize t))
24 (declare #.*standard-optimize-settings*)
25 (multiple-value-bind (start count) (normalize-input-range in start count)
26 (declare (fixnum start count))
29 (out-c (make-array count :element-type 'double-float))
30 (in-cf (fftw-malloc (* count +size-double+)))
31 (out-cf (fftw-malloc (* count +size-double+))))
32 (declare (fixnum in-pos))
35 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
38 (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-r2hc+
39 (logior +fftw-estimate+ +fftw-destroy-input+))))
41 (fftw-destroy-plan p))
44 (let ((factor (coerce (/ 1 count) 'double-float)))
45 (declare (double-float factor))
48 (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i)))))
51 (setf (aref out-c i) (cffi:mem-aref out-cf :double i))))
58 (defun make-fftw-r2r-1d-multi (count type)
59 (declare #.*standard-optimize-settings*)
60 (declare (fixnum count type))
61 (let* ((out-c (make-array count :element-type 'double-float))
62 (in-cf (fftw-malloc (* count +size-double+)))
63 (out-cf (fftw-malloc (* count +size-double+)))
64 (plan (fftw-plan-r2r-1d count in-cf out-cf type
65 (logior +fftw-estimate+ +fftw-destroy-input+))))
66 (make-instance 'fftw-r2r-1d-multi
67 :in-count count :out-c out-c :plan plan
68 :in-cf in-cf :out-cf out-cf)))
70 (defun fftw-r2r-1d-multi (multi in &key (start 0) (count nil) (normalize t))
71 (declare #.*standard-optimize-settings*)
72 (multiple-value-bind (start count) (normalize-input-range in start count)
73 (declare (fixnum start count))
74 (unless (equal (in-count multi) count)
75 (error "Different plan and vector lengths."))
77 (let ((in-cf (in-cf multi))
78 (out-cf (out-cf multi))
81 (declare (fixnum in-pos))
85 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
88 (fftw-execute (plan multi))
90 (let ((factor (coerce (/ 1 count) 'double-float)))
91 (declare (double-float factor))
94 (setf (aref out-c i) (* factor (cffi:mem-aref out-cf :double i)))))
97 (setf (aref out-c i) (cffi:mem-aref out-cf :double i))))
102 (defun fftw-hc2r-1d (in)
103 (declare #.*standard-optimize-settings*)
104 (let* ((count (length in))
105 (out-c (make-array count :element-type 'double-float))
106 (in-cf (fftw-malloc (* count +size-double+)))
107 (out-cf (fftw-malloc (* count +size-double+))))
108 (declare (fixnum count))
112 (setf (cffi:mem-aref in-cf :double i) (aref in i)))
114 (let ((p (fftw-plan-r2r-1d count in-cf out-cf +fftw-hc2r+
115 (logior +fftw-estimate+ +fftw-destroy-input+))))
117 (fftw-destroy-plan p))
121 (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))
126 (defun fftw-r2c-1d (in &key (start 0) (count nil) (normalize t))
127 (declare #.*standard-optimize-settings*)
128 (multiple-value-bind (start count) (normalize-input-range in start count)
129 (declare (fixnum start count))
131 (let* ((out-n (1+ (floor count 2)))
132 (in-cf (fftw-malloc (* count +size-double+)))
133 (out-cf (fftw-malloc (* 2 out-n +size-double+)))
135 (declare (fixnum in-pos out-n))
139 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
142 (let ((p (fftw-plan-dft-r2c-1d count in-cf out-cf
143 (logior +fftw-estimate+ +fftw-destroy-input+))))
145 (fftw-destroy-plan p))
147 (let ((out (make-array out-n :element-type '(or rational complex)))
149 (declare (fixnum opos))
151 (let ((factor (coerce (/ 1 count) 'double-float)))
152 (declare (double-float factor))
155 (let ((impos (1+ opos)))
156 (declare (fixnum impos))
157 (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor)
158 (* (cffi:mem-aref out-cf :double impos) factor)))
162 (let ((impos (1+ opos)))
163 (declare (fixnum impos))
164 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
165 (cffi:mem-aref out-cf :double impos)))
171 (defun make-fftw-r2c-1d-multi (count &key (complex-output t))
172 (declare #.*standard-optimize-settings*)
173 (declare (fixnum count))
174 (let* ((out-n (1+ (floor count 2)))
175 (in-cf (fftw-malloc (* count +size-double+)))
176 (out-cf (fftw-malloc (* 2 out-n +size-double+)))
177 (plan (fftw-plan-dft-r2c-1d count in-cf out-cf
178 (logior +fftw-estimate+ +fftw-destroy-input+))))
179 (declare (fixnum out-n))
181 (unless complex-output
184 (make-instance 'fftw-r2c-1d-multi
185 :in-count count :out-count out-n :plan plan
188 :complex-output complex-output
191 (make-array out-n :element-type '(or double-float complex) :initial-element 0d0)
192 (make-array out-n :element-type 'double-float)))))
194 (defgeneric destroy-fftw-multi (multi))
195 (defmethod destroy-fftw-multi ((multi fftw-multi))
196 (declare #.*standard-optimize-settings*)
198 (fftw-destroy-plan (plan multi))
199 (setf (plan multi) nil))
200 (dolist (slot '(in-cf out-cf))
201 (when (slot-value multi slot)
202 (fftw-free (slot-value multi slot))
203 (setf (slot-value multi slot) nil)))
204 (dolist (slot '(plan-even plan-odd))
205 (when (and (slot-exists-p multi slot)
206 (slot-value multi slot))
207 (fftw-destroy-plan (slot-value multi slot))
208 (setf (slot-value multi slot) nil))))
211 (defun fftw-r2c-1d-multi (multi in &key (start 0) (normalize t))
212 (declare #.*standard-optimize-settings*)
213 (multiple-value-bind (start count) (normalize-input-range in start (in-count multi))
214 (declare (fixnum start count))
216 (unless (equal (length in) (in-count multi))
217 (error "In count of multi plan doesn't equal length of in vector."))
219 (let ((in-cf (in-cf multi))
220 (out-cf (out-cf multi))
224 (in-count (in-count multi))
225 (out-count (out-count multi)))
226 (declare (fixnum in-pos opos in-count out-count))
228 (dotimes (i in-count)
230 (setf (cffi:mem-aref in-cf :double i) (coerce (aref in in-pos) 'double-float))
233 (fftw-execute (plan multi))
237 ((complex-output multi)
239 (let ((factor (coerce (/ 1 count) 'double-float)))
240 (declare (double-float factor))
241 (dotimes (i out-count)
243 (let ((impos (1+ opos)))
244 (declare (fixnum impos))
245 (setf (aref out i) (complex (* (cffi:mem-aref out-cf :double opos) factor)
246 (* (cffi:mem-aref out-cf :double impos) factor)))
248 (dotimes (i out-count)
250 (let ((impos (1+ opos)))
251 (declare (fixnum impos))
252 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
253 (cffi:mem-aref out-cf :double impos)))
257 (let ((o2 (* out-count 2))
258 (factor (coerce (/ 1 count) 'double-float)))
260 (double-float factor))
263 (setf (aref out i) (* factor (aref out i))))))))
266 (defun fftw-c2r-1d (in)
267 (declare #.*standard-optimize-settings*)
268 (let* ((count (length in))
269 (out-n (if (complex-is-even-p in)
271 (1+ (* 2 (1- count)))))
272 (out-c (make-array out-n :element-type 'double-float))
273 (in-cf (fftw-malloc (* 2 count +size-double+)))
274 (out-cf (fftw-malloc (* 2 count +size-double+)))
276 (declare (fixnum count out-n pos))
281 (let ((c (aref in i)))
282 (declare (complex c))
283 (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float))
285 (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float))
288 (let ((p (fftw-plan-dft-c2r-1d out-n in-cf out-cf
289 (logior +fftw-estimate+ +fftw-destroy-input+))))
291 (fftw-destroy-plan p))
294 (setf (aref out-c i) (cffi:mem-aref out-cf :double i)))
300 (defun make-fftw-c2r-1d-multi (count)
301 (declare #.*standard-optimize-settings*)
302 (declare (fixnum count))
303 (let* ((out-count-even (* 2 (1- count)))
304 (out-count-odd (1+ (* 2 (1- count))))
305 (in-cf (fftw-malloc (* 2 count +size-double+)))
306 (out-cf (fftw-malloc (* 2 (max out-count-odd out-count-even) +size-double+))) ;; use larger size
307 (plan-odd (fftw-plan-dft-c2r-1d out-count-odd in-cf out-cf
308 (logior +fftw-estimate+ +fftw-destroy-input+)))
309 (plan-even (fftw-plan-dft-c2r-1d out-count-even in-cf out-cf +fftw-measure+)))
310 (declare (fixnum out-count-even out-count-odd))
312 (make-instance 'fftw-c2r-1d-multi
318 :out-count-even out-count-even
319 :out-count-odd out-count-odd
320 :out-even (make-array out-count-even :element-type 'double-float)
321 :out-odd (make-array out-count-odd :element-type 'double-float))))
323 (defun fftw-c2r-1d-multi (multi in)
324 (declare #.*standard-optimize-settings*)
325 (let* ((is-even (complex-is-even-p in))
326 (out (if is-even (out-even multi) (out-odd multi)))
327 (out-count (if is-even (out-count-even multi) (out-count-odd multi)))
328 (plan (if is-even (plan-even multi) (plan-odd multi)))
329 (count (in-count multi))
330 (in-cf (in-cf multi))
331 (out-cf (out-cf multi))
333 (declare (fixnum count out-count pos))
338 (let ((c (aref in i)))
339 (setf (cffi:mem-aref in-cf :double pos) (coerce (realpart c) 'double-float))
341 (setf (cffi:mem-aref in-cf :double pos) (coerce (imagpart c) 'double-float))
346 (dotimes (i out-count)
348 (setf (aref out i) (cffi:mem-aref out-cf :double i)))
352 (defun fftw-c-1d (in &key (direction +fftw-forward+)
353 (start 0) (count nil) (normalize t))
354 (declare #.*standard-optimize-settings*)
355 (multiple-value-bind (start count) (normalize-input-range in start count)
356 (declare (fixnum start count))
359 (out (make-array out-n :element-type '(or rational complex)))
360 (in-cf (fftw-malloc (* 2 count +size-double+)))
361 (out-cf (fftw-malloc (* 2 out-n +size-double+)))
364 (declare (fixnum out-n ipos opos))
370 (let ((c (aref in ipos)))
372 (setf (cffi:mem-aref in-cf :double opos) (coerce (realpart c) 'double-float))
374 (setf (cffi:mem-aref in-cf :double opos) (coerce (imagpart c) 'double-float))
377 (let ((p (fftw-plan-dft-1d count in-cf out-cf direction
378 (logior +fftw-estimate+ +fftw-destroy-input+))))
380 (fftw-destroy-plan p))
384 (let ((factor (coerce (/ 1 count) 'double-float)))
385 (declare (double-float factor))
388 (let ((impos (1+ opos)))
389 (declare (fixnum impos))
390 (setf (aref out i) (complex (* factor (cffi:mem-aref out-cf :double opos))
391 (* factor (cffi:mem-aref out-cf :double impos)))))
395 (let ((impos (1+ opos)))
396 (declare (fixnum impos))
397 (setf (aref out i) (complex (cffi:mem-aref out-cf :double opos)
398 (cffi:mem-aref out-cf :double impos))))