summaryrefslogtreecommitdiffstats
path: root/dft.scm
diff options
context:
space:
mode:
Diffstat (limited to 'dft.scm')
-rw-r--r--dft.scm195
1 files changed, 195 insertions, 0 deletions
diff --git a/dft.scm b/dft.scm
new file mode 100644
index 0000000..29180d0
--- /dev/null
+++ b/dft.scm
@@ -0,0 +1,195 @@
+;;;"dft.scm" Discrete Fourier Transform
+;Copyright (C) 1999, 2003, 2006 Aubrey Jaffer
+;
+;Permission to copy this software, to modify it, to redistribute it,
+;to distribute modified versions, and to use it for any purpose is
+;granted, subject to the following restrictions and understandings.
+;
+;1. Any copy made of this software must include this copyright notice
+;in full.
+;
+;2. I have made no warranty or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3. In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+;;;; For one-dimensional power-of-two length see:
+;;; Introduction to Algorithms (MIT Electrical
+;;; Engineering and Computer Science Series)
+;;; by Thomas H. Cormen, Charles E. Leiserson (Contributor),
+;;; Ronald L. Rivest (Contributor)
+;;; MIT Press; ISBN: 0-262-03141-8 (July 1990)
+
+;;; Flipped polarity of exponent to agree with
+;;; http://en.wikipedia.org/wiki/Discrete_Fourier_transform
+
+(require 'array)
+(require 'logical)
+(require 'subarray)
+
+;;@code{(require 'dft)} or
+;;@code{(require 'Fourier-transform)}
+;;@ftindex dft, Fourier-transform
+;;
+;;@code{fft} and @code{fft-1} compute the Fast-Fourier-Transforms
+;;(O(n*log(n))) of arrays whose dimensions are all powers of 2.
+;;
+;;@code{sft} and @code{sft-1} compute the Discrete-Fourier-Transforms
+;;for all combinations of dimensions (O(n^2)).
+
+(define (dft:sft1d! new ara n dir)
+ (define scl (if (negative? dir) (/ 1.0 n) 1))
+ (define pi2i/n (/ (* 0-8i (atan 1) dir) n))
+ (do ((k (+ -1 n) (+ -1 k)))
+ ((negative? k) new)
+ (let ((sum 0))
+ (do ((j (+ -1 n) (+ -1 j)))
+ ((negative? j) (array-set! new sum k))
+ (set! sum (+ sum (* (exp (* pi2i/n j k))
+ (array-ref ara j)
+ scl)))))))
+
+(define (dft:fft1d! new ara n dir)
+ (define scl (if (negative? dir) (/ 1.0 n) 1))
+ (define lgn (integer-length (+ -1 n)))
+ (define pi2i (* 0-8i (atan 1) dir))
+ (do ((k 0 (+ 1 k)))
+ ((>= k n))
+ (array-set! new (* (array-ref ara k) scl) (reverse-bit-field k 0 lgn)))
+ (do ((s 1 (+ 1 s))
+ (m (expt 2 1) (expt 2 (+ 1 s))))
+ ((> s lgn) new)
+ (let ((w_m (exp (/ pi2i m)))
+ (m/2-1 (+ (quotient m 2) -1)))
+ (do ((j 0 (+ 1 j))
+ (w 1 (* w w_m)))
+ ((> j m/2-1))
+ (do ((k j (+ m k))
+ (k+m/2 (+ j m/2-1 1) (+ m k m/2-1 1)))
+ ((>= k n))
+ (let ((t (* w (array-ref new k+m/2)))
+ (u (array-ref new k)))
+ (array-set! new (+ u t) k)
+ (array-set! new (- u t) k+m/2)))))))
+
+;;; Row-major order is suboptimal for Scheme.
+;;; N are copied into and operated on in place
+;;; A[a, *, c] --> N1[c, a, *]
+;;; N1[c, *, b] --> N2[b, c, *]
+;;; N2[b, *, a] --> N3[a, b, *]
+
+(define (dft:rotate-indexes idxs)
+ (define ridxs (reverse idxs))
+ (cons (car ridxs) (reverse (cdr ridxs))))
+
+(define (dft:dft prot ara dir transform-1d)
+ (define (ranker ara rdx dims)
+ (define ndims (dft:rotate-indexes dims))
+ (if (negative? rdx)
+ ara
+ (let ((new (apply make-array prot ndims))
+ (rdxlen (car (last-pair ndims))))
+ (define x1d
+ (cond (transform-1d)
+ ((eqv? rdxlen (expt 2 (integer-length (+ -1 rdxlen))))
+ dft:fft1d!)
+ (else dft:sft1d!)))
+ (define (ramap rdims inds)
+ (cond ((null? rdims)
+ (x1d (apply subarray new (dft:rotate-indexes inds))
+ (apply subarray ara inds)
+ rdxlen dir))
+ ((null? inds)
+ (do ((i (+ -1 (car rdims)) (+ -1 i)))
+ ((negative? i))
+ (ramap (cddr rdims)
+ (cons #f (cons i inds)))))
+ (else
+ (do ((i (+ -1 (car rdims)) (+ -1 i)))
+ ((negative? i))
+ (ramap (cdr rdims) (cons i inds))))))
+ (if (= 1 (length dims))
+ (x1d new ara rdxlen dir)
+ (ramap (reverse dims) '()))
+ (ranker new (+ -1 rdx) ndims))))
+ (ranker ara (+ -1 (array-rank ara)) (array-dimensions ara)))
+
+;;@args array prot
+;;@args array
+;;@var{array} is an array of positive rank. @code{sft} returns an
+;;array of type @2 (defaulting to @1) of complex numbers comprising
+;;the @dfn{Discrete Fourier Transform} of @var{array}.
+(define (sft ara . prot)
+ (dft:dft (if (null? prot) ara (car prot)) ara 1 dft:sft1d!))
+
+;;@args array prot
+;;@args array
+;;@var{array} is an array of positive rank. @code{sft-1} returns an
+;;array of type @2 (defaulting to @1) of complex numbers comprising
+;;the inverse Discrete Fourier Transform of @var{array}.
+(define (sft-1 ara . prot)
+ (dft:dft (if (null? prot) ara (car prot)) ara -1 dft:sft1d!))
+
+(define (dft:check-dimensions ara name)
+ (for-each (lambda (n)
+ (if (not (eqv? n (expt 2 (integer-length (+ -1 n)))))
+ (slib:error name "array length not power of 2" n)))
+ (array-dimensions ara)))
+
+;;@args array prot
+;;@args array
+;;@var{array} is an array of positive rank whose dimensions are all
+;;powers of 2. @code{fft} returns an array of type @2 (defaulting to
+;;@1) of complex numbers comprising the Discrete Fourier Transform of
+;;@var{array}.
+(define (fft ara . prot)
+ (dft:check-dimensions ara 'fft)
+ (dft:dft (if (null? prot) ara (car prot)) ara 1 dft:fft1d!))
+
+;;@args array prot
+;;@args array
+;;@var{array} is an array of positive rank whose dimensions are all
+;;powers of 2. @code{fft-1} returns an array of type @2 (defaulting
+;;to @1) of complex numbers comprising the inverse Discrete Fourier
+;;Transform of @var{array}.
+(define (fft-1 ara . prot)
+ (dft:check-dimensions ara 'fft-1)
+ (dft:dft (if (null? prot) ara (car prot)) ara -1 dft:fft1d!))
+
+;;@code{dft} and @code{dft-1} compute the discrete Fourier transforms
+;;using the best method for decimating each dimension.
+
+;;@args array prot
+;;@args array
+;;@0 returns an array of type @2 (defaulting to @1) of complex
+;;numbers comprising the Discrete Fourier Transform of @var{array}.
+(define (dft ara . prot)
+ (dft:dft (if (null? prot) ara (car prot)) ara 1 #f))
+
+;;@args array prot
+;;@args array
+;;@0 returns an array of type @2 (defaulting to @1) of
+;;complex numbers comprising the inverse Discrete Fourier Transform of
+;;@var{array}.
+(define (dft-1 ara . prot)
+ (dft:dft (if (null? prot) ara (car prot)) ara -1 #f))
+
+;;@noindent
+;;@code{(fft-1 (fft @var{array}))} will return an array of values close to
+;;@var{array}.
+;;
+;;@example
+;;(fft '#(1 0+i -1 0-i 1 0+i -1 0-i)) @result{}
+;;
+;;#(0.0 0.0 0.0+628.0783185208527e-18i 0.0
+;; 0.0 0.0 8.0-628.0783185208527e-18i 0.0)
+;;
+;;(fft-1 '#(0 0 0 0 0 0 8 0)) @result{}
+;;
+;;#(1.0 -61.23031769111886e-18+1.0i -1.0 61.23031769111886e-18-1.0i
+;; 1.0 -61.23031769111886e-18+1.0i -1.0 61.23031769111886e-18-1.0i)
+;;@end example