1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
;;;"fft.scm" Fast Fourier Transform
;Copyright (C) 1999, 2003 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.
;;;; 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)
;;; http://www.astro.virginia.edu/~eww6n/math/DiscreteFourierTransform.html
;;; differs in the direction of rotation of the complex unit vectors.
(require 'array)
(require 'logical)
;;@code{(require 'fft)}
;;@ftindex fft
(define (fft:shuffle&scale new ara n scale)
(define lgn (integer-length (+ -1 n)))
(if (not (eqv? n (expt 2 lgn)))
(slib:error 'fft "array length not power of 2" n))
(do ((k 0 (+ 1 k)))
((>= k n) new)
(array-set! new (* (array-ref ara k) scale) (reverse-bit-field k 0 lgn))))
(define (dft! ara n dir)
(define lgn (integer-length (+ -1 n)))
(define pi2i (* 0+8i (atan 1)))
(do ((s 1 (+ 1 s)))
((> s lgn) ara)
(let* ((m (expt 2 s))
(w_m (exp (* dir (/ 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 n))
(let* ((k+m/2 (+ k m/2-1 1))
(t (* w (array-ref ara k+m/2)))
(u (array-ref ara k)))
(array-set! ara (+ u t) k)
(array-set! ara (- u t) k+m/2)))))))
;;@args array
;;@var{array} is an array of @code{(expt 2 n)} numbers. @code{fft}
;;returns an array of complex numbers comprising the
;;@dfn{Discrete Fourier Transform} of @var{array}.
(define (fft ara)
(define n (car (array-dimensions ara)))
(define new (apply make-array ara (array-dimensions ara)))
(dft! (fft:shuffle&scale new ara n 1) n 1))
;;@args array
;;@code{fft-1} returns an array of complex numbers comprising the
;;inverse Discrete Fourier Transform of @var{array}.
(define (fft-1 ara)
(define n (car (array-dimensions ara)))
(define new (apply make-array ara (array-dimensions ara)))
(dft! (fft:shuffle&scale new ara n (/ n)) n -1))
;;@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
|
