From dcc05d7fd50cf08deda3dbc389b024f9243f1e98 Mon Sep 17 00:00:00 2001 From: Sumant Oemrawsingh <soemraws@xs4all.nl> Date: Sun, 1 Nov 2009 04:28:57 +0100 Subject: [PATCH] Added attempt to clean FFT interface A hopefully nicer interface to the FFT functions are given, where the precise FFT function is deduced from the supplied vector. While the interface is nicer, the code is ugly, hacky and not fool proof yet. Cleanup required. --- .../fft-interface-example.lisp | 39 +++++ fast-fourier-transforms/fft-interface.lisp | 137 ++++++++++++++++++ gsll.asd | 1 + 3 files changed, 177 insertions(+) create mode 100644 fast-fourier-transforms/fft-interface-example.lisp create mode 100644 fast-fourier-transforms/fft-interface.lisp diff --git a/fast-fourier-transforms/fft-interface-example.lisp b/fast-fourier-transforms/fft-interface-example.lisp new file mode 100644 index 00000000..693c0ff8 --- /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 00000000..9a33f58a --- /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 cf377b1e..7fdd9185 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) -- GitLab