regression/fft.sml

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