Commit 9d22c860 authored by Liam Healy's avatar Liam Healy
Browse files

Unify radix2 forward transforms over real and complex types

The individual real and complex defmfuns for forward transforms have
been removed from fft-real.lisp and fft-complex.lisp and unified in
fft-interface.lisp as forward-fourier-transform-radix2 and
forward-fourier-transform-nonradix2.  This is untested and is meant to
experiment with unification of the functions; if successful, the other
functions should be changed as similarly, and radix2/non-radix2 switch
incorporated if possible.
parent 0e3f8c2a
......@@ -8,18 +8,6 @@
;; Power of 2 functions
(defmfun fft-complex-radix2-forward
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_forward")
(((c-pointer vector) :pointer) (stride sizet) (n sizet))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"Forward FFT for a complex radix-2 vector")
(defmfun fft-complex-radix2-backward
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft" :type "_radix2_backward")
......@@ -221,20 +209,4 @@
:documentation
"Inverse FFT for a complex vector")
(defmfun fft-complex-transform
((vector vector) direction &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)))))
("gsl_fft" :type "_transform")
(((c-pointer vector) :pointer) (stride sizet) (n sizet)
((mpointer wavetable) :pointer) ((mpointer workspace) :pointer)
(direction fft-direction))
:definition :generic
:element-types :complex
:inputs (vector)
:outputs (vector)
:return (vector)
:documentation
"FFT in the given direction for a complex vector")
......@@ -5,7 +5,7 @@
;; Utility to determine if a vector can use a radix-2 transform.
(defun number-is-radix2 (num)
(= num (expt 2 (floor (log num 2)))))
(= 1 (logcount num)))
;; Generalised ways of making wavetables
(export 'make-fft-wavetable)
......@@ -61,6 +61,43 @@
;;; 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
......@@ -69,20 +106,13 @@
(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)))))))
(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)
......
......@@ -6,20 +6,6 @@
(in-package :gsl)
;; Power of 2 functions
(defmfun fft-real-radix2-transform
((vector vector) &key (stride 1) (n (expt 2 (floor (log (size vector) 2)))))
("gsl_fft_real" :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 real radix-2 vector")
;; Mixed Radix general-N functions
(defmobject fft-real-wavetable
......@@ -79,23 +65,6 @@
"This function allocates a workspace for a real float transform of length
n.")
(defmfun fft-real-transform
((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)))))
("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")
(defmfun fft-real-unpack
((vector vector)
&key (stride 1) (n (size vector))
......
;; Helpers that define a single GSL function interface
;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
;; Time-stamp: <2009-04-13 22:09:35EDT defmfun-single.lisp>
;; Time-stamp: <2009-11-01 14:20:53EST defmfun-single.lisp>
;; $Id: $
(in-package :gsl)
......@@ -76,24 +76,25 @@
(cl-argument-types arglist c-arguments)
(set-difference ; find all the unused variables
(arglist-plain-and-categories arglist nil)
(remove-duplicates (union
(if mapdown
(apply 'union
(mapcar 'variables-used-in-c-arguments c-arguments))
(variables-used-in-c-arguments c-arguments))
;; Forms in :before, :after are checked for used variables
(stupid-code-walk-find-variables
(cons
'values
(append before after
(callback-symbol-set
callback-dynamic cbinfo (first callback-dynamic-variables))
;; &optional/&key/&aux defaults are checked
(let ((auxstart (after-llk arglist)))
(when auxstart
(apply
'append
(mapcar 'rest (remove-if 'atom auxstart)))))))))))
(remove-duplicates
(union
(if mapdown
(apply 'union
(mapcar 'variables-used-in-c-arguments c-arguments))
(variables-used-in-c-arguments c-arguments))
;; Forms in :before, :after are checked for used variables
(stupid-code-walk-find-variables
(cons
'values
(append before after
(callback-symbol-set
callback-dynamic cbinfo (first callback-dynamic-variables))
;; &optional/&key/&aux defaults are checked
(let ((auxstart (after-llk arglist)))
(when auxstart
(apply
'append
(mapcar 'rest (remove-if 'atom auxstart)))))))))))
(first callback-dynamic-variables))
,@(when documentation (list documentation))
,(funcall body-maker name arglist gsl-name c-arguments key-args))
......
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