-----------------------------------------------------------------------------
-- |
-- Module      :  Matrix.Cholesky
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains an implementation of Levinson-Durbin recursion.
--
-----------------------------------------------------------------------------

module Matrix.Levinson (levinson) where

import Data.Array
import Data.Complex

-- * Functions

-- Section 6.3.3 in Kay, formulas 6.46--6.48

-- TODO: rho is typing as complex, but it is real
-- TODO: add stepdown function
-- TODO: some applications may want all model estimations from [1..p]

-- | levinson takes an array, r, of autocorrelation values, and a
-- model order, p, and returns an array, a, of the model estimate and
-- rho, the noise power.

levinson :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ r
	                                    -> a                   -- ^ p
					    -> (Array a (Complex b),b) -- ^ (a,rho)

levinson r p = (array (1,p) [ (k, a!(p,k)) | k <- [1..p] ], realPart (rho!p))
    where a   = array ((1,1),(p,p)) [ ((k,i), ak k i) | k <- [1..p], i <- [1..k] ]
	  rho = array (1,p) [ (k, rhok k) | k <- [1..p] ]
	  ak 1 1             = -r!1 / r!0
	  ak k i | k==i      = -(r!k + sum [ a!(k-1,l) * r!(k-l) | l <- [1..(k-1)] ]) / rho!(k-1)
		 | otherwise = a!(k-1,i) + a!(k,k) * (conjugate (a!(k-1,k-i)))
	  rhok 1 = (1 - (abs (a!(1,1)))^(2::Int)) * r!0
	  rhok k = (1 - (abs (a!(k,k)))^(2::Int)) * rho!(k-1)

-- r = array (0,2) [ (0, (2.0 :+ 0.0)), (1, ((-1.0) :+ 1.0)), (2, (0.0 :+ 0.0)) ]
-- a = fst (levinson r 2)

-- verify = a == array (1,2) [(1,1.0 :+ (-1.0)),(2,0.0 :+ (-1.0))]