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
|
(*fft.sml*)
(*by Torben Mogensen (torbenm@diku.dk)*)
val pi = 3.14159265358979
fun pr (s : string) : unit = print(s)
exception Impossible
fun impossible s = impossible0 (s ^ "\n")
and impossible0 s = (pr ("\nimpossible: " ^ s); raise Impossible)
fun zipWith f ([],[]) = []
| zipWith f ((a::b),(c::d)) = f (a,c) :: zipWith f (b,d)
| zipWith f _ = impossible "zipWith"
fun zip ([],[]) = []
| zip ((a::b),(c::d)) = (a,c) :: zip (b,d)
| zip _ = impossible "zip"
fun ~+ ((x:real,y:real),(v,w)) = (x+v,y+w)
fun ~- ((x:real,y:real),(v,w)) = (x-v,y-w)
fun ~* ((x:real,y:real),(v,w)) = (x*v-y*w,x*w+y*v)
fun evens [] = []
| evens (x::y::l) = x :: evens l
| evens _ = impossible "evens"
fun odds [] = []
| odds (x::y::l) = y :: odds l
| odds _ = impossible "odds"
fun fmul (c,pin,[]) = []
| fmul (c,pin,(a::b))
= ~*((Math.cos(c),Math.sin(c)), a) :: fmul (c+pin,pin,b)
fun cp [] = []
| cp (a::b) = a :: cp b
fun fft ([(a,b)], 1) = [(a+0.0,b+0.0)]
| fft (x, n2)
= let val n = n2 div 2
val a = fft (evens x, n)
val cb = fmul (0.0,pi/(real n),fft (odds x, n))
in
let val l1 = zipWith ~+ (a,cb)
val l2 = zipWith ~- (a,cb)
in (*resetRegions a; resetRegions cb;*) l1 @ l2
end
end
local val a = 16807.0 and m = 2147483678.0
in
fun nextrand seed =
let val t = a*seed
in t - m * real(floor (t/m)) end
end
fun mkList(tr as (seed,0,acc)) = tr
| mkList(seed,n,acc) = mkList(nextrand seed, n-1, seed::acc)
val n = 256 * 256
fun run () = (pr "\nfft by Torben Mogensen (torbenm@diku.dk)\n\nfft'ing... ";
let val r = fft (zip (#3(mkList(7.0,n,[])),
#3(mkList(8.0,n,[]))), n) in
pr " done\n" end);
val _ = run ()
|
