Commit 48f7e66a authored by Sumant Oemrawsingh's avatar Sumant Oemrawsingh
Browse files

Some clean-up of real/half-complex FFT

Half-complex FFT functions were split off into their own file.
Wavetable and workspace names now have an fft- prepended to their name, for
clarity.
parent a4e52ac3
......@@ -111,7 +111,7 @@
;; Mixed Radix general-N functions
(defmobject complex-wavetable
(defmobject fft-complex-wavetable
"gsl_fft_complex_wavetable" ((n sizet))
"structure that holds the factorization and trigonometric lookup tables for
the mixed radix complex fft algorithm"
......@@ -132,7 +132,7 @@
functions. The same wavetable can be used for both forward and backward (or
inverse) transforms of a given length.")
(defmobject complex-wavetable-float
(defmobject fft-complex-wavetable-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"
......@@ -154,14 +154,14 @@
inverse) transforms of a given length.")
(defmobject complex-workspace
(defmobject fft-complex-workspace
"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 complex-workspace-float
(defmobject fft-complex-workspace-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"
......@@ -172,10 +172,10 @@
(defmfun fft-complex-forward
((vector vector) &key (stride 1) (n (size vector))
(wavetable (eltcase (complex single-float) (make-complex-wavetable-float (size vector))
(complex double-float) (make-complex-wavetable (size vector))))
(workspace (eltcase (complex single-float) (make-complex-workspace-float (size vector))
(complex double-float) (make-complex-workspace (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))
......@@ -189,10 +189,10 @@
(defmfun fft-complex-backward
((vector vector) &key (stride 1) (n (size vector))
(wavetable (eltcase (complex single-float) (make-complex-wavetable-float (size vector))
(complex double-float) (make-complex-wavetable (size vector))))
(workspace (eltcase (complex single-float) (make-complex-workspace-float (size vector))
(complex double-float) (make-complex-workspace (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 "_backward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet)
((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
......@@ -206,10 +206,10 @@
(defmfun fft-complex-inverse
((vector vector) &key (stride 1) (n (size vector))
(wavetable (eltcase (complex single-float) (make-complex-wavetable-float (size vector))
(complex double-float) (make-complex-wavetable (size vector))))
(workspace (eltcase (complex single-float) (make-complex-workspace-float (size vector))
(complex double-float) (make-complex-workspace (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 "_inverse")
(((c-pointer vector) :pointer) (stride sizet) (n sizet)
((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
......
;; Functions for fast fourier transforms on real data.
;; Sumant Oemrawsingh, Sat Oct 31 2009 - 20:12
;; /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-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")
(((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")
(defmfun fft-half-complex-radix2-unpack
((h-complex vector)
&key (stride 1) (n (expt 2 (floor (log (size h-complex) 2))))
(complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n)
t (make-marray '(complex double-float) :dimensions n))))
("gsl_fft_halfcomplex" :type "_radix2_unpack")
(((c-pointer h-complex) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (h-complex complex-output)
:outputs (complex-output)
:return (complex-output)
:documentation
"Convert an array of half-complex coefficients as returned by
real-fft-radix2-transform, into an ordinary complex array.")
;; Mixed Radix general-N functions
(defmobject fft-half-complex-wavetable
"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-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-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)))))
("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-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)))))
("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")
(defmfun fft-half-complex-unpack
((h-complex vector)
&key (stride 1) (n (size h-complex))
(complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n)
t (make-marray '(complex double-float) :dimensions n))))
("gsl_fft_halfcomplex" :type "_unpack")
(((c-pointer h-complex) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (h-complex complex-output)
:outputs (complex-output)
:return (complex-output)
:documentation
"This function converts an array of half-complex
coefficients as returned by fft-real-transform, into an ordinary complex
array. It fills in the complex array using the symmetry z_k = z_{n-k}^* to
reconstruct the redundant elements.")
......@@ -20,52 +20,9 @@
:documentation
"Forward FFT for a real radix-2 vector")
;; /usr/include/gsl/gsl_fft_halfcomplex.h
;; /usr/include/gsl/gsl_fft_halfcomplex_float.h
(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")
(defmfun fft-half-complex-radix2-unpack
((h-complex vector)
&key (stride 1) (n (expt 2 (floor (log (size h-complex) 2))))
(complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n)
t (make-marray '(complex double-float) :dimensions n))))
("gsl_fft_halfcomplex" :type "_radix2_unpack")
(((c-pointer h-complex) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (h-complex complex-output)
:outputs (complex-output)
:return (complex-output)
:documentation
"Convert an array of half-complex coefficients as returned by
real-fft-radix2-transform, into an ordinary complex array.")
;; Mixed Radix general-N functions
(defmobject real-wavetable
(defmobject fft-real-wavetable
"gsl_fft_real_wavetable" ((n sizet))
"structure that holds the factorization and trigonometric lookup tables for
the mixed radix real fft algorithm"
......@@ -86,28 +43,7 @@
functions. The appropriate type of wavetable must be used for forward real
or inverse half-complex transforms.")
(defmobject half-complex-wavetable
"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 real-wavetable-float
(defmobject fft-real-wavetable-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"
......@@ -128,36 +64,14 @@
functions. The appropriate type of wavetable must be used for forward real
or inverse half-complex transforms.")
(defmobject half-complex-wavetable-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.")
(defmobject real-workspace
(defmobject fft-real-workspace
"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 real-workspace-float
(defmobject fft-real-workspace-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"
......@@ -167,10 +81,10 @@
(defmfun fft-real-transform
((vector vector) &key (stride 1) (n (size vector))
(wavetable (eltcase single-float (make-real-wavetable-float (size vector))
double-float (make-real-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-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)))))
("gsl_fft_real" :type "_transform")
(((c-pointer vector) :pointer) (stride sizet) (n sizet)
((mpointer wavetable) :pointer) ((mpointer workspace) :pointer))
......@@ -180,42 +94,7 @@
:outputs (vector)
:return (vector)
:documentation
"FFT in the given direction for a real vector")
(defmfun fft-half-complex-transform
((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)))))
("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
"FFT in the given direction for a real vector")
(defmfun fft-half-complex-unpack
((h-complex vector)
&key (stride 1) (n (size h-complex))
(complex-output (eltcase single-float (make-marray '(complex single-float) :dimensions n)
t (make-marray '(complex double-float) :dimensions n))))
("gsl_fft_halfcomplex" :type "_unpack")
(((c-pointer h-complex) :pointer) ((c-pointer complex-output) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :float
:inputs (h-complex complex-output)
:outputs (complex-output)
:return (complex-output)
:documentation
"This function converts an array of half-complex
coefficients as returned by fft-real-transform, into an ordinary complex
array. It fills in the complex array using the symmetry z_k = z_{n-k}^* to
reconstruct the redundant elements.")
"Forward FFT for a real vector")
(defmfun fft-real-unpack
((vector vector)
......@@ -233,4 +112,3 @@
"This function converts a single real array into an equivalent complex
array (with imaginary part set to zero), suitable for fft-complex
routines.")
......@@ -138,6 +138,7 @@
:components
((:file "fft-complex")
(:file "fft-real")
(:file "fft-half-complex" :depends-on ("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