fftwf-fft.lisp

   1 (declaim (optimize (speed 2) (debug 3) (safety 3)))
   2 #.(require :vol)
   3 (defpackage :fftwf-fft
   4   (:use :cl :sb-alien :sb-c-call)
   5   (:export #:ft3-csf
   6            #:ift3-csf
   7            #:ift2-csf
   8            #:ft2-csf))
   9 (in-package :fftwf-fft)
  10
  11 (load-shared-object "/usr/lib/libfftw3f.so")
  12
  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
  17
  18 #+nil
  19 (progn
  20   (load-shared-object "/usr/lib/libfftw3f_threads.so")
  21
  22   (define-alien-routine ("fftwf_init_threads" init-threads)
  23       int)
  24   (define-alien-routine ("fftwf_plan_with_nthreads" plan-with-nthreads)
  25       void
  26     (nthreads int))
  27
  28   (defun init-ft ()
  29     (init-threads)
  30     (plan-with-nthreads 4)))
  31
  32 ;; to clean up completely call void fftw_cleanup_threads(void)
  33
  34
  35 (define-alien-routine ("fftwf_execute" execute)
  36     void
  37   (plan (* int)))
  38
  39 (defconstant +forward+ 1)
  40 (defconstant +backward+ -1)
  41 (defconstant +estimate+ (ash 1 6))
  42
  43
  44 (define-alien-routine ("fftwf_plan_dft_3d" plan-dft-3d)
  45     (* int)
  46   (n0 int)
  47   (n1 int)
  48   (n2 int)
  49   (in (* double-float))
  50   (out (* double-float))
  51   (sign int)
  52   (flags unsigned-int))
  53
  54 (define-alien-routine ("fftwf_plan_dft_2d" plan-dft-2d)
  55     (* int)
  56   (n0 int)
  57   (n1 int)
  58   (in (* double-float))
  59   (out (* double-float))
  60   (sign int)
  61   (flags unsigned-int))
  62
  63 (defun ft3-csf (in &key (forward t))
  64   (declare ((simple-array (complex single-float) 3) in)
  65            (boolean forward)
  66            (values (simple-array (complex single-float) 3) &optional))
  67   (let ((dims (array-dimensions in)))
  68    (destructuring-bind (z y x)
  69        dims
  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
  73                                (sb-sys:vector-sap
  74                                 (sb-ext:array-storage-vector in))
  75                                (sb-sys:vector-sap
  76                                 (sb-ext:array-storage-vector out))
  77                                (if forward
  78                                    +forward+
  79                                    +backward+)
  80                                +estimate+)))
  81            (execute p)))
  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))))))
  86        out))))
  87
  88 (defmacro ift3-csf (in)
  89   `(ft3 ,in :forward nil))
  90
  91 (defun ft2-csf (in &key (forward t))
  92   (declare ((simple-array (complex single-float) 2) in)
  93            (boolean forward)
  94            (values (simple-array (complex single-float) 2) &optional))
  95   (let ((dims (array-dimensions in)))
  96     (destructuring-bind (y x)
  97         dims
  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
 101                                 (sb-sys:vector-sap
 102                                  (sb-ext:array-storage-vector in))
 103                                 (sb-sys:vector-sap
 104                                  (sb-ext:array-storage-vector out))
 105                                 (if forward
 106                                     +forward+
 107                                    +backward+)
 108                                 +estimate+)))
 109             (execute p)))
 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))))))
 114         out))))
 115
 116 (defmacro ift2-csf (in)
 117   `(ft2 ,in :forward nil))
 118
 119
 120 #+nil
 121 (progn
 122   (time
 123    (let* ((nx 256)
 124           (ny nx)
 125           (nz ny)
 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)))))))