| 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 () |