-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Estimation.Spectral.MA
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains one algorithm for MA parameter estimation.  It
-- is 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.
--
-----------------------------------------------------------------------------


module DSP.Estimation.Spectral.MA (ma_durbin) where

import Data.Array
import Data.Complex

import DSP.Estimation.Spectral.AR

-- * Functions

-------------------------------------------------------------------------------
-- ma_durbin x q l
-------------------------------------------------------------------------------

-- Section 8.4 in Kay

-- | Computes an MA(q) model estimate from x using the Durbin's method
-- where l is the order of the AR process used in the algorithm

ma_durbin :: (Ix a, Integral a, RealFloat b) => Array a (Complex b)  -- ^ x
          -> a                        -- ^ q
          -> a                        -- ^ l
          -> (Array a (Complex b), b) -- ^ (a,rho)

ma_durbin :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> a -> (Array a (Complex b), b)
ma_durbin Array a (Complex b)
x a
q a
l = (Array a (Complex b)
b, b
sig2)
    where (Array a (Complex b)
b,b
_)       = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> (Array a (Complex b), b)
ar_yw Array a (Complex b)
a' a
q
 	  a' :: Array a (Complex b)
a'          = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (a
0,a
l) ((a
0,Complex b
1) forall a. a -> [a] -> [a]
: [ (a
i, Array a (Complex b)
a''forall i e. Ix i => Array i e -> i -> e
!a
i) | a
i <- [a
1..a
l] ])
          (Array a (Complex b)
a'', b
sig2) = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> (Array a (Complex b), b)
ar_yw Array a (Complex b)
x a
l