1 (declaim (optimize (speed 2) (debug 3) (safety 3)))
4 (:use
:cl
:sb-alien
:sb-c-call
)
9 (in-package :fftwf-fft
)
11 (load-shared-object "/usr/lib/libfftw3f.so")
13 ;; multithreading for fftw is just a matter of two initializing
14 ;; function calls, see:
15 ;; http://www.fftw.org/fftw3_doc/Usage-of-Multi_002dthreaded-FFTW.html#Usage-of-Multi_002dthreaded-FFTW
16 ;; but it doesn't get significantly faster
20 (load-shared-object "/usr/lib/libfftw3f_threads.so")
22 (define-alien-routine ("fftwf_init_threads" init-threads
)
24 (define-alien-routine ("fftwf_plan_with_nthreads" plan-with-nthreads
)
30 (plan-with-nthreads 4)))
32 ;; to clean up completely call void fftw_cleanup_threads(void)
35 (define-alien-routine ("fftwf_execute" execute
)
39 (defconstant +forward
+ 1)
40 (defconstant +backward
+ -
1)
41 (defconstant +estimate
+ (ash 1 6))
44 (define-alien-routine ("fftwf_plan_dft_3d" plan-dft-3d
)
50 (out (* double-float
))
54 (define-alien-routine ("fftwf_plan_dft_2d" plan-dft-2d
)
59 (out (* double-float
))
63 (defun ft3-csf (in &key
(forward t
))
64 (declare ((simple-array (complex single-float
) 3) in
)
66 (values (simple-array (complex single-float
) 3) &optional
))
67 (let ((dims (array-dimensions in
)))
68 (destructuring-bind (z y x
)
70 (let* ((out (make-array dims
:element-type
'(complex single-float
))))
71 (sb-sys:with-pinned-objects
(in out
)
72 (let ((p (plan-dft-3d z y x
74 (sb-ext:array-storage-vector in
))
76 (sb-ext:array-storage-vector out
))
82 (when forward
;; normalize if forward
83 (let ((1/n
(/ 1s0
(* x y z
))))
84 (vol:do-box
(k j i
0 z
0 y
0 x
)
85 (setf (aref out k j i
) (* 1/n
(aref out k j i
))))))
88 (defmacro ift3-csf
(in)
89 `(ft3 ,in
:forward nil
))
91 (defun ft2-csf (in &key
(forward t
))
92 (declare ((simple-array (complex single-float
) 2) in
)
94 (values (simple-array (complex single-float
) 2) &optional
))
95 (let ((dims (array-dimensions in
)))
96 (destructuring-bind (y x
)
98 (let* ((out (make-array dims
:element-type
'(complex single-float
))))
99 (sb-sys:with-pinned-objects
(in out
)
100 (let ((p (plan-dft-2d y x
102 (sb-ext:array-storage-vector in
))
104 (sb-ext:array-storage-vector out
))
110 (when forward
;; normalize if forward
111 (let ((1/n
(/ 1s0
(* x y
))))
112 (vol:do-rectangle
(j i
0 y
0 x
)
113 (setf (aref out j i
) (* 1/n
(aref out j i
))))))
116 (defmacro ift2-csf
(in)
117 `(ft2 ,in
:forward nil
))
126 (a (vol:convert3-ub8
/csf-complex
127 (vol:draw-sphere-ub8
20d0 nz ny nx
))))
128 (vol:write-pgm
"/home/martin/tmp/fftwf.pgm"
129 (vol:normalize2-csf
/ub8-abs
130 (vol:cross-section-xz-csf
(ft3-csf a
)))))))