Commit dcc05d7f authored by Sumant Oemrawsingh's avatar Sumant Oemrawsingh
Browse files

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.
parent 48f7e66a
;; 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))
;; 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)))))))
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment