-----------------------------------------------------------------------------
-- |
-- 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 :: forall b i.
(RealFloat b, Integral i, Ix i) =>
Array i (Complex b) -> i -> i -> Array i (Complex b)
arma_mywe Array i (Complex b)
x i
p i
q = Array i (Complex b)
a'
    where r :: Array i (Complex b)
r = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (i
qforall a. Num a => a -> a -> a
-i
2forall a. Num a => a -> a -> a
*i
pforall a. Num a => a -> a -> a
+i
1,i
qforall a. Num a => a -> a -> a
+i
p) [ (i
k, forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
rxx_u Array i (Complex b)
x i
k) | i
k <- [(i
qforall a. Num a => a -> a -> a
-i
2forall a. Num a => a -> a -> a
*i
pforall a. Num a => a -> a -> a
+i
1)..(i
qforall a. Num a => a -> a -> a
+i
p)] ]
 	  a' :: Array i (Complex b)
a' = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (i
1,i
p) [ (i
k, Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
p,i
k)) | i
k <- [i
1..i
p] ]
 	  a :: Array (i, i) (Complex b)
a = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array ((i
1,i
1),(i
p,i
p)) [ ((i
k,i
i), i -> i -> Complex b
ak i
k i
i) | i
k <- [i
1..i
p], i
i <- [i
1..i
k] ]
 	  b :: Array (i, i) (Complex b)
b = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array ((i
1,i
1),(i
pforall a. Num a => a -> a -> a
-i
1,i
pforall a. Num a => a -> a -> a
-i
1)) [ ((i
k,i
i), i -> i -> Complex b
bk i
k i
i) | i
k <- [i
1..(i
pforall a. Num a => a -> a -> a
-i
1)], i
i <- [i
1..i
k] ]
 	  rho :: Array i (Complex b)
rho = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (i
1,i
pforall a. Num a => a -> a -> a
-i
1) [ (i
k, i -> Complex b
rhok i
k) | i
k <- [i
1..(i
pforall a. Num a => a -> a -> a
-i
1)] ]
 	  ak :: i -> i -> Complex b
ak i
1 i
1             = -Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!(i
qforall a. Num a => a -> a -> a
+i
1) forall a. Fractional a => a -> a -> a
/ Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!i
q
	  ak i
k i
i | i
iforall a. Eq a => a -> a -> Bool
==i
k      = -(Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!(i
qforall a. Num a => a -> a -> a
+i
k) forall a. Num a => a -> a -> a
+ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1,i
l) forall a. Num a => a -> a -> a
* Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!(i
qforall a. Num a => a -> a -> a
+i
kforall a. Num a => a -> a -> a
-i
l) | i
l <- [i
1..(i
kforall a. Num a => a -> a -> a
-i
1)] ] ) forall a. Fractional a => a -> a -> a
/ Array i (Complex b)
rhoforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1)
 		 | Bool
otherwise = Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1,i
i) forall a. Num a => a -> a -> a
+ Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
k,i
k) forall a. Num a => a -> a -> a
* Array (i, i) (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1,i
kforall a. Num a => a -> a -> a
-i
i)
 	  bk :: i -> i -> Complex b
bk i
1 i
1             = -Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!(i
qforall a. Num a => a -> a -> a
-i
1) forall a. Fractional a => a -> a -> a
/ Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!i
q
 	  bk i
k i
i | i
iforall a. Eq a => a -> a -> Bool
==i
k      = -(Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!(i
qforall a. Num a => a -> a -> a
-i
k) forall a. Num a => a -> a -> a
+ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array (i, i) (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1,i
l) forall a. Num a => a -> a -> a
* Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!(i
qforall a. Num a => a -> a -> a
-i
kforall a. Num a => a -> a -> a
-i
l) | i
l <- [i
1..(i
kforall a. Num a => a -> a -> a
-i
1)] ] ) forall a. Fractional a => a -> a -> a
/ Array i (Complex b)
rhoforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1)
 		 | Bool
otherwise = Array (i, i) (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1,i
i) forall a. Num a => a -> a -> a
+ Array (i, i) (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!(i
k,i
k) forall a. Num a => a -> a -> a
* Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1,i
kforall a. Num a => a -> a -> a
-i
i)
 	  rhok :: i -> Complex b
rhok i
1 = (Complex b
1 forall a. Num a => a -> a -> a
- Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
1,i
1) forall a. Num a => a -> a -> a
* Array (i, i) (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!(i
1,i
1)) forall a. Num a => a -> a -> a
* Array i (Complex b)
rforall i e. Ix i => Array i e -> i -> e
!i
q
 	  rhok i
k = (Complex b
1 forall a. Num a => a -> a -> a
- Array (i, i) (Complex b)
aforall i e. Ix i => Array i e -> i -> e
!(i
k,i
k) forall a. Num a => a -> a -> a
* Array (i, i) (Complex b)
bforall i e. Ix i => Array i e -> i -> e
!(i
k,i
k)) forall a. Num a => a -> a -> a
* Array i (Complex b)
rhoforall i e. Ix i => Array i e -> i -> e
!(i
kforall a. Num a => a -> a -> a
-i
1)