-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Pisarenko
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains an implementation of Pisarenko Harmonic
-- Decomposition for a single real sinusoid.  For this case, eigenvalues
-- do not need to be computed.
--
-----------------------------------------------------------------------------

-- This implmentation is based off of a Matlab version by Peter
-- Kootsookos (p.kootsookos@ieee.org).

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

rss :: forall a b. (Ix a, Integral a, Num b) => Array a b -> a -> b
rss 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's method for a single sinusoid

pisarenko :: (Ix a, Integral a, Floating b) => Array a b -- ^ x
	  -> b -- ^ w

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