Import Upstream version 3a4upstream/3a4

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