fftw-fft.lisp

   1 (declaim (optimize (speed 2) (debug 3) (safety 3)))
   2 #.(require :vol)
   3 (defpackage :fftw-fft
   4   (:use :cl :sb-alien :sb-c-call)
   5   (:export #:ft3-cdf
   6                  #:ift3-cdf
   7                  #:ift2-cdf
   8                  #:ft2-cdf))
   9 (in-package :fftw-fft)
  10
  11 (load-shared-object "/usr/lib/libfftw3.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 didn't seem to get faster
  17 #+nil
  18 (progn
  19   (load-shared-object "/usr/lib/libfftw3_threads.so")
  20
  21   (define-alien-routine ("fftw_init_threads" init-threads)
  22       int)
  23   (define-alien-routine ("fftw_plan_with_nthreads" plan-with-nthreads)
  24       void
  25     (nthreads int))
  26
  27   (defun init-ft ()
  28     (init-threads)
  29     (plan-with-nthreads 4)))
  30
  31 ;; to clean up completely call void fftw_cleanup_threads(void)
  32
  33
  34 (define-alien-routine ("fftw_execute" execute)
  35     void
  36   (plan (* int)))
  37
  38 (defconstant +forward+ 1)
  39 (defconstant +backward+ -1)
  40 (defconstant +estimate+ (ash 1 6))
  41
  42
  43 (define-alien-routine ("fftw_plan_dft_3d" plan-dft-3d)
  44     (* int)
  45   (n0 int)
  46   (n1 int)
  47   (n2 int)
  48   (in (* double-float))
  49   (out (* double-float))
  50   (sign int)
  51   (flags unsigned-int))
  52
  53 (define-alien-routine ("fftw_plan_dft_2d" plan-dft-2d)
  54     (* int)
  55   (n0 int)
  56   (n1 int)
  57   (in (* double-float))
  58   (out (* double-float))
  59   (sign int)
  60   (flags unsigned-int))
  61
  62 (defun ft3-cdf (in &key (forward t))
  63   (declare ((simple-array (complex double-float) 3) in)
  64            (boolean forward)
  65            (values (simple-array (complex double-float) 3) &optional))
  66   (let ((dims (array-dimensions in)))
  67    (destructuring-bind (z y x)
  68        dims
  69      (let* ((out (make-array dims :element-type '(complex double-float))))
  70        (sb-sys:with-pinned-objects (in out)
  71          (let ((p (plan-dft-3d z y x
  72                                (sb-sys:vector-sap
  73                                 (sb-ext:array-storage-vector in))
  74                                (sb-sys:vector-sap
  75                                 (sb-ext:array-storage-vector out))
  76                                (if forward
  77                                    +forward+
  78                                    +backward+)
  79                                +estimate+)))
  80            (execute p)))
  81        (when forward ;; normalize if forward
  82          (let ((1/n (/ 1d0 (* x y z))))
  83            (vol:do-box (k j i 0 z 0 y 0 x)
  84              (setf (aref out k j i) (* 1/n (aref out k j i))))))
  85        out))))
  86
  87 (defmacro ift3-cdf (in)
  88   `(ft3 ,in :forward nil))
  89
  90 (defun ft2-cdf (in &key (forward t))
  91   (declare ((simple-array (complex double-float) 2) in)
  92            (boolean forward)
  93            (values (simple-array (complex double-float) 2) &optional))
  94   (let ((dims (array-dimensions in)))
  95     (destructuring-bind (y x)
  96         dims
  97       (let* ((out (make-array dims :element-type '(complex double-float))))
  98         (sb-sys:with-pinned-objects (in out)
  99           (let ((p (plan-dft-2d y x
 100                                 (sb-sys:vector-sap
 101                                  (sb-ext:array-storage-vector in))
 102                                 (sb-sys:vector-sap
 103                                  (sb-ext:array-storage-vector out))
 104                                 (if forward
 105                                     +forward+
 106                                    +backward+)
 107                                 +estimate+)))
 108             (execute p)))
 109         (when forward ;; normalize if forward
 110           (let ((1/n (/ 1d0 (* x y))))
 111             (vol:do-rectangle (j i 0 y 0 x)
 112               (setf (aref out j i) (* 1/n (aref out j i))))))
 113         out))))
 114
 115 (defmacro ift2-cdf (in)
 116   `(ft2 ,in :forward nil))
 117
 118 #+nil
 119 (progn
 120   (time
 121    (let* ((nx 256)
 122           (ny nx)
 123           (nz ny)
 124           (a (vol:convert3-ub8/cdf-complex
 125               (vol:draw-sphere-ub8 20d0 nz ny nx))))
 126      (vol:write-pgm "/home/martin/tmp/fftwf.pgm"
 127                 (vol:normalize2-cdf/ub8-abs
 128                  (vol:cross-section-xz-cdf (ft3-cdf a)))))))