view Lib/fftw-3.2.1/genfft/gen_mdct.ml @ 1:e86e9c111b29

Updates stuff that potentially fixes the memory leak and also makes it work on Windows and Linux (Need to test). Still have to fix fftw include for linux in Jucer.

(*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 * Copyright (c) 2003, 2007-8 Matteo Frigo
 * Copyright (c) 2003, 2007-8 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *)

(* generation of trigonometric transforms *)

open Util
open Genutil
open C


let usage = "Usage: " ^ Sys.argv.(0) ^ " -n <number>"

let uistride = ref Stride_variable
let uostride = ref Stride_variable
let uivstride = ref Stride_variable
let uovstride = ref Stride_variable
let normalization = ref 1

type mode =
  | MDCT
  | MDCT_MP3
  | MDCT_VORBIS
  | MDCT_WINDOW
  | MDCT_WINDOW_SYM
  | IMDCT
  | IMDCT_MP3
  | IMDCT_VORBIS
  | IMDCT_WINDOW
  | IMDCT_WINDOW_SYM
  | NONE

let mode = ref NONE

let speclist = [
  "-with-istride",
  Arg.String(fun x -> uistride := arg_to_stride x),
  " specialize for given input stride";

  "-with-ostride",
  Arg.String(fun x -> uostride := arg_to_stride x),
  " specialize for given output stride";

  "-with-ivstride",
  Arg.String(fun x -> uivstride := arg_to_stride x),
  " specialize for given input vector stride";

  "-with-ovstride",
  Arg.String(fun x -> uovstride := arg_to_stride x),
  " specialize for given output vector stride";

  "-normalization",
  Arg.String(fun x -> normalization := int_of_string x),
  " normalization integer to divide by";

  "-mdct",
  Arg.Unit(fun () -> mode := MDCT),
  " generate an MDCT codelet";

  "-mdct-mp3",
  Arg.Unit(fun () -> mode := MDCT_MP3),
  " generate an MDCT codelet with MP3 windowing";

  "-mdct-window",
  Arg.Unit(fun () -> mode := MDCT_WINDOW),
  " generate an MDCT codelet with window array";

  "-mdct-window-sym",
  Arg.Unit(fun () -> mode := MDCT_WINDOW_SYM),
  " generate an MDCT codelet with symmetric window array";

  "-imdct",
  Arg.Unit(fun () -> mode := IMDCT),
  " generate an IMDCT codelet";

  "-imdct-mp3",
  Arg.Unit(fun () -> mode := IMDCT_MP3),
  " generate an IMDCT codelet with MP3 windowing";

  "-imdct-window",
  Arg.Unit(fun () -> mode := IMDCT_WINDOW),
  " generate an IMDCT codelet with window array";

  "-imdct-window-sym",
  Arg.Unit(fun () -> mode := IMDCT_WINDOW_SYM),
  " generate an IMDCT codelet with symmetric window array";
]

let unity_window n i = Complex.one

(* MP3 window(k) = sin(pi/(2n) * (k + 1/2)) *)
let mp3_window n k = 
  Complex.imag (Complex.exp (8 * n) (2*k + 1))

(* Vorbis window(k) = sin(pi/2 * (mp3_window(k))^2)
    ... this is transcendental, though, so we can't do it with our
        current Complex.exp function *)

let window_array n w =
    array n (fun i ->
      let stride = C.SInteger 1
      and klass = Unique.make () in
      let refr = C.array_subscript w stride i in
      let kr = Variable.make_constant klass refr in
      load_r (kr, kr))

let load_window w n i = w i
let load_window_sym w n i = w (if (i < n) then i else (2*n - 1 - i))

(* fixme: use same locations for input and output so that it works in-place? *)

(* Note: only correct for even n! *)
let load_array_mdct window n rarr iarr locations =
  let twon = 2 * n in
  let arr = load_array_c twon 
      (locative_array_c twon rarr iarr locations "BUG") in
  let arrw = fun i -> Complex.times (window n i) (arr i) in
  array n
    ((Complex.times Complex.half) @@
     (fun i ->
       if (i < n/2) then
	 Complex.uminus (Complex.plus [arrw (i + n + n/2); 
				       arrw (n + n/2 - 1 - i)])
       else
	 Complex.plus [arrw (i - n/2); 
		       Complex.uminus (arrw (n + n/2 - 1 - i))]))

let store_array_mdct window n rarr iarr locations arr =
  store_array_r n (locative_array_c n rarr iarr locations "BUG") arr

let load_array_imdct window n rarr iarr locations =
  load_array_c n (locative_array_c n rarr iarr locations "BUG")

let store_array_imdct window n rarr iarr locations arr =
  let n2 = n/2 in
  let threen2 = 3*n2 in
  let arr2 = fun i ->
    if (i < n2) then
      arr (i + n2)
    else if (i < threen2) then
      Complex.uminus (arr (threen2 - 1 - i))
    else
      Complex.uminus (arr (i - threen2))
  in
  let arr2w = fun i -> Complex.times (window n i) (arr2 i) in
  let twon = 2 * n in
  store_array_r twon (locative_array_c twon rarr iarr locations "BUG") arr2w

let window_param = function
    MDCT_WINDOW -> true
  | MDCT_WINDOW_SYM -> true
  | IMDCT_WINDOW -> true
  | IMDCT_WINDOW_SYM -> true
  | _ -> false

let generate n mode =
  let iarray = "I"
  and oarray = "O"
  and istride = "istride"
  and ostride = "ostride" 
  and window = "W" 
  and name = !Magic.codelet_name in

  let vistride = either_stride (!uistride) (C.SVar istride)
  and vostride = either_stride (!uostride) (C.SVar ostride)
  in

  let sivs = stride_to_string "ovs" !uovstride in
  let sovs = stride_to_string "ivs" !uivstride in

  let (transform, load_input, store_output) = match mode with
  | MDCT -> Trig.dctIV, load_array_mdct unity_window,
      store_array_mdct unity_window
  | MDCT_MP3 -> Trig.dctIV, load_array_mdct mp3_window,
      store_array_mdct unity_window
  | MDCT_WINDOW -> Trig.dctIV, load_array_mdct
	(load_window (window_array (2 * n) window)),
      store_array_mdct unity_window
  | MDCT_WINDOW_SYM -> Trig.dctIV, load_array_mdct
	(load_window_sym (window_array n window)),
      store_array_mdct unity_window
  | IMDCT -> Trig.dctIV, load_array_imdct unity_window,
      store_array_imdct unity_window
  | IMDCT_MP3 -> Trig.dctIV, load_array_imdct unity_window,
      store_array_imdct mp3_window
  | IMDCT_WINDOW -> Trig.dctIV, load_array_imdct unity_window,
      store_array_imdct (load_window (window_array (2 * n) window))
  | IMDCT_WINDOW_SYM -> Trig.dctIV, load_array_imdct unity_window,
      store_array_imdct (load_window_sym (window_array n window))
  | _ -> failwith "must specify transform kind"
  in
    
  let locations = unique_array_c (2*n) in
  let input = 
    load_input n
      (C.array_subscript iarray vistride)
      (C.array_subscript "BUG" vistride)
      locations
  in
  let output = (Complex.times (Complex.inverse_int !normalization)) 
    @@ (transform n input) in
  let odag =
    store_output n
      (C.array_subscript oarray vostride)
      (C.array_subscript "BUG" vostride)
      locations 
      output
  in
  let annot = standard_optimizer odag in

  let tree =
    Fcn ("void", name,
	 ([Decl (C.constrealtypep, iarray);
	   Decl (C.realtypep, oarray)]
	  @ (if stride_fixed !uistride then [] 
               else [Decl (C.stridetype, istride)])
	  @ (if stride_fixed !uostride then [] 
	       else [Decl (C.stridetype, ostride)])
	  @ (choose_simd []
	       (if stride_fixed !uivstride then [] else 
	       [Decl ("int", sivs)]))
	  @ (choose_simd []
	       (if stride_fixed !uovstride then [] else 
	       [Decl ("int", sovs)]))
	  @ (if (not (window_param mode)) then [] 
	       else [Decl (C.constrealtypep, window)])
	 ),
	 add_constants (Asch annot))

  in
  (unparse tree) ^ "\n"


let main () =
  begin
    parse speclist usage;
    print_string (generate (check_size ()) !mode);
  end

let _ = main()