Commit e090c4d8 authored by Raymond Toy's avatar Raymond Toy
Browse files

Add test

The test also picked up an error in convert-inverse-rfft where it
wasn't handling FFTs of even length.
parent fd5147d7
......@@ -65,7 +65,7 @@
(declare (type (simple-array (complex single-float) (*)) x))
(let ((res (make-array n :element-type 'single-float)))
(setf (aref res 0) (realpart (aref x 0)))
(loop for j from 1 below (length x)
(loop for j from 1 below (if (evenp n) (1- (length x)) (length x))
for k from 1 by 2
do
(let ((z (aref x j)))
......@@ -124,3 +124,61 @@
;; allocated the correct size of the arrays.
(error "rfft1b failed with code ~A" ier))
inv)))
(defun test-rfft (n &key verbose)
"Test the rfft and inverse-rfft routines using a simple ramp of
length n"
;; For testing, use a simple ramp: 1, 2, 3, ..., N.
;; With some care, we can derive the FFT for this is
;;
;; X[0] = (N + 1) / 2
;; X[n] = -1/2 - cot(%pi*n/N)/2, n = 1, N
(let* ((x (make-array n :element-type 'single-float)))
(loop for k from 0 below n
do
(setf (aref x k) (float (+ k 1) 1f0)))
(let* ((xfrm (rfft x))
(xfrm-len (length xfrm))
(expected (make-array xfrm-len :element-type '(complex single-float)))
(noise-pwr 0.0)
(signal-pwr 0.0))
(declare (single-float noise-pwr signal-pwr))
(setf (aref expected 0) (complex (* 0.5 (+ n 1)) 0.0))
(loop for k from 1 below xfrm-len
with omega = (coerce (/ pi n) 'single-float)
do
(setf (aref expected k)
(complex -0.5 (/ -0.5 (tan (* omega k))))))
(when verbose
(format t "Forward transform; actual vs expected~%")
(loop for k from 0 below xfrm-len
do
(format t "~4d: ~A ~A~%" k (aref xfrm k) (aref expected k))))
(incf noise-pwr (expt (abs (- (aref xfrm 0) (* 0.5 (+ n 1)))) 2))
(incf signal-pwr (expt (* 0.5 (+ n 1)) 2))
(loop for k from 1 below (length xfrm)
do
(incf noise-pwr (expt (abs (- (aref xfrm k)
(aref expected k)))
2))
(incf signal-pwr (expt (abs (aref expected k)) 2)))
(let ((inv (inverse-rfft xfrm n))
(inv-noise-pwr 0.0)
(inv-signal-pwr 0.0))
(when verbose
(format t "Inverse transform; actual vs expected~%")
(loop for k from 0 below n
do
(format t "~4d: ~A ~A~%" k (aref inv k) (+ k 1))))
(loop for k from 0 below n
do
(incf inv-noise-pwr (expt (- (aref inv k) (+ k 1)) 2))
(incf inv-signal-pwr (expt (coerce (+ k 1) 'single-float) 2)))
(flet ((db (s n)
(if (zerop n)
1000.0
(* 10 (log (/ s n) 10)))))
(values (db signal-pwr noise-pwr) (db inv-signal-pwr inv-noise-pwr) noise-pwr signal-pwr inv-noise-pwr inv-signal-pwr))))))
\ No newline at end of file
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