cannam@167: (*
cannam@167:  * Copyright (c) 1997-1999 Massachusetts Institute of Technology
cannam@167:  * Copyright (c) 2003, 2007-14 Matteo Frigo
cannam@167:  * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
cannam@167:  *
cannam@167:  * This program is free software; you can redistribute it and/or modify
cannam@167:  * it under the terms of the GNU General Public License as published by
cannam@167:  * the Free Software Foundation; either version 2 of the License, or
cannam@167:  * (at your option) any later version.
cannam@167:  *
cannam@167:  * This program is distributed in the hope that it will be useful,
cannam@167:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@167:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@167:  * GNU General Public License for more details.
cannam@167:  *
cannam@167:  * You should have received a copy of the GNU General Public License
cannam@167:  * along with this program; if not, write to the Free Software
cannam@167:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@167:  *
cannam@167:  *)
cannam@167: 
cannam@167: (* generation of trigonometric transforms *)
cannam@167: 
cannam@167: open Util
cannam@167: open Genutil
cannam@167: open C
cannam@167: 
cannam@167: 
cannam@167: let usage = "Usage: " ^ Sys.argv.(0) ^ " -n <number>"
cannam@167: 
cannam@167: let uistride = ref Stride_variable
cannam@167: let uostride = ref Stride_variable
cannam@167: let uivstride = ref Stride_variable
cannam@167: let uovstride = ref Stride_variable
cannam@167: let normalization = ref 1
cannam@167: 
cannam@167: type mode =
cannam@167:   | MDCT
cannam@167:   | MDCT_MP3
cannam@167:   | MDCT_VORBIS
cannam@167:   | MDCT_WINDOW
cannam@167:   | MDCT_WINDOW_SYM
cannam@167:   | IMDCT
cannam@167:   | IMDCT_MP3
cannam@167:   | IMDCT_VORBIS
cannam@167:   | IMDCT_WINDOW
cannam@167:   | IMDCT_WINDOW_SYM
cannam@167:   | NONE
cannam@167: 
cannam@167: let mode = ref NONE
cannam@167: 
cannam@167: let speclist = [
cannam@167:   "-with-istride",
cannam@167:   Arg.String(fun x -> uistride := arg_to_stride x),
cannam@167:   " specialize for given input stride";
cannam@167: 
cannam@167:   "-with-ostride",
cannam@167:   Arg.String(fun x -> uostride := arg_to_stride x),
cannam@167:   " specialize for given output stride";
cannam@167: 
cannam@167:   "-with-ivstride",
cannam@167:   Arg.String(fun x -> uivstride := arg_to_stride x),
cannam@167:   " specialize for given input vector stride";
cannam@167: 
cannam@167:   "-with-ovstride",
cannam@167:   Arg.String(fun x -> uovstride := arg_to_stride x),
cannam@167:   " specialize for given output vector stride";
cannam@167: 
cannam@167:   "-normalization",
cannam@167:   Arg.String(fun x -> normalization := int_of_string x),
cannam@167:   " normalization integer to divide by";
cannam@167: 
cannam@167:   "-mdct",
cannam@167:   Arg.Unit(fun () -> mode := MDCT),
cannam@167:   " generate an MDCT codelet";
cannam@167: 
cannam@167:   "-mdct-mp3",
cannam@167:   Arg.Unit(fun () -> mode := MDCT_MP3),
cannam@167:   " generate an MDCT codelet with MP3 windowing";
cannam@167: 
cannam@167:   "-mdct-window",
cannam@167:   Arg.Unit(fun () -> mode := MDCT_WINDOW),
cannam@167:   " generate an MDCT codelet with window array";
cannam@167: 
cannam@167:   "-mdct-window-sym",
cannam@167:   Arg.Unit(fun () -> mode := MDCT_WINDOW_SYM),
cannam@167:   " generate an MDCT codelet with symmetric window array";
cannam@167: 
cannam@167:   "-imdct",
cannam@167:   Arg.Unit(fun () -> mode := IMDCT),
cannam@167:   " generate an IMDCT codelet";
cannam@167: 
cannam@167:   "-imdct-mp3",
cannam@167:   Arg.Unit(fun () -> mode := IMDCT_MP3),
cannam@167:   " generate an IMDCT codelet with MP3 windowing";
cannam@167: 
cannam@167:   "-imdct-window",
cannam@167:   Arg.Unit(fun () -> mode := IMDCT_WINDOW),
cannam@167:   " generate an IMDCT codelet with window array";
cannam@167: 
cannam@167:   "-imdct-window-sym",
cannam@167:   Arg.Unit(fun () -> mode := IMDCT_WINDOW_SYM),
cannam@167:   " generate an IMDCT codelet with symmetric window array";
cannam@167: ]
cannam@167: 
cannam@167: let unity_window n i = Complex.one
cannam@167: 
cannam@167: (* MP3 window(k) = sin(pi/(2n) * (k + 1/2)) *)
cannam@167: let mp3_window n k = 
cannam@167:   Complex.imag (Complex.exp (8 * n) (2*k + 1))
cannam@167: 
cannam@167: (* Vorbis window(k) = sin(pi/2 * (mp3_window(k))^2)
cannam@167:     ... this is transcendental, though, so we can't do it with our
cannam@167:         current Complex.exp function *)
cannam@167: 
cannam@167: let window_array n w =
cannam@167:     array n (fun i ->
cannam@167:       let stride = C.SInteger 1
cannam@167:       and klass = Unique.make () in
cannam@167:       let refr = C.array_subscript w stride i in
cannam@167:       let kr = Variable.make_constant klass refr in
cannam@167:       load_r (kr, kr))
cannam@167: 
cannam@167: let load_window w n i = w i
cannam@167: let load_window_sym w n i = w (if (i < n) then i else (2*n - 1 - i))
cannam@167: 
cannam@167: (* fixme: use same locations for input and output so that it works in-place? *)
cannam@167: 
cannam@167: (* Note: only correct for even n! *)
cannam@167: let load_array_mdct window n rarr iarr locations =
cannam@167:   let twon = 2 * n in
cannam@167:   let arr = load_array_c twon 
cannam@167:       (locative_array_c twon rarr iarr locations "BUG") in
cannam@167:   let arrw = fun i -> Complex.times (window n i) (arr i) in
cannam@167:   array n
cannam@167:     ((Complex.times Complex.half) @@
cannam@167:      (fun i ->
cannam@167:        if (i < n/2) then
cannam@167: 	 Complex.uminus (Complex.plus [arrw (i + n + n/2); 
cannam@167: 				       arrw (n + n/2 - 1 - i)])
cannam@167:        else
cannam@167: 	 Complex.plus [arrw (i - n/2); 
cannam@167: 		       Complex.uminus (arrw (n + n/2 - 1 - i))]))
cannam@167: 
cannam@167: let store_array_mdct window n rarr iarr locations arr =
cannam@167:   store_array_r n (locative_array_c n rarr iarr locations "BUG") arr
cannam@167: 
cannam@167: let load_array_imdct window n rarr iarr locations =
cannam@167:   load_array_c n (locative_array_c n rarr iarr locations "BUG")
cannam@167: 
cannam@167: let store_array_imdct window n rarr iarr locations arr =
cannam@167:   let n2 = n/2 in
cannam@167:   let threen2 = 3*n2 in
cannam@167:   let arr2 = fun i ->
cannam@167:     if (i < n2) then
cannam@167:       arr (i + n2)
cannam@167:     else if (i < threen2) then
cannam@167:       Complex.uminus (arr (threen2 - 1 - i))
cannam@167:     else
cannam@167:       Complex.uminus (arr (i - threen2))
cannam@167:   in
cannam@167:   let arr2w = fun i -> Complex.times (window n i) (arr2 i) in
cannam@167:   let twon = 2 * n in
cannam@167:   store_array_r twon (locative_array_c twon rarr iarr locations "BUG") arr2w
cannam@167: 
cannam@167: let window_param = function
cannam@167:     MDCT_WINDOW -> true
cannam@167:   | MDCT_WINDOW_SYM -> true
cannam@167:   | IMDCT_WINDOW -> true
cannam@167:   | IMDCT_WINDOW_SYM -> true
cannam@167:   | _ -> false
cannam@167: 
cannam@167: let generate n mode =
cannam@167:   let iarray = "I"
cannam@167:   and oarray = "O"
cannam@167:   and istride = "istride"
cannam@167:   and ostride = "ostride" 
cannam@167:   and window = "W" 
cannam@167:   and name = !Magic.codelet_name in
cannam@167: 
cannam@167:   let vistride = either_stride (!uistride) (C.SVar istride)
cannam@167:   and vostride = either_stride (!uostride) (C.SVar ostride)
cannam@167:   in
cannam@167: 
cannam@167:   let sivs = stride_to_string "ovs" !uovstride in
cannam@167:   let sovs = stride_to_string "ivs" !uivstride in
cannam@167: 
cannam@167:   let (transform, load_input, store_output) = match mode with
cannam@167:   | MDCT -> Trig.dctIV, load_array_mdct unity_window,
cannam@167:       store_array_mdct unity_window
cannam@167:   | MDCT_MP3 -> Trig.dctIV, load_array_mdct mp3_window,
cannam@167:       store_array_mdct unity_window
cannam@167:   | MDCT_WINDOW -> Trig.dctIV, load_array_mdct
cannam@167: 	(load_window (window_array (2 * n) window)),
cannam@167:       store_array_mdct unity_window
cannam@167:   | MDCT_WINDOW_SYM -> Trig.dctIV, load_array_mdct
cannam@167: 	(load_window_sym (window_array n window)),
cannam@167:       store_array_mdct unity_window
cannam@167:   | IMDCT -> Trig.dctIV, load_array_imdct unity_window,
cannam@167:       store_array_imdct unity_window
cannam@167:   | IMDCT_MP3 -> Trig.dctIV, load_array_imdct unity_window,
cannam@167:       store_array_imdct mp3_window
cannam@167:   | IMDCT_WINDOW -> Trig.dctIV, load_array_imdct unity_window,
cannam@167:       store_array_imdct (load_window (window_array (2 * n) window))
cannam@167:   | IMDCT_WINDOW_SYM -> Trig.dctIV, load_array_imdct unity_window,
cannam@167:       store_array_imdct (load_window_sym (window_array n window))
cannam@167:   | _ -> failwith "must specify transform kind"
cannam@167:   in
cannam@167:     
cannam@167:   let locations = unique_array_c (2*n) in
cannam@167:   let input = 
cannam@167:     load_input n
cannam@167:       (C.array_subscript iarray vistride)
cannam@167:       (C.array_subscript "BUG" vistride)
cannam@167:       locations
cannam@167:   in
cannam@167:   let output = (Complex.times (Complex.inverse_int !normalization)) 
cannam@167:     @@ (transform n input) in
cannam@167:   let odag =
cannam@167:     store_output n
cannam@167:       (C.array_subscript oarray vostride)
cannam@167:       (C.array_subscript "BUG" vostride)
cannam@167:       locations 
cannam@167:       output
cannam@167:   in
cannam@167:   let annot = standard_optimizer odag in
cannam@167: 
cannam@167:   let tree =
cannam@167:     Fcn ("void", name,
cannam@167: 	 ([Decl (C.constrealtypep, iarray);
cannam@167: 	   Decl (C.realtypep, oarray)]
cannam@167: 	  @ (if stride_fixed !uistride then [] 
cannam@167:                else [Decl (C.stridetype, istride)])
cannam@167: 	  @ (if stride_fixed !uostride then [] 
cannam@167: 	       else [Decl (C.stridetype, ostride)])
cannam@167: 	  @ (choose_simd []
cannam@167: 	       (if stride_fixed !uivstride then [] else 
cannam@167: 	       [Decl ("int", sivs)]))
cannam@167: 	  @ (choose_simd []
cannam@167: 	       (if stride_fixed !uovstride then [] else 
cannam@167: 	       [Decl ("int", sovs)]))
cannam@167: 	  @ (if (not (window_param mode)) then [] 
cannam@167: 	       else [Decl (C.constrealtypep, window)])
cannam@167: 	 ),
cannam@167: 	 finalize_fcn (Asch annot))
cannam@167: 
cannam@167:   in
cannam@167:   (unparse tree) ^ "\n"
cannam@167: 
cannam@167: 
cannam@167: let main () =
cannam@167:   begin
cannam@167:     parse speclist usage;
cannam@167:     print_string (generate (check_size ()) !mode);
cannam@167:   end
cannam@167: 
cannam@167: let _ = main()