-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Estimation.Spectral.ARMA
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains a few algorithms for ARMA parameter estimation.
-- Algorithms are taken from Steven M. Kay, _Modern Spectral Estimation:
-- Theory and Application_, which is one of the standard texts on the
-- subject.  When possible, variable conventions are the same in the code
-- as they are found in the text.
--
-- BROKEN: DO NOT USE
--
-----------------------------------------------------------------------------

module DSP.Estimation.Spectral.ARMA (arma_mywe) where

import Data.Array
import Data.Complex

import DSP.Correlation
-- import DSP.Estimation.Spectral.MA

-- import Matrix.LU


-- * Functions

-- | THIS DOES NOT WORK

arma_mywe ::
   (RealFloat b, Integral i, Ix i) =>
   Array i (Complex b) -> i -> i -> Array i (Complex b)
arma_mywe x p q = a'
    where r = array (q-2*p+1,q+p) [ (k, rxx_u x k) | k <- [(q-2*p+1)..(q+p)] ]
 	  a' = array (1,p) [ (k, a!(p,k)) | k <- [1..p] ]
 	  a = array ((1,1),(p,p)) [ ((k,i), ak k i) | k <- [1..p], i <- [1..k] ]
 	  b = array ((1,1),(p-1,p-1)) [ ((k,i), bk k i) | k <- [1..(p-1)], i <- [1..k] ]
 	  rho = array (1,p-1) [ (k, rhok k) | k <- [1..(p-1)] ]
 	  ak 1 1             = -r!(q+1) / r!q
	  ak k i | i==k      = -(r!(q+k) + sum [ a!(k-1,l) * r!(q+k-l) | l <- [1..(k-1)] ] ) / rho!(k-1)
 		 | otherwise = a!(k-1,i) + a!(k,k) * b!(k-1,k-i)
 	  bk 1 1             = -r!(q-1) / r!q
 	  bk k i | i==k      = -(r!(q-k) + sum [ b!(k-1,l) * r!(q-k-l) | l <- [1..(k-1)] ] ) / rho!(k-1)
 		 | otherwise = b!(k-1,i) + b!(k,k) * a!(k-1,k-i)
 	  rhok 1 = (1 - a!(1,1) * b!(1,1)) * r!q
 	  rhok k = (1 - a!(k,k) * b!(k,k)) * rho!(k-1)