-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.IIR.Prony
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- General case of Prony's Method where K > p+q
--
-- References: L&I, Sect 8.1; P&B, Sect 7.5; P&M, Sect 8.5.2
--
-- Notation follows L&I
--
-----------------------------------------------------------------------------

-- TODO: Handle rank deficiencies of G3 gracefully.  Can/should we
-- generate a (K/2+1) by (K/2+1) G2, and set p=q=rank(G2)?  Need SVD to
-- compute rank, though.

module DSP.Filter.IIR.Prony (prony) where

import Data.Array

import Matrix.Matrix
import Matrix.LU

{------------------------------------------------------------------------------

Case 1: K=p+q

a = array (0,p)
b = array (0,q)

g1 : q+1 by p+1
g2 : p   by p+1
g3 : p   by p

We do not define G1 and G2, but

mg2 = array ((1,1),(p,p+1)) [ ((i,j), g!(p+i+1-j)) | j <- [1..p+1], i <- [1..p] ]

prony p q g = (a,b)
    where mg3 = array ((1,1),(p,p)) [ ((i,j), g!(p+i-j)) | j <- [1..p], i <- [1..p] ]
          g1  = array (1,p) [ (i, g!(p+i)) | i <- [1..p] ]
          a'  = solve mg3 (fmap negate g1)
          a   = array (0,p) $ (0,1) : [ (i,a'!i) | i <- [1..p] ]
          b   = listArray (0,q) [ sum [ a!j * g!(i-j) | j <- [0..(min i p)] ] | i <- [0..q] ]

Test case, pg 422

g = listArray (0,6) [ 1, 18, 9, 2, 1, 2/9, 1/9 ] :: Array Int Double

------------------------------------------------------------------------------}

-- Case 2: K>p+q

-- a = array (0,p)
-- b = array (0,q)

-- g1 : q+1 by p+1
-- g2 : K-q by p+1
-- g3 : K-q by p

-- We need gi for the q<p cases because these generate zero elements in
-- G3, and this is the easiest way to take care of that.

-- mg1 = array ((1,1),(q+1,p+1)) [ ((i,j), gi (i-j)) | j <- [1..p+1], i <- [1..q+1] ]
-- mg2 = array ((1,1),(k-q,p+1)) [ ((i,j), gi (q+i-j+1)) | j <- [1..p+1], i <- [1..k-q] ]

-- | Implementation of Prony's method

prony :: Int -- ^ p
      -> Int -- ^ q
      -> Array Int Double -- ^ g[n]
      -> (Array Int Double, Array Int Double) -- ^ (b,a)

prony :: Int
-> Int -> Array Int Double -> (Array Int Double, Array Int Double)
prony Int
p Int
q Array Int Double
g = (Array Int Double
b,Array Int Double
a)
    where k :: Int
k   = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int Double
g
	  gi :: Int -> Double
gi Int
i | Int
i forall a. Ord a => a -> a -> Bool
< Int
0     = Double
0
	       | Int
i forall a. Ord a => a -> a -> Bool
> Int
k     = Double
0
	    | Bool
otherwise = Array Int Double
gforall i e. Ix i => Array i e -> i -> e
!Int
i
	  mg3 :: Array (Int, Int) Double
mg3 = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array ((Int
1,Int
1),(Int
kforall a. Num a => a -> a -> a
-Int
q,Int
p)) [ ((Int
i,Int
j), Int -> Double
gi (Int
qforall a. Num a => a -> a -> a
+Int
iforall a. Num a => a -> a -> a
-Int
j)) | Int
j <- [Int
1..Int
p], Int
i <- [Int
1..Int
kforall a. Num a => a -> a -> a
-Int
q] ]
	  g1 :: Array Int Double
g1  = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int
1,Int
kforall a. Num a => a -> a -> a
-Int
q) [ (Int
i, Int -> Double
gi (Int
qforall a. Num a => a -> a -> a
+Int
i)) | Int
i <- [Int
1..Int
kforall a. Num a => a -> a -> a
-Int
q] ]
	  a' :: Array Int Double
a'  = Array (Int, Int) Double -> Array Int Double -> Array Int Double
solve (forall i j k a.
(Ix i, Ix j, Ix k, Num a) =>
Array (i, j) a -> Array (j, k) a -> Array (i, k) a
mm_mult (forall i j a.
(Ix i, Ix j, Num a) =>
Array (i, j) a -> Array (j, i) a
m_trans Array (Int, Int) Double
mg3) Array (Int, Int) Double
mg3) (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Num a => a -> a
negate (forall i j a.
(Ix i, Ix j, Num a) =>
Array (i, j) a -> Array j a -> Array i a
mv_mult (forall i j a.
(Ix i, Ix j, Num a) =>
Array (i, j) a -> Array (j, i) a
m_trans Array (Int, Int) Double
mg3) Array Int Double
g1))
          a :: Array Int Double
a   = forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int
0,Int
p) forall a b. (a -> b) -> a -> b
$ (Int
0,Double
1) forall a. a -> [a] -> [a]
: [ (Int
i,Array Int Double
a'forall i e. Ix i => Array i e -> i -> e
!Int
i) | Int
i <- [Int
1..Int
p] ]
	  b :: Array Int Double
b   = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
q) [ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array Int Double
aforall i e. Ix i => Array i e -> i -> e
!Int
j forall a. Num a => a -> a -> a
* Int -> Double
gi (Int
iforall a. Num a => a -> a -> a
-Int
j) | Int
j <- [Int
0..(forall a. Ord a => a -> a -> a
min Int
i Int
p)] ] | Int
i <- [Int
0..Int
q] ]