Skip to content
Snippets Groups Projects
Commit b4f09c28 authored by Liam Healy's avatar Liam Healy
Browse files

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.
parent 9d22c860
No related branches found
No related tags found
No related merge requests found
......@@ -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>
......
;; 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")
;; 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))
;; 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))
......
;; 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)))
......
;; 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)
......
;; 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
......
;; 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)))))
;; 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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment