module DSP.Filter.IIR.Prony (prony) where
import Data.Array
import Matrix.Matrix
import Matrix.LU
prony :: Int
-> Int
-> Array Int Double
-> (Array Int Double, Array Int Double)
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] ]