diff --git a/fast-fourier-transforms/fft-interface-example.lisp b/fast-fourier-transforms/fft-interface-example.lisp new file mode 100644 index 0000000000000000000000000000000000000000..693c0ff8289e252fc3c53cfc84b36407fdba1099 --- /dev/null +++ b/fast-fourier-transforms/fft-interface-example.lisp @@ -0,0 +1,39 @@ +;; Example FFT: transform a pulse (using the "clean" fft interface) +;; Sumant Oemrawsingh, Sat Oct 31 2009 - 00:24 + +;;; 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 fft-pulse-test (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)) + (forward-fourier-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-interface.lisp b/fast-fourier-transforms/fft-interface.lisp new file mode 100644 index 0000000000000000000000000000000000000000..9a33f58a72b25142eda408a16dda7ed29d3df2e6 --- /dev/null +++ b/fast-fourier-transforms/fft-interface.lisp @@ -0,0 +1,137 @@ +;; A dirty hack to attempt to make a cleaner interface to the FFTs. +;; Sumant Oemrawsingh, Sat Oct 31 2009 - 23:48 + +(in-package :gsl) + +;; Utility to determine if a vector can use a radix-2 transform. +(defun number-is-radix2 (num) + (= num (expt 2 (floor (log num 2))))) + +;; Generalised ways of making wavetables +(export 'make-fft-wavetable) +(defun make-fft-wavetable (element-type dimension &optional (half-complex nil)) + "Make a wavetable for an FFT of the given element type and length. T can be + given as an optional third argument if the wavetable is meant for a Fourier + 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))) + ((eql element-type 'double-float) + (if half-complex + (make-fft-real-wavetable dimension) + (make-fft-half-complex-wavetable dimension))) + ((equal element-type '(complex single-float)) + (make-fft-complex-wavetable-float dimension)) + ((equal element-type '(complex double-float)) + (make-fft-complex-wavetable 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)) + ((eql element-type 'double-float) + (make-fft-real-workspace dimension)) + ((equal element-type '(complex single-float)) + (make-fft-complex-workspace-float dimension)) + ((equal element-type '(complex double-float)) + (make-fft-complex-workspace dimension)))) + + +;; An environment to allow more efficient FFTs of the same type and length +(export 'with-fourier-transform-environment) +(defmacro with-fourier-transform-environment + ((wavetable workspace element-type dimension &optional (half-complex nil)) + &body body) + "Create an environment where all FFTs will be performed on vectors of the + same type and with the same length. This allows to calculculate and reuse the + wavetable and workspace only once. + + The first and second arguments will be bound to the wavetable and workspace, + the third argument is the element type of the vectors to be FFT'd and the + fourth argument indicates the length of the vectors to which FFTs will be + applied. Optionally, T can be given as a fifth argument if the element type + of the vectors is real, but must be considered as half-complex." + `(let ((,wavetable (make-fft-wavetable ,element-type ,dimension ,half-complex)) + (,workspace (make-fft-workspace ,element-type ,dimension))) + ,@body)) + + +;;; General FFT functions + +(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) + (cond ((subtypep el-type 'float) + (fft-real-radix2-transform vector)) + ((subtypep el-type 'complex) + (fft-complex-radix2-forward vector))) + (progn (unless wavetable + (setf wavetable (make-fft-wavetable el-type len nil))) + (unless workspace + (setf workspace (make-fft-workspace el-type len))) + (cond ((subtypep el-type 'float) + (fft-real-transform vector :wavetable wavetable + :workspace workspace)) + ((subtypep el-type 'complex) + (fft-complex-forward 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 + 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. + + Note that if the vector is of type real, it is assumed that it contains + half-complex data." + (let ((el-type (element-type vector)) + (len (size vector))) + (if (number-is-radix2 len) + (cond ((subtypep el-type 'float) + (fft-half-complex-radix2-backward vector)) + ((subtypep el-type 'complex) + (fft-complex-radix2-backward vector))) + (progn (unless wavetable + (setf wavetable (make-fft-wavetable el-type len t))) + (unless workspace + (setf workspace (make-fft-workspace el-type len))) + (cond ((subtypep el-type 'float) + (fft-half-complex-backward vector :wavetable wavetable + :workspace workspace)) + ((subtypep el-type 'complex) + (fft-complex-backward vector :wavetable wavetable + :workspace workspace))))))) + +(export 'inverse-fourier-transform) +(defun inverse-fourier-transform (vector &key workspace wavetable) + "Generalised inverse FFT. Perform an inverse 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. + + Note that if the vector is of type real, it is assumed that it contains + half-complex data." + (let ((el-type (element-type vector)) + (len (size vector))) + (if (number-is-radix2 len) + (cond ((subtypep el-type 'float) + (fft-half-complex-radix2-inverse vector)) + ((subtypep el-type 'complex) + (fft-complex-radix2-inverse vector))) + (progn (unless wavetable + (setf wavetable (make-fft-wavetable el-type len t))) + (unless workspace + (setf workspace (make-fft-workspace el-type len))) + (cond ((subtypep el-type 'float) + (fft-half-complex-inverse vector :wavetable wavetable + :workspace workspace)) + ((subtypep el-type 'complex) + (fft-complex-inverse vector :wavetable wavetable + :workspace workspace))))))) diff --git a/gsll.asd b/gsll.asd index cf377b1e389e04037e731c24d1014699c9e214fb..7fdd91852b2f048a11046899b6ce348af223aac8 100644 --- a/gsll.asd +++ b/gsll.asd @@ -139,6 +139,7 @@ ((:file "fft-complex") (:file "fft-real") (:file "fft-half-complex" :depends-on ("fft-real")) + (:file "fft-interface" :depends-on ("fft-complex" "fft-half-complex" "fft-real")) (:file "fft-example" :depends-on ("fft-complex" "fft-real")))) (:module random :depends-on (init data)