From b4f09c286dd0077238c1195bcf3b476ffa3dda82 Mon Sep 17 00:00:00 2001
From: Liam Healy <liam@thinkpad.local>
Date: Sun, 1 Nov 2009 23:05:30 -0500
Subject: [PATCH] Unify forward FFT into new file and function

The function #'forward-fourier-transform selects among the 12 GSL
forward FFT functions, depending on whether the length of the vector
is a power of 2, whether the :half-complex flag is set, and the type
of the vector.  This exported function, along with its supporting
internal functions, are in the file forward.lisp.  The dimension
argument for fft-pulse-test is now mandatory.  All mobjects have been
renamed so that they end with either -single-float or -double-float,
consistent with the naming scheme in the rest of GSLL.  Eliminated
(evidently) obsolete example file fft-example.lisp.
---
 documentation/index.html                      |   7 +-
 fast-fourier-transforms/fft-complex.lisp      |  48 ++-----
 fast-fourier-transforms/fft-example.lisp      |  52 -------
 fast-fourier-transforms/fft-half-complex.lisp |  54 ++-----
 .../fft-interface-example.lisp                |   7 +-
 fast-fourier-transforms/fft-interface.lisp    |  82 ++---------
 fast-fourier-transforms/fft-real.lisp         |  18 +--
 fast-fourier-transforms/forward.lisp          | 136 ++++++++++++++++++
 gsll.asd                                      |   5 +-
 9 files changed, 202 insertions(+), 207 deletions(-)
 delete mode 100644 fast-fourier-transforms/fft-example.lisp
 create mode 100644 fast-fourier-transforms/forward.lisp

diff --git a/documentation/index.html b/documentation/index.html
index f958e1aa..e1f895c9 100644
--- a/documentation/index.html
+++ b/documentation/index.html
@@ -384,6 +384,11 @@ eigen-gen eigen-genv
 monte-carlo-plain monte-carlo-miser monte-carlo-vegas
 ode-stepper ode-evolution standard-control y-control
 yp-control scaled-control
+fft-real-wavetable-double-float fft-real-wavetable-single-float
+fft-real-workspace-double-float fft-real-workspace-single-float
+fft-complex-wavetable-double-float fft-complex-wavetable-single-float
+fft-complex-workspace-double-float fft-complex-workspace-single-float
+fft-half-complex-wavetable-double-float fft-half-complex-wavetable-single-float
 </pre>
 <p>
 An instance may be created with a function whose name is "make-"
@@ -429,7 +434,7 @@ and arrays used internally or for function return.
 <!-- Created: Feb 25 2005 -->
 <!-- hhmts start -->
 <small>
-Time-stamp: <2009-10-28 22:49:53EDT index.html>
+Time-stamp: <2009-11-01 22:50:22EST index.html>
 </small>
 <!-- hhmts end -->
  </div>
diff --git a/fast-fourier-transforms/fft-complex.lisp b/fast-fourier-transforms/fft-complex.lisp
index 96eff060..135629be 100644
--- a/fast-fourier-transforms/fft-complex.lisp
+++ b/fast-fourier-transforms/fft-complex.lisp
@@ -1,5 +1,6 @@
 ;; Functions for fast fourier transforms on complex data.
 ;; Sumant Oemrawsingh, Sat Oct 24 2009 - 12:55
+;; Time-stamp: <2009-11-01 22:46:09EST fft-complex.lisp>
 
 ;; /usr/include/gsl/gsl_fft_complex.h
 ;; /usr/include/gsl/gsl_fft_complex_float.h
@@ -99,7 +100,7 @@
 
 ;; Mixed Radix general-N functions
 
-(defmobject fft-complex-wavetable
+(defmobject fft-complex-wavetable-double-float
     "gsl_fft_complex_wavetable" ((n sizet))
   "structure that holds the factorization and trigonometric lookup tables for
   the mixed radix complex fft algorithm"
@@ -120,7 +121,7 @@
   functions. The same wavetable can be used for both forward and backward (or
   inverse) transforms of a given length.")
 
-(defmobject fft-complex-wavetable-float
+(defmobject fft-complex-wavetable-single-float
     "gsl_fft_complex_wavetable_float" ((n sizet))
   "structure that holds the factorization and trigonometric lookup tables for
   the mixed radix complex float fft algorithm"
@@ -141,15 +142,14 @@
   functions. The same wavetable can be used for both forward and backward (or
   inverse) transforms of a given length.")
 
-
-(defmobject fft-complex-workspace
+(defmobject fft-complex-workspace-double-float
     "gsl_fft_complex_workspace" ((n sizet))
   "Structure that holds the additional working space required for the
   intermediate steps of the mixed radix complex fft algoritms"
   :documentation
   "This function allocates a workspace for a complex transform of length n.")
 
-(defmobject fft-complex-workspace-float
+(defmobject fft-complex-workspace-single-float
     "gsl_fft_complex_workspace_float" ((n sizet))
   "Structure that holds the additional working space required for the
   intermediate steps of the mixed radix complex float fft algoritms"
@@ -157,30 +157,14 @@
   "This function allocates a workspace for a complex float transform of length
   n.")
 
-
-(defmfun fft-complex-forward
-    ((vector vector) &key (stride 1) (n (size vector))
-                     (wavetable (eltcase (complex single-float) (make-fft-complex-wavetable-float (size vector))
-                                         (complex double-float) (make-fft-complex-wavetable (size vector))))
-                     (workspace (eltcase (complex single-float) (make-fft-complex-workspace-float (size vector))
-                                         (complex double-float) (make-fft-complex-workspace (size vector)))))
-  ("gsl_fft" :type "_forward")
-  (((c-pointer vector) :pointer) (stride sizet) (n sizet)
-   ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
-  :definition :generic
-  :element-types :complex
-  :inputs (vector)
-  :outputs (vector)
-  :return (vector)
-  :documentation
-  "Forward FFT for a complex vector")
-
 (defmfun fft-complex-backward
     ((vector vector) &key (stride 1) (n (size vector))
-                     (wavetable (eltcase (complex single-float) (make-fft-complex-wavetable-float (size vector))
-                                         (complex double-float) (make-fft-complex-wavetable (size vector))))
-                     (workspace (eltcase (complex single-float) (make-fft-complex-workspace-float (size vector))
-                                         (complex double-float) (make-fft-complex-workspace (size vector)))))
+     (wavetable
+      (eltcase (complex single-float) (make-fft-complex-wavetable-single-float (size vector))
+	       (complex double-float) (make-fft-complex-wavetable-double-float (size vector))))
+     (workspace
+      (eltcase (complex single-float) (make-fft-complex-workspace-single-float (size vector))
+	       (complex double-float) (make-fft-complex-workspace-double-float (size vector)))))
   ("gsl_fft" :type "_backward")
   (((c-pointer vector) :pointer) (stride sizet) (n sizet)
    ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
@@ -194,10 +178,10 @@
 
 (defmfun fft-complex-inverse
     ((vector vector) &key (stride 1) (n (size vector))
-                     (wavetable (eltcase (complex single-float) (make-fft-complex-wavetable-float (size vector))
-                                         (complex double-float) (make-fft-complex-wavetable (size vector))))
-                     (workspace (eltcase (complex single-float) (make-fft-complex-workspace-float (size vector))
-                                         (complex double-float) (make-fft-complex-workspace (size vector)))))
+                     (wavetable (eltcase (complex single-float) (make-fft-complex-wavetable-single-float (size vector))
+                                         (complex double-float) (make-fft-complex-wavetable-double-float (size vector))))
+                     (workspace (eltcase (complex single-float) (make-fft-complex-workspace-single-float (size vector))
+                                         (complex double-float) (make-fft-complex-workspace-double-float (size vector)))))
   ("gsl_fft" :type "_inverse")
   (((c-pointer vector) :pointer) (stride sizet) (n sizet)
    ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
@@ -208,5 +192,3 @@
   :return (vector)
   :documentation
   "Inverse FFT for a complex vector")
-
-
diff --git a/fast-fourier-transforms/fft-example.lisp b/fast-fourier-transforms/fft-example.lisp
deleted file mode 100644
index a9a489b3..00000000
--- a/fast-fourier-transforms/fft-example.lisp
+++ /dev/null
@@ -1,52 +0,0 @@
-;; Example FFT: transform a pulse
-;; Sumant Oemrawsingh, Sat Oct 31 2009 - 15:31
-
-;;; Here is an example program modelled after the example given in Section
-;;; 15.3 of the GSL Manual, which computes the FFT of a short pulse. To make
-;;; the resulting fourier transform real the pulse is defined for equal
-;;; positive and negative times (-10 ... 10), where the negative times wrap
-;;; around the end of the array.
-;;; 
-;;; The output array from the example in Section 15.3 of the GSL Manual can be
-;;; reproduced with:
-;;; (fft-pulse-test '(complex double-float) :dimension 128)
-;;;
-;;; This example program also yields the same output array as the example
-;;; program in Section 15.4 of the GSL Manual:
-;;; (fft-pulse-test '(complex double-float) :dimension 630)
-
-(in-package :gsl)
-
-(defun generate-pulse (element-type &key dimension)
-  (assert (and (integerp dimension) (> dimension 20)))
-  (let ((pulse (make-marray element-type :dimensions dimension))
-        (init-value (coerce 1 element-type)))
-    (setf (maref pulse 0) init-value)
-    (loop for i from 1 to 10
-          do (setf (maref pulse i) init-value
-                   (maref pulse (- dimension i)) init-value))
-    pulse))
-
-(defun fft-pulse-test (element-type &key dimension)
-  (let ((pulse (generate-pulse element-type :dimension dimension))
-        (is-radix2 (= dimension (expt 2 (floor (log dimension 2))))))
-    (cond
-      ((subtypep element-type 'complex)
-       (if is-radix2
-         (fft-complex-radix2-forward pulse)
-         (fft-complex-forward pulse)))
-      ((subtypep element-type 'float)
-       (if is-radix2
-         (fft-real-radix2-transform pulse)
-         (fft-real-transform pulse))))))
-
-(save-test
-  fft
-  (fft-pulse-test '(complex single-float) 128)
-  (fft-pulse-test '(complex single-float) 630)
-  (fft-pulse-test '(complex double-float) 128)
-  (fft-pulse-test '(complex double-float) 630)
-  (fft-pulse-test 'single-float 128)
-  (fft-pulse-test 'single-float 630)
-  (fft-pulse-test 'double-float 128)
-  (fft-pulse-test 'double-float 630))
diff --git a/fast-fourier-transforms/fft-half-complex.lisp b/fast-fourier-transforms/fft-half-complex.lisp
index 28a14e08..66a12ba4 100644
--- a/fast-fourier-transforms/fft-half-complex.lisp
+++ b/fast-fourier-transforms/fft-half-complex.lisp
@@ -1,5 +1,6 @@
 ;; Functions for fast fourier transforms on real data.
 ;; Sumant Oemrawsingh, Sat Oct 31 2009 - 20:12
+;; Time-stamp: <2009-11-01 23:04:46EST fft-half-complex.lisp>
 
 ;; /usr/include/gsl/gsl_fft_halfcomplex.h
 ;; /usr/include/gsl/gsl_fft_halfcomplex_float.h
@@ -8,18 +9,6 @@
 
 ;; Power of 2 functions
 
-(defmfun fft-half-complex-radix2-transform
-    ((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
-  ("gsl_fft_halfcomplex" :type "_radix2_transform")
-  (((c-pointer vector) :pointer) (stride sizet) (n sizet))
-  :definition :generic
-  :element-types :float
-  :inputs (vector)
-  :outputs (vector)
-  :return (vector)
-  :documentation
-  "Forward FFT for a half-complex radix-2 vector")
-
 (defmfun fft-half-complex-radix2-backward
     ((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
   ("gsl_fft_halfcomplex" :type "_radix2_backward")
@@ -62,7 +51,7 @@
 
 ;; Mixed Radix general-N functions
 
-(defmobject fft-half-complex-wavetable
+(defmobject fft-half-complex-wavetable-double-float
     "gsl_fft_halfcomplex_wavetable" ((n sizet))
   "structure that holds the factorization and trigonometric lookup tables for
   the mixed radix halfcomplex fft algorithm"
@@ -83,7 +72,7 @@
   functions. The appropriate type of wavetable must be used for forward real
   or inverse half-complex transforms.")
 
-(defmobject fft-half-complex-wavetable-float
+(defmobject fft-half-complex-wavetable-single-float
     "gsl_fft_halfcomplex_wavetable_float" ((n sizet))
   "structure that holds the factorization and trigonometric lookup tables for
   the mixed radix real float fft algorithm"
@@ -104,29 +93,14 @@
   functions. The appropriate type of wavetable must be used for forward real
   or inverse half-complex transforms.")
 
-(defmfun fft-half-complex-transform
-    ((vector vector) &key (stride 1) (n (size vector))
-                     (wavetable (eltcase single-float (make-fft-half-complex-wavetable-float (size vector))
-                                         double-float (make-fft-half-complex-wavetable (size vector))))
-                     (workspace (eltcase single-float (make-fft-real-workspace-float (size vector))
-                                         double-float (make-fft-real-workspace (size vector)))))
-  ("gsl_fft_halfcomplex" :type "_transform")
-  (((c-pointer vector) :pointer) (stride sizet) (n sizet)
-   ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
-  :definition :generic
-  :element-types :float
-  :inputs (vector)
-  :outputs (vector)
-  :return (vector)
-  :documentation
-  "Forward FFT for a half-complex vector")
-
 (defmfun fft-half-complex-backward
     ((vector vector) &key (stride 1) (n (size vector))
-                     (wavetable (eltcase single-float (make-half-complex-wavetable-float (size vector))
-                                         double-float (make-half-complex-wavetable (size vector))))
-                     (workspace (eltcase single-float (make-real-workspace-float (size vector))
-                                         double-float (make-real-workspace (size vector)))))
+     (wavetable
+      (eltcase single-float (make-half-complex-wavetable-single-float (size vector))
+	       double-float (make-half-complex-wavetable-double-float (size vector))))
+     (workspace
+      (eltcase single-float (make-real-workspace-single-float (size vector))
+	       double-float (make-real-workspace-double-float (size vector)))))
   ("gsl_fft_halfcomplex" :type "_backward")
   (((c-pointer vector) :pointer) (stride sizet) (n sizet)
    ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
@@ -140,10 +114,12 @@
 
 (defmfun fft-half-complex-inverse
     ((vector vector) &key (stride 1) (n (size vector))
-                     (wavetable (eltcase single-float (make-half-complex-wavetable-float (size vector))
-                                         double-float (make-half-complex-wavetable (size vector))))
-                     (workspace (eltcase single-float (make-real-workspace-float (size vector))
-                                         double-float (make-real-workspace (size vector)))))
+     (wavetable
+      (eltcase single-float (make-half-complex-wavetable-single-float (size vector))
+	       double-float (make-half-complex-wavetable-double-float (size vector))))
+     (workspace
+      (eltcase single-float (make-real-workspace-single-float (size vector))
+	       double-float (make-real-workspace-double-float (size vector)))))
   ("gsl_fft_halfcomplex" :type "_inverse")
   (((c-pointer vector) :pointer) (stride sizet) (n sizet)
    ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
diff --git a/fast-fourier-transforms/fft-interface-example.lisp b/fast-fourier-transforms/fft-interface-example.lisp
index 693c0ff8..ea5371b0 100644
--- a/fast-fourier-transforms/fft-interface-example.lisp
+++ b/fast-fourier-transforms/fft-interface-example.lisp
@@ -1,5 +1,6 @@
 ;; Example FFT: transform a pulse (using the "clean" fft interface)
 ;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24
+;; Time-stamp: <2009-11-01 23:04:58EST fft-interface-example.lisp>
 
 ;;; Here is an example program modelled after the example given in Section
 ;;; 15.3 of the GSL Manual, which computes the FFT of a short pulse. To make
@@ -9,15 +10,15 @@
 ;;; 
 ;;; The output array from the example in Section 15.3 of the GSL Manual can be
 ;;; reproduced with:
-;;; (fft-pulse-test '(complex double-float) :dimension 128)
+;;; (fft-pulse-test '(complex double-float) 128)
 ;;;
 ;;; This example program also yields the same output array as the example
 ;;; program in Section 15.4 of the GSL Manual:
-;;; (fft-pulse-test '(complex double-float) :dimension 630)
+;;; (fft-pulse-test '(complex double-float) 630)
 
 (in-package :gsl)
 
-(defun fft-pulse-test (element-type &key dimension)
+(defun fft-pulse-test (element-type dimension)
   (assert (and (integerp dimension) (> dimension 20)))
   (let ((pulse (make-marray element-type :dimensions dimension))
         (init-value (coerce 1 element-type)))
diff --git a/fast-fourier-transforms/fft-interface.lisp b/fast-fourier-transforms/fft-interface.lisp
index 12028430..1267a565 100644
--- a/fast-fourier-transforms/fft-interface.lisp
+++ b/fast-fourier-transforms/fft-interface.lisp
@@ -1,12 +1,9 @@
 ;; A dirty hack to attempt to make a cleaner interface to the FFTs.
 ;; Sumant Oemrawsingh, Sat Oct 31 2009 - 23:48
+;; Time-stamp: <2009-11-01 22:47:41EST fft-interface.lisp>
 
 (in-package :gsl)
 
-;; Utility to determine if a vector can use a radix-2 transform.
-(defun number-is-radix2 (num)
-  (= 1 (logcount num)))
-
 ;; Generalised ways of making wavetables
 (export 'make-fft-wavetable)
 (defun make-fft-wavetable (element-type dimension &optional (half-complex nil))
@@ -15,29 +12,29 @@
   transform on a half-complex vector."
   (cond ((eql element-type 'single-float)
          (if half-complex
-           (make-fft-half-complex-wavetable-float dimension)
-           (make-fft-real-wavetable-float dimension)))
+           (make-fft-half-complex-wavetable-single-float dimension)
+           (make-fft-real-wavetable-single-float dimension)))
         ((eql element-type 'double-float)
          (if half-complex
-           (make-fft-real-wavetable dimension)
-           (make-fft-half-complex-wavetable dimension)))
+           (make-fft-real-wavetable-double-float dimension)
+           (make-fft-half-complex-wavetable-double-float dimension)))
         ((equal element-type '(complex single-float))
-         (make-fft-complex-wavetable-float dimension))
+         (make-fft-complex-wavetable-single-float dimension))
         ((equal element-type '(complex double-float))
-         (make-fft-complex-wavetable dimension))))
+         (make-fft-complex-wavetable-double-float dimension))))
 
 ;; Generalised ways of making workspaces
 (export 'make-fft-workspace)
 (defun make-fft-workspace (element-type dimension)
   "Make a wavetable for an FFT of the given element type and length."
   (cond ((eql element-type 'single-float)
-         (make-fft-real-workspace-float dimension))
+         (make-fft-real-workspace-single-float dimension))
         ((eql element-type 'double-float)
-         (make-fft-real-workspace dimension))
+         (make-fft-real-workspace-double-float dimension))
         ((equal element-type '(complex single-float))
-         (make-fft-complex-workspace-float dimension))
+         (make-fft-complex-workspace-single-float dimension))
         ((equal element-type '(complex double-float))
-         (make-fft-complex-workspace dimension))))
+         (make-fft-complex-workspace-double-float dimension))))
 
 
 ;; An environment to allow more efficient FFTs of the same type and length
@@ -61,59 +58,6 @@
 
 ;;; General FFT functions
 
-;; Power of 2 functions
-(defmfun forward-fourier-transform-radix2
-    ((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
-  (double-float "gsl_fft_real_radix2_transform"
-   single-float "gsl_fft_real_float_radix2_transform"
-   complex-double-float "gsl_fft_complex_radix2_forward"
-   complex-single-float "gsl_fft_complex_float_radix2_forward")
-  (((c-pointer vector) :pointer) (stride sizet) (n sizet))
-  :definition :generic
-  :element-types :float-complex
-  :inputs (vector)
-  :outputs (vector)
-  :return (vector)
-  :documentation
-  "Forward FFT for a real radix-2 vector")
-
-(defmfun forward-fourier-transform-nonradix2
-    ((vector vector) &key (stride 1) (n (size vector))
-     (wavetable (eltcase single-float (make-fft-real-wavetable-float (size vector))
-			 double-float (make-fft-real-wavetable (size vector))))
-     (workspace (eltcase single-float (make-fft-real-workspace-float (size vector))
-			 double-float (make-fft-real-workspace (size vector)))))
-  (double-float "gsl_fft_real_transform"
-   single-float "gsl_fft_real_float_transform"
-   complex-double-float "gsl_fft_complex_forward"
-   complex-single-float "gsl_fft_complex_float_forward")
-  ("gsl_fft_real" :type "_transform")
-  (((c-pointer vector) :pointer) (stride sizet) (n sizet)
-   ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
-  :definition :generic
-  :element-types :float
-  :inputs (vector)
-  :outputs (vector)
-  :return (vector)
-  :documentation
-  "Forward FFT for a real vector")
-
-(export 'forward-fourier-transform)
-(defun forward-fourier-transform (vector &key wavetable workspace)
-  "Generalised forward FFT. Perform a forward FFT on the given vector. If the
-  length of the vector is not a power of 2, and the user has a suitable
-  wavetable and/or workspace, these can be supplied as keyword arguments."
-  (let ((el-type (element-type vector))
-        (len (size vector)))
-    (if (number-is-radix2 len)
-	(forward-fourier-transform-radix2 vector)
-	(progn (unless wavetable
-		 (setf wavetable (make-fft-wavetable el-type len nil)))
-	       (unless workspace
-		 (setf workspace (make-fft-workspace el-type len)))
-	       (forward-fourier-transform-radix2
-		vector :wavetable wavetable :workspace workspace)))))
-
 (export 'backward-fourier-transform)
 (defun backward-fourier-transform (vector &key wavetable workspace)
   "Generalised backward FFT. Perform a backward FFT on the given vector. If
@@ -124,7 +68,7 @@
   half-complex data."
   (let ((el-type (element-type vector))
         (len (size vector)))
-    (if (number-is-radix2 len)
+    (if (power-of-2-p len)
       (cond ((subtypep el-type 'float)
              (fft-half-complex-radix2-backward vector))
             ((subtypep el-type 'complex)
@@ -150,7 +94,7 @@
   half-complex data."
   (let ((el-type (element-type vector))
         (len (size vector)))
-    (if (number-is-radix2 len)
+    (if (power-of-2-p len)
       (cond ((subtypep el-type 'float)
              (fft-half-complex-radix2-inverse vector))
             ((subtypep el-type 'complex)
diff --git a/fast-fourier-transforms/fft-real.lisp b/fast-fourier-transforms/fft-real.lisp
index ca85e3aa..33875e15 100644
--- a/fast-fourier-transforms/fft-real.lisp
+++ b/fast-fourier-transforms/fft-real.lisp
@@ -1,5 +1,6 @@
 ;; Functions for fast fourier transforms on real data.
 ;; Sumant Oemrawsingh, Sun Oct 25 2009 - 16:35
+;; Time-stamp: <2009-11-01 22:41:14EST fft-real.lisp>
 
 ;; /usr/include/gsl/gsl_fft_real.h
 ;; /usr/include/gsl/gsl_fft_real_float.h
@@ -8,7 +9,7 @@
 
 ;; Mixed Radix general-N functions
 
-(defmobject fft-real-wavetable
+(defmobject fft-real-wavetable-double-float
     "gsl_fft_real_wavetable" ((n sizet))
   "structure that holds the factorization and trigonometric lookup tables for
   the mixed radix real fft algorithm"
@@ -29,7 +30,7 @@
   functions. The appropriate type of wavetable must be used for forward real
   or inverse half-complex transforms.")
 
-(defmobject fft-real-wavetable-float
+(defmobject fft-real-wavetable-single-float
     "gsl_fft_real_wavetable_float" ((n sizet))
   "structure that holds the factorization and trigonometric lookup tables for
   the mixed radix real float fft algorithm"
@@ -50,14 +51,14 @@
   functions. The appropriate type of wavetable must be used for forward real
   or inverse half-complex transforms.")
 
-(defmobject fft-real-workspace
+(defmobject fft-real-workspace-double-float
     "gsl_fft_real_workspace" ((n sizet))
   "Structure that holds the additional working space required for the
   intermediate steps of the mixed radix real fft algoritms"
   :documentation
   "This function allocates a workspace for a real transform of length n.")
 
-(defmobject fft-real-workspace-float
+(defmobject fft-real-workspace-single-float
     "gsl_fft_real_workspace_float" ((n sizet))
   "Structure that holds the additional working space required for the
   intermediate steps of the mixed radix real float fft algoritms"
@@ -66,10 +67,11 @@
   n.")
 
 (defmfun fft-real-unpack
-         ((vector vector)
-          &key (stride 1) (n (size vector))
-          (complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n)
-                                        t (make-marray '(complex double-float) :dimensions n))))
+    ((vector vector)
+     &key (stride 1) (n (size vector))
+     (complex-output
+      (eltcase single-float (make-marray '(complex single-float) :dimensions n)
+	       t (make-marray '(complex double-float) :dimensions n))))
   ("gsl_fft_real" :type "_unpack")
   (((c-pointer vector) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet))
   :definition :generic
diff --git a/fast-fourier-transforms/forward.lisp b/fast-fourier-transforms/forward.lisp
new file mode 100644
index 00000000..cc130824
--- /dev/null
+++ b/fast-fourier-transforms/forward.lisp
@@ -0,0 +1,136 @@
+;; FFT transforms
+;; Sumant Oemrawsingh, Sat Oct 31 2009 - 23:48
+;; Time-stamp: <2009-11-01 22:43:18EST forward.lisp>
+
+(in-package :gsl)
+
+;; /usr/include/gsl/gsl_fft_real.h
+;; /usr/include/gsl/gsl_fft_real_float.h
+;; /usr/include/gsl/gsl_fft_complex.h
+;; /usr/include/gsl/gsl_fft_complex_float.h
+;; /usr/include/gsl/gsl_fft_halfcomplex.h
+;; /usr/include/gsl/gsl_fft_halfcomplex_float.h
+
+;;;;****************************************************************************
+;;;; Forward transformation
+;;;;****************************************************************************
+
+;; Forward transformations where the size is a power of 2
+(defmfun forward-fourier-transform-radix2 ((vector vector) &key (stride 1))
+  (double-float "gsl_fft_real_radix2_transform"
+   single-float "gsl_fft_real_float_radix2_transform"
+   complex-double-float "gsl_fft_complex_radix2_forward"
+   complex-single-float "gsl_fft_complex_float_radix2_forward")
+  (((c-pointer vector) :pointer) (stride sizet) ((size vector) sizet))
+  :definition :generic
+  :element-types :float-complex
+  :inputs (vector)
+  :outputs (vector)
+  :return (vector)
+  :export nil
+  :index forward-fourier-transform
+  :documentation
+  "Forward FFT on a vector for which (floor length stride) is a power of 2.")
+
+(defmfun forward-fourier-transform-nonradix2
+    ((vector vector) &key (stride 1)
+     (wavetable
+      (eltcase single-float (make-fft-real-wavetable-single-float (size vector))
+	       double-float (make-fft-real-wavetable-double-float (size vector))
+	       (complex single-float)
+	       (make-fft-complex-wavetable-single-float (size vector))
+	       (complex double-float)
+	       (make-fft-complex-wavetable-double-float (size vector))))
+     (workspace
+      (eltcase single-float (make-fft-real-workspace-single-float (size vector))
+	       double-float (make-fft-real-workspace-double-float (size vector))
+	       (complex single-float)
+	       (make-fft-complex-workspace-single-float (size vector))
+	       (complex double-float)
+	       (make-fft-complex-workspace-double-float (size vector)))))
+  (double-float "gsl_fft_real_transform"
+		single-float "gsl_fft_real_float_transform"
+		complex-double-float "gsl_fft_complex_forward"
+		complex-single-float "gsl_fft_complex_float_forward")
+  (((c-pointer vector) :pointer) (stride sizet) ((size vector) sizet)
+   ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
+  :definition :generic
+  :element-types :float-complex
+  :inputs (vector)
+  :outputs (vector)
+  :return (vector)
+  :export nil
+  :index forward-fourier-transform
+  :documentation
+  "Forward FFT on a vector for which (floor length stride) is not a power of 2.")
+
+;;;;****************************************************************************
+;;;; Half-complex  
+;;;;****************************************************************************
+
+(defmfun forward-fourier-transform-halfcomplex-radix2
+    ((vector vector) &key (stride 1))
+  ("gsl_fft_halfcomplex" :type "_radix2_transform")
+  (((c-pointer vector) :pointer) (stride sizet) ((size vector) sizet))
+  :definition :generic
+  :element-types :float
+  :inputs (vector)
+  :outputs (vector)
+  :return (vector)
+  :export nil
+  :index forward-fourier-transform
+  :documentation
+  "Forward FFT on a vector for which (floor length stride) is a power of 2,
+   in half complex form.")
+
+(defmfun forward-fourier-transform-halfcomplex-nonradix2
+    ((vector vector) &key (stride 1)
+     (wavetable
+      (eltcase single-float (make-fft-real-wavetable-single-float (size vector))
+	       double-float (make-fft-real-wavetable-double-float (size vector))))
+     (workspace
+      (eltcase single-float (make-fft-real-workspace-single-float (size vector))
+	       double-float (make-fft-real-workspace-double-float (size vector)))))
+  ("gsl_fft_halfcomplex" :type "_transform")
+  (((c-pointer vector) :pointer) (stride sizet) ((size vector) sizet)
+   ((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
+  :definition :generic
+  :element-types :float
+  :inputs (vector)
+  :outputs (vector)
+  :return (vector)
+  :export nil
+  :index forward-fourier-transform
+  :documentation
+  "Forward FFT on a vector for which (floor length stride) is not a power of 2,
+   in half complex form.")
+
+;;;;****************************************************************************
+;;;; Unified function forward-fourier-transform
+;;;;****************************************************************************
+
+;; Utility to determine if a vector can use a radix-2 transform.
+(defun power-of-2-p (num)
+  "The integer is a power of 2."
+  (= 1 (logcount num)))
+
+;;; This selects among the 12 forward FFT functions that GSL defines.
+(export 'forward-fourier-transform)
+(defun forward-fourier-transform
+    (vector &rest args &key half-complex (stride 1) &allow-other-keys)
+  "Perform a forward fast Fourier transform on the given vector. If
+  the length of the vector is not a power of 2, and the user has a
+  suitable wavetable and/or workspace, these can be supplied as
+  keyword arguments.  If the (real) vector is in half-complex form,
+  then the key argument :half-complex should be non-NIL."
+  (let ((pass-on-args (copy-list args)))
+    (remf pass-on-args :half-complex)
+    (if (power-of-2-p (floor (size vector) stride))
+	(if half-complex
+	    (apply 'forward-fourier-transform-halfcomplex-radix2
+		   vector pass-on-args)
+	    (apply 'forward-fourier-transform-radix2 vector pass-on-args))
+	(if half-complex
+	    (apply 'forward-fourier-transform-halfcomplex-nonradix2
+		   vector pass-on-args)
+	    (apply 'forward-fourier-transform-nonradix2 vector pass-on-args)))))
diff --git a/gsll.asd b/gsll.asd
index 7fdd9185..30a7b8db 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2009-09-19 18:54:29EDT gsll.asd>
+;; Time-stamp: <2009-11-01 23:01:16EST gsll.asd>
 
 (when (asdf:find-system :fsbv nil)
   (pushnew :fsbv *features*))
@@ -139,8 +139,9 @@
             ((:file "fft-complex")
              (:file "fft-real")
              (:file "fft-half-complex" :depends-on ("fft-real"))
+             (:file "forward")
              (:file "fft-interface" :depends-on ("fft-complex" "fft-half-complex" "fft-real"))
-             (:file "fft-example" :depends-on ("fft-complex" "fft-real"))))
+	     (:file "fft-interface-example")))
    (:module random
 	    :depends-on (init data)
 	    :components
-- 
GitLab