module DSP.Estimation.Frequency.Pisarenko (pisarenko) where
import DSP.Basic((^!))
import Data.Array
rss :: (Ix a, Integral a, Num b) =>
Array a b
-> a
-> b
Array a b
x a
k = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array a b
xforall i e. Ix i => Array i e -> i -> e
!(a
iforall a. Num a => a -> a -> a
+a
k) forall a. Num a => a -> a -> a
* Array a b
xforall i e. Ix i => Array i e -> i -> e
!a
i | a
i <- [a
0..(a
nforall a. Num a => a -> a -> a
-a
1forall a. Num a => a -> a -> a
-a
k)] ]
where n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a b
x) forall a. Num a => a -> a -> a
+ a
1
pisarenko :: (Ix a, Integral a, Floating b) => Array a b
-> b
pisarenko :: forall a b. (Ix a, Integral a, Floating b) => Array a b -> b
pisarenko Array a b
x = forall a. Floating a => a -> a
acos (b
alpha forall a. Fractional a => a -> a -> a
/ b
2)
where alpha :: b
alpha = (b
rss2 forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt (b
rss2forall a. Num a => a -> Int -> a
^!Int
2 forall a. Num a => a -> a -> a
+ b
8forall a. Num a => a -> a -> a
*b
rss1forall a. Num a => a -> Int -> a
^!Int
2)) forall a. Fractional a => a -> a -> a
/ (b
2forall a. Num a => a -> a -> a
*(b
rss1 forall a. Num a => a -> a -> a
+ b
eps))
rss1 :: b
rss1 = forall a b. (Ix a, Integral a, Num b) => Array a b -> a -> b
rss Array a b
x a
1
rss2 :: b
rss2 = forall a b. (Ix a, Integral a, Num b) => Array a b -> a -> b
rss Array a b
x a
2
eps :: b
eps = b
1.0e-15