Commit 5d2b771c authored by Liam Healy's avatar Liam Healy
Browse files

Merge branch 'master' of ssh://repo.or.cz/srv/git/gsll

parents f1a11505 53fcec99
......@@ -126,16 +126,18 @@
(export 'backward-fourier-transform)
(defun backward-fourier-transform
(vector &rest args
&key decimation-in-frequency (stride 1) &allow-other-keys)
&key decimation-in-frequency (stride 1) non-radix-2 &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."
keyword arguments. If the vector is real, it is assumed to be in
half-complex form. If the length of the vector is a power of 2, use
of a non-radix-2 transform can be forced."
(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))
(remf pass-on-args :non-radix-2)
(if (and (not non-radix-2) (power-of-2-p (floor (size vector) stride)))
(if (subtypep (element-type vector) 'real)
(apply 'backward-fourier-transform-halfcomplex-radix2
vector pass-on-args)
......
......@@ -146,32 +146,33 @@
(setf (grid:gref vector i) (coerce (/ (grid:gref vector i) length) element-type)))
vector))
(defun test-real-fft-noise (vector &key (stride 1))
(defun test-real-fft-noise (vector &key (stride 1) non-radix-2)
"Test forward and inverse FFT for a real vector, and return both results in unpacked form."
(let* ((forward
(forward-fourier-transform (realpart-vector vector :stride stride :init-offset 0)
:stride stride))
:stride stride :non-radix-2 non-radix-2))
(inverse
(forward-fourier-transform (copy forward) :half-complex t :stride stride)))
(forward-fourier-transform (copy forward) :half-complex t :stride stride
:non-radix-2 non-radix-2)))
(values (unpack forward :unpack-type 'complex :stride stride)
(unpack (vector/length inverse :stride stride) :stride stride))))
(defun test-complex-fft-noise (vector &key (stride 1))
(defun test-complex-fft-noise (vector &key (stride 1) non-radix-2)
"Test forward, inverse and backward FFT for a complex vector and return all three results."
(let ((forward
(forward-fourier-transform
(copy-with-stride vector :stride stride :init-offset 2000)
:stride stride)))
:stride stride :non-radix-2 non-radix-2)))
(values forward
(inverse-fourier-transform (copy-with-stride forward :init-offset 0)
:stride stride)
(inverse-fourier-transform (copy-with-stride forward :stride stride :init-offset 0)
:stride stride :non-radix-2 non-radix-2)
;; in backward-fourier-transform, we could use copy instead of ;;
;; copy-with-stride
(backward-fourier-transform (copy-with-stride forward :init-offset
2000)
:stride stride))))
:stride stride :non-radix-2 non-radix-2))))
(defun test-fft-noise (element-type size &key (stride 1))
(defun test-fft-noise (element-type size &key (stride 1) non-radix-2)
"A test of real forward and complex forward, reverse, and inverse FFT
on random noise. Returns the result of the DFT forward Fourier transformation
and the forward FFT, which should be the same, and the original vector and the
......@@ -188,7 +189,8 @@
:init-offset 3000))))
(if (subtypep element-type 'complex)
(multiple-value-bind (forward inverse backward)
(test-complex-fft-noise random-vector :stride stride)
(test-complex-fft-noise random-vector :stride stride
:non-radix-2 non-radix-2)
(values
dft-random-vector ; DFT forward result for reference
forward ; FFT forward result; should check
......@@ -196,7 +198,8 @@
inverse ; The inverse FFT applied to the forward result
(vector/length backward :stride stride)))
(multiple-value-bind (forward inverse)
(test-real-fft-noise random-vector :stride stride)
(test-real-fft-noise random-vector :stride stride
:non-radix-2 non-radix-2)
(values dft-random-vector forward random-vector inverse)))))
;; (test-fft-noise '(complex double-float) 2 :stride 2)
......
......@@ -136,16 +136,19 @@
(export 'forward-fourier-transform)
(defun forward-fourier-transform
(vector &rest args
&key half-complex decimation-in-frequency (stride 1) &allow-other-keys)
&key half-complex decimation-in-frequency (stride 1) non-radix-2 &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."
then the key argument :half-complex should be non-NIL. If the
length of the vector is a power of 2, use of a non-radix-2 transform
can be forced."
(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))
(remf pass-on-args :non-radix-2)
(if (and (not non-radix-2) (power-of-2-p (floor (size vector) stride)))
(if half-complex
(apply 'forward-fourier-transform-halfcomplex-radix2
vector pass-on-args)
......
......@@ -126,16 +126,18 @@
(export 'inverse-fourier-transform)
(defun inverse-fourier-transform
(vector &rest args
&key decimation-in-frequency (stride 1) &allow-other-keys)
&key decimation-in-frequency (stride 1) non-radix-2 &allow-other-keys)
"Perform a inverse 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."
in half-complex form. If the length of the vector is a power of 2,
use of a non-radix-2 transform can be forced."
(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))
(remf pass-on-args :non-radix-2)
(if (and (not non-radix-2) (power-of-2-p (floor (size vector) stride)))
(if (subtypep (element-type vector) 'real)
(apply 'inverse-fourier-transform-halfcomplex-radix2
vector pass-on-args)
......
......@@ -22,6 +22,7 @@
(defvar *allowed-ticks* 1000000)
#+nil
(defun off-stride (array stride)
"Make a grid that is all of the original array except that the
stride points are excluded."
......@@ -32,7 +33,19 @@
unless (zerop (mod i stride))
collect (grid:gref* array i))))
(defmacro fft-complex-result-check (form element-component-type)
(defun fft-complex-off-stride-check (vector stride offset)
(loop with size = (size vector)
for i from 0 below size
unless (zerop (mod i stride))
unless (lisp-unit:numerical-equal
(complex (+ offset (* 2 i))
(+ offset (* 2 i) 1))
(grid:gref vector i))
return nil
finally
(return t)))
(defmacro fft-complex-result-check (form element-component-type stride)
"T if all FFT tests pass."
`(multiple-value-bind
(dft fft original fft-roundtrip fft-reverse)
......@@ -43,15 +56,29 @@
'+sgl-epsilon+
'+dbl-epsilon+)))
(lisp-unit:*measure* :infinity))
(when (> ,stride 1)
(lisp-unit:assert-true (fft-complex-off-stride-check fft ,stride 2000))
(lisp-unit:assert-true (fft-complex-off-stride-check fft-roundtrip ,stride 0))
(lisp-unit:assert-true (fft-complex-off-stride-check fft-reverse ,stride 2000))
(setf dft (grid:stride dft ,stride)
fft (grid:stride fft ,stride)
original (grid:stride original ,stride)
fft-roundtrip (grid:stride fft-roundtrip ,stride)
fft-reverse (grid:stride fft-reverse ,stride)))
(lisp-unit:assert-norm-equal dft fft)
(lisp-unit:assert-norm-equal original fft-roundtrip)
(lisp-unit:assert-norm-equal original fft-reverse))))
(defmacro fft-real-result-check (form element-type)
(defmacro fft-real-result-check (form element-type stride)
"T if all FFT tests pass."
`(multiple-value-bind
(dft fft original fft-roundtrip)
(progn (reset-urand) ,form)
(when (> ,stride 1)
(setf dft (grid:stride dft ,stride)
fft (grid:stride fft ,stride)
original (grid:stride original ,stride)
fft-roundtrip (grid:stride fft-roundtrip ,stride)))
(let ((lisp-unit:*epsilon*
(* *allowed-ticks*
,(if (eq element-type 'single-float)
......@@ -64,15 +91,21 @@
(eval-when (:compile-toplevel :load-toplevel :execute)
(defun fft-test-forms (size stride)
`((fft-real-result-check
(test-fft-noise 'double-float ,size :stride ,stride) double-float)
(test-fft-noise 'double-float ,size :stride ,stride :non-radix-2 t) double-float ,stride)
(fft-real-result-check
(test-fft-noise 'single-float ,size :stride ,stride) single-float)
(test-fft-noise 'single-float ,size :stride ,stride :non-radix-2 t) single-float ,stride)
(fft-complex-result-check
(test-fft-noise
'(complex double-float) ,size :stride ,stride) double-float)
'(complex double-float) ,size :stride ,stride :non-radix-2 t) double-float ,stride)
(fft-complex-result-check
(test-fft-noise
'(complex single-float) ,size :stride ,stride) single-float))))
'(complex single-float) ,size :stride ,stride :non-radix-2 t) single-float ,stride)
(when (power-of-2-p (floor ,size ,stride))
(fft-real-result-check
(test-fft-noise 'double-float ,size :stride ,stride :non-radix-2 nil) double-float ,stride)
(fft-complex-result-check
(test-fft-noise
'(complex double-float) ,size :stride ,stride :non-radix-2 nil) double-float ,stride)))))
(defmacro all-fft-test-forms (size-max stride-max &optional additional-single-stride)
`(lisp-unit:define-test fast-fourier-transform
......@@ -86,5 +119,5 @@
(fft-test-forms size 1)))))
;;; (all-fft-test-forms 9 3 (64 99))
;(all-fft-test-forms 20 3)
;;; Tests commented out because they come out not so good.
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