Levinson-Durvin method (estimation of AR coefficients)

levinson.ml

(** [levinson r] computes AR coefficients by using Levinson-Durvin method.
 @param r values of autocorrelation function: [[r(0); r(1); r(2); ...; r(M)]]
 @return [([ar(1); ar(2); ...; ar(M)], sigma2)] where [ar(i)] is the [i]-th
         coefficient of AR([M]) and [sigma2] is variance of errors.
 *)
let levinson r =
  let rec loop r ar sigma2 = function
    | [] -> (List.rev ar, sigma2)
    | rm :: rest ->
       let rev_ar = List.rev ar in
       let s = List.fold_left2 (fun acc x y -> acc +. x *. y) 0. rev_ar r in
       let arm = (rm -. s) /. sigma2 in
       let ar_rest = List.map2 (fun x y -> x -. arm *. y) ar rev_ar in
       let sigma2 = sigma2 *. (1. -. arm *. arm) in
       loop (rm :: r) (arm :: ar_rest) sigma2 rest
  in
  match r with
  | [] -> invalid_arg "levinson: length r < 0"
  | r0 :: r -> loop [] [] r0 r