Commit ce3177ad authored by Liam Healy's avatar Liam Healy
Browse files

Fast Fourier transform complete

The fast Fourier transform port now complete except for the tests.
The backward, inverse, and selectable direction transforms are now
unified in their respective files, a new file has been created for the
wavetable and workspace mobjects.  The output of the example functions
looks plausible.
parent 7192e225
;; Backward FFT
;; Sumant Oemrawsingh, Sat Oct 24 2009 - 12:55
;; Time-stamp: <2009-11-07 10:46:20EST backward.lisp>
(in-package :gsl)
;; /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
;;;;****************************************************************************
;;;; Backward transformation
;;;;****************************************************************************
;;; Backward transformations where the size is a power of 2
(defmfun backward-fourier-transform-radix2
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_backward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:export nil
:index backward-fourier-transform
:documentation
"Backward FFT on a vector for which (floor length stride) is a power of 2.")
(defmfun backward-fourier-transform-nonradix2
((vector vector) &key (stride 1) (n (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))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:export nil
:index backward-fourier-transform
:documentation
"Backward FFT on a complex vector for which (floor length stride) is
not a power of 2.")
;;;;****************************************************************************
;;;; Half-complex
;;;;****************************************************************************
(defmfun backward-fourier-transform-halfcomplex-radix2
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft_halfcomplex" :type "_radix2_backward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (vector)
:outputs (vector)
:return (vector)
:export nil
:index backward-fourier-transform
:documentation
"Backward FFT on a vector for which (floor length stride) is a power of 2,
in half complex form.")
(defmfun backward-fourier-transform-halfcomplex-nonradix2
((vector vector) &key (stride 1) (n (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))
:definition :generic
:element-types :float
:inputs (vector)
:outputs (vector)
:return (vector)
:export nil
:index backward-fourier-transform
:documentation
"Backward FFT on a vector for which (floor length stride) is not a power of 2,
in half complex form.")
;;;;****************************************************************************
;;;; Decimation-in-frequency backward FFT
;;;;****************************************************************************
(defmfun backward-fourier-transform-dif-radix2
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_dif_backward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:export nil
:index backward-fourier-transform
:documentation
"Backward decimation-in-frequency FFT on a vector for which (floor
length stride) is a power of 2.")
;;;;****************************************************************************
;;;; Unified function backward-fourier-transform
;;;;****************************************************************************
;;; This selects among the 10 backward FFT functions that GSL defines.
(export 'backward-fourier-transform)
(defun backward-fourier-transform
(vector &rest args
&key decimation-in-frequency (stride 1) &allow-other-keys)
"Perform a backward 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 vector is real, it is assumed to be
in half-complex form."
(let ((pass-on-args (copy-list args)))
(remf pass-on-args :half-complex)
(remf pass-on-args :decimation-in-frequency)
(if (power-of-2-p (floor (size vector) stride))
(if (subtypep (element-type vector) 'real)
(apply 'backward-fourier-transform-halfcomplex-radix2
vector pass-on-args)
(if decimation-in-frequency
(apply 'backward-fourier-transform-dif-radix2 vector pass-on-args)
(apply 'backward-fourier-transform-radix2 vector pass-on-args)))
(if (subtypep (element-type vector) 'real)
(apply 'backward-fourier-transform-halfcomplex-nonradix2
vector pass-on-args)
(apply 'backward-fourier-transform-nonradix2 vector pass-on-args)))))
;; Functions for fast fourier transforms on complex data.
;; Sumant Oemrawsingh, Sat Oct 24 2009 - 12:55
;; Time-stamp: <2009-11-03 22:21:07EST fft-complex.lisp>
;; /usr/include/gsl/gsl_fft_complex.h
;; /usr/include/gsl/gsl_fft_complex_float.h
(in-package :gsl)
;; Power of 2 functions
(defmfun fft-complex-radix2-backward
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_backward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Backward FFT for a complex radix-2 vector")
(defmfun fft-complex-radix2-inverse
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_inverse")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Inverse FFT for a complex radix-2 vector")
(defmfun fft-complex-radix2-transform
((vector vector) direction &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_transform")
(((c-pointer vector) :pointer) (stride sizet) (n sizet) (direction fft-direction))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"FFT in the given direction for a complex radix-2 vector")
;; These are decimation-in-frequency versions of the radix-2 FFT functions.
;; I don't know what they are and how they are supposed to work, so I don't
;; know if they work.
(defmfun fft-complex-radix2-dif-backward
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_dif_backward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Decimation-in-frequency version of the backward FFT for a complex radix-2 vector")
(defmfun fft-complex-radix2-dif-inverse
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_inverse")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Decimation-in-frequency version of the inverse FFT for a complex radix-2 vector")
(defmfun fft-complex-radix2-dif-transform
((vector vector) direction &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_dif_transform")
(((c-pointer vector) :pointer) (stride sizet) (n sizet) (direction fft-direction))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Decimation-in-frequency version of the FFT in the given direction for a
complex radix-2 vector")
;; Mixed Radix general-N functions
(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"
:documentation
"This function prepares a trigonometric lookup table for a complex FFT of
length n. The function returns a pointer to the newly allocated
gsl_fft_complex_wavetable if no errors were detected, and a null pointer in
the case of error. The length n is factorized into a product of
subtransforms, and the factors and their trigonometric coefficients are
stored in the wavetable. The trigonometric coefficients are computed using
direct calls to sin and cos, for accuracy. Recursion relations could be used
to compute the lookup table faster, but if an application performs many FFTs
of the same length then this computation is a one-off overhead which does
not affect the final throughput.
The wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT
functions. The same wavetable can be used for both forward and backward (or
inverse) transforms of a given length.")
(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"
:documentation
"This function prepares a trigonometric lookup table for a complex float FFT
of length n. The function returns a pointer to the newly allocated
gsl_fft_complex_wavetable if no errors were detected, and a null pointer in
the case of error. The length n is factorized into a product of
subtransforms, and the factors and their trigonometric coefficients are
stored in the wavetable. The trigonometric coefficients are computed using
direct calls to sin and cos, for accuracy. Recursion relations could be used
to compute the lookup table faster, but if an application performs many FFTs
of the same length then this computation is a one-off overhead which does
not affect the final throughput.
The wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT
functions. The same wavetable can be used for both forward and backward (or
inverse) transforms of a given length.")
(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-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"
:documentation
"This function allocates a workspace for a complex float transform of length
n.")
(defmfun fft-complex-backward
((vector vector) &key (stride 1) (n (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))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Backward FFT for a complex vector")
(defmfun fft-complex-inverse
((vector vector) &key (stride 1) (n (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))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Inverse FFT for a complex vector")
;; Functions for fast fourier transforms on real data.
;; Sumant Oemrawsingh, Sat Oct 31 2009 - 20:12
;; Time-stamp: <2009-11-03 23:02:17EST fft-half-complex.lisp>
;; /usr/include/gsl/gsl_fft_halfcomplex.h
;; /usr/include/gsl/gsl_fft_halfcomplex_float.h
(in-package :gsl)
;; Power of 2 functions
(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")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Backward FFT for a half-complex radix-2 vector")
(defmfun fft-half-complex-radix2-inverse
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft_halfcomplex" :type "_radix2_inverse")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Inverse FFT for a half-complex radix-2 vector")
;; Mixed Radix general-N functions
(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"
:documentation
"These functions prepare trigonometric lookup tables for an FFT of size n
real elements. The functions return a pointer to the newly allocated struct
if no errors were detected, and a null pointer in the case of error. The
length n is factorized into a product of subtransforms, and the factors and
their trigonometric coefficients are stored in the wavetable. The
trigonometric coefficients are computed using direct calls to sin and cos,
for accuracy. Recursion relations could be used to compute the lookup table
faster, but if an application performs many FFTs of the same length then
computing the wavetable is a one-off overhead which does not affect the
final throughput.
The wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT
functions. The appropriate type of wavetable must be used for forward real
or inverse half-complex transforms.")
(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"
:documentation
"These functions prepare trigonometric lookup tables for an FFT of size n
real float elements. The functions return a pointer to the newly allocated
struct if no errors were detected, and a null pointer in the case of error.
The length n is factorized into a product of subtransforms, and the factors
and their trigonometric coefficients are stored in the wavetable. The
trigonometric coefficients are computed using direct calls to sin and cos,
for accuracy. Recursion relations could be used to compute the lookup table
faster, but if an application performs many FFTs of the same length then
computing the wavetable is a one-off overhead which does not affect the
final throughput.
The wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT
functions. The appropriate type of wavetable must be used for forward real
or inverse half-complex transforms.")
(defmfun fft-half-complex-backward
((vector vector) &key (stride 1) (n (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))
:definition :generic
:element-types :float
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Forward FFT for a half-complex vector")
(defmfun fft-half-complex-inverse
((vector vector) &key (stride 1) (n (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))
:definition :generic
:element-types :float
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Forward FFT for a half-complex vector")
;; 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)
;; 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-single-float dimension)
(make-fft-real-wavetable-single-float dimension)))
((eql element-type 'double-float)
(if half-complex
(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-single-float dimension))
((equal element-type '(complex double-float))
(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-single-float dimension))
((eql element-type 'double-float)
(make-fft-real-workspace-double-float dimension))
((equal element-type '(complex single-float))
(make-fft-complex-workspace-single-float dimension))
((equal element-type '(complex double-float))
(make-fft-complex-workspace-double-float 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 '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 (power-of-2-p 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 (power-of-2-p 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)))))))
;; Functions for fast fourier transforms on real data.
;; Sumant Oemrawsingh, Sun Oct 25 2009 - 16:35
;; Time-stamp: <2009-11-03 23:02:17EST fft-real.lisp>
;; /usr/include/gsl/gsl_fft_real.h
;; /usr/include/gsl/gsl_fft_real_float.h
(in-package :gsl)
;; Mixed Radix general-N functions
(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"
:documentation
"These functions prepare trigonometric lookup tables for an FFT of size n
real elements. The functions return a pointer to the newly allocated struct
if no errors were detected, and a null pointer in the case of error. The
length n is factorized into a product of subtransforms, and the factors and
their trigonometric coefficients are stored in the wavetable. The
trigonometric coefficients are computed using direct calls to sin and cos,
for accuracy. Recursion relations could be used to compute the lookup table
faster, but if an application performs many FFTs of the same length then
computing the wavetable is a one-off overhead which does not affect the
final throughput.
The wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT
functions. The appropriate type of wavetable must be used for forward real
or inverse half-complex transforms.")
(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"
:documentation
"These functions prepare trigonometric lookup tables for an FFT of size n
real float elements. The functions return a pointer to the newly allocated
struct if no errors were detected, and a null pointer in the case of error.
The length n is factorized into a product of subtransforms, and the factors
and their trigonometric coefficients are stored in the wavetable. The