Commit | Line | Data |
---|---|---|
7f918cf1 CE |
1 | (*fft.sml*) |
2 | ||
3 | (*by Torben Mogensen (torbenm@diku.dk)*) | |
4 | ||
5 | ||
6 | val pi = 3.14159265358979 | |
7 | ||
8 | fun pr (s : string) : unit = print(s) | |
9 | exception Impossible | |
10 | fun impossible s = impossible0 (s ^ "\n") | |
11 | and impossible0 s = (pr ("\nimpossible: " ^ s); raise Impossible) | |
12 | ||
13 | fun zipWith f ([],[]) = [] | |
14 | | zipWith f ((a::b),(c::d)) = f (a,c) :: zipWith f (b,d) | |
15 | | zipWith f _ = impossible "zipWith" | |
16 | ||
17 | fun zip ([],[]) = [] | |
18 | | zip ((a::b),(c::d)) = (a,c) :: zip (b,d) | |
19 | | zip _ = impossible "zip" | |
20 | ||
21 | fun ~+ ((x:real,y:real),(v,w)) = (x+v,y+w) | |
22 | ||
23 | fun ~- ((x:real,y:real),(v,w)) = (x-v,y-w) | |
24 | ||
25 | fun ~* ((x:real,y:real),(v,w)) = (x*v-y*w,x*w+y*v) | |
26 | ||
27 | fun evens [] = [] | |
28 | | evens (x::y::l) = x :: evens l | |
29 | | evens _ = impossible "evens" | |
30 | ||
31 | fun odds [] = [] | |
32 | | odds (x::y::l) = y :: odds l | |
33 | | odds _ = impossible "odds" | |
34 | ||
35 | fun fmul (c,pin,[]) = [] | |
36 | | fmul (c,pin,(a::b)) | |
37 | = ~*((Math.cos(c),Math.sin(c)), a) :: fmul (c+pin,pin,b) | |
38 | ||
39 | fun cp [] = [] | |
40 | | cp (a::b) = a :: cp b | |
41 | ||
42 | fun fft ([(a,b)], 1) = [(a+0.0,b+0.0)] | |
43 | | fft (x, n2) | |
44 | = let val n = n2 div 2 | |
45 | val a = fft (evens x, n) | |
46 | val cb = fmul (0.0,pi/(real n),fft (odds x, n)) | |
47 | in | |
48 | let val l1 = zipWith ~+ (a,cb) | |
49 | val l2 = zipWith ~- (a,cb) | |
50 | in (*resetRegions a; resetRegions cb;*) l1 @ l2 | |
51 | end | |
52 | end | |
53 | ||
54 | local val a = 16807.0 and m = 2147483678.0 | |
55 | in | |
56 | fun nextrand seed = | |
57 | let val t = a*seed | |
58 | in t - m * real(floor (t/m)) end | |
59 | end | |
60 | ||
61 | ||
62 | fun mkList(tr as (seed,0,acc)) = tr | |
63 | | mkList(seed,n,acc) = mkList(nextrand seed, n-1, seed::acc) | |
64 | ||
65 | val n = 256 * 256 | |
66 | ||
67 | fun run () = (pr "\nfft by Torben Mogensen (torbenm@diku.dk)\n\nfft'ing... "; | |
68 | let val r = fft (zip (#3(mkList(7.0,n,[])), | |
69 | #3(mkList(8.0,n,[]))), n) in | |
70 | pr " done\n" end); | |
71 | ||
72 | val _ = run () |