-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Estimation.Frequency.WLP
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains a few algorithms for weighted linear predictors
-- for estimating the frequency of a complex sinusoid in noise.
--
-----------------------------------------------------------------------------

module DSP.Estimation.Frequency.WLP (wlp, lrp, kay, lw, ckq,) where

import DSP.Basic((^!))
import Data.Array
import Data.Complex

-- | The weighted linear predictor form of the frequency estimator

wlp :: (Ix a, Integral a, RealFloat b) => Array a b -- ^ window
    -> Array a (Complex b) -- ^ z
    -> b -- ^ w

wlp :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b -> Array a (Complex b) -> b
wlp Array a b
w Array a (Complex b)
z = forall a. RealFloat a => Complex a -> a
phase (forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ (Array a b
wforall i e. Ix i => Array i e -> i -> e
!a
t forall a. a -> a -> Complex a
:+ b
0) forall a. Num a => a -> a -> a
* Array a (Complex b)
zforall i e. Ix i => Array i e -> i -> e
!a
t forall a. Num a => a -> a -> a
* forall a. Num a => Complex a -> Complex a
conjugate (Array a (Complex b)
zforall i e. Ix i => Array i e -> i -> e
!(a
tforall a. Num a => a -> a -> a
-a
1)) | a
t <- [a
1..(a
nforall a. Num a => a -> a -> a
-a
1)] ])
    where n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
z) forall a. Num a => a -> a -> a
+ a
1

-- | WLP using Lank, Reed, and Pollon's window

lrp :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z
    -> b -- ^ w

lrp :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> b
lrp = forall a b.
(Integral a, Ix a, RealFloat b) =>
(b -> a -> b -> b) -> Array a (Complex b) -> b
processArray (\b
n a
_ b
_ -> forall a. Fractional a => a -> a
recip (b
nforall a. Num a => a -> a -> a
-b
1))

-- | WLP using kay's window

kay :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z
    -> b -- ^ w

kay :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> b
kay = forall a b.
(Integral a, Ix a, RealFloat b) =>
(b -> a -> b -> b) -> Array a (Complex b) -> b
processArray (\b
n a
_ b
t -> forall b. Fractional b => b -> b -> b
kayWin b
n b
t)

kayWin :: Fractional b => b -> b -> b
kayWin :: forall b. Fractional b => b -> b -> b
kayWin b
n b
t = b
6forall a. Num a => a -> a -> a
*b
tforall a. Num a => a -> a -> a
*(b
nforall a. Num a => a -> a -> a
-b
t) forall b. Fractional b => b -> b -> b
/ (b
nforall a. Num a => a -> a -> a
*(b
nforall a. Num a => a -> Int -> a
^!Int
2forall a. Num a => a -> a -> a
-b
1))

-- | WLP using Lovell and Williamson's window

lw :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z
    -> b -- ^ w

lw :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> b
lw Array a (Complex b)
z =
   forall a b.
(Integral a, Ix a, RealFloat b) =>
(b -> a -> b -> b) -> Array a (Complex b) -> b
processArray
      (\ b
n a
ti b
t -> forall b. Fractional b => b -> b -> b
kayWin b
n b
t forall b. Fractional b => b -> b -> b
/
          (forall a. RealFloat a => Complex a -> a
magnitude (Array a (Complex b)
zforall i e. Ix i => Array i e -> i -> e
!a
ti) forall a. Num a => a -> a -> a
* forall a. RealFloat a => Complex a -> a
magnitude (Array a (Complex b)
zforall i e. Ix i => Array i e -> i -> e
!(a
tiforall a. Num a => a -> a -> a
-a
1)))) Array a (Complex b)
z

-- | WLP using Clarkson, Kootsookos, and Quinn's window

ckq :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z
    -> b -- ^ rho
    -> b -- ^ sigma
    -> b -- ^ w

ckq :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> b -> b -> b
ckq Array a (Complex b)
z b
rho b
sig =
   forall a b.
(Integral a, Ix a, RealFloat b) =>
(b -> a -> b -> b) -> Array a (Complex b) -> b
processArray
      (\b
n ->
         let num :: b -> b
num b
t = forall a. Floating a => a -> a
sinh (b
n forall a. Num a => a -> a -> a
* b
th) forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sinh (b
t forall a. Num a => a -> a -> a
* b
th) forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
sinh ((b
nforall a. Num a => a -> a -> a
-b
t) forall a. Num a => a -> a -> a
* b
th)
             den :: b
den = (b
nforall a. Num a => a -> a -> a
-b
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh (b
n forall a. Num a => a -> a -> a
* b
th) forall a. Num a => a -> a -> a
-
                      b
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh (b
0.5 forall a. Num a => a -> a -> a
* b
n forall a. Num a => a -> a -> a
* b
th)
                        forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh (b
0.5 forall a. Num a => a -> a -> a
* (b
nforall a. Num a => a -> a -> a
-b
1) forall a. Num a => a -> a -> a
* b
th) forall b. Fractional b => b -> b -> b
/ forall a. Floating a => a -> a
sinh (b
0.5 forall a. Num a => a -> a -> a
* b
th)
             sigRho2 :: b
sigRho2 = (b
sig forall b. Fractional b => b -> b -> b
/ b
rho) forall a. Num a => a -> Int -> a
^! Int
2
             th :: b
th = forall a. Floating a => a -> a
log (b
1 forall a. Num a => a -> a -> a
+ b
sigforall a. Num a => a -> Int -> a
^!Int
2 forall b. Fractional b => b -> b -> b
/ b
rhoforall a. Num a => a -> Int -> a
^!Int
2 forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt ((b
sigRho2forall a. Num a => a -> a -> a
+b
1) forall a. Num a => a -> a -> a
* b
sigRho2))
         in  \a
_ b
t -> b -> b
num b
t forall b. Fractional b => b -> b -> b
/ b
den) Array a (Complex b)
z

processArray ::
   (Integral a, Ix a, RealFloat b) =>
   (b -> a -> b -> b) -> Array a (Complex b) -> b
processArray :: forall a b.
(Integral a, Ix a, RealFloat b) =>
(b -> a -> b -> b) -> Array a (Complex b) -> b
processArray b -> a -> b -> b
f Array a (Complex b)
z =
  let n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
z) forall a. Num a => a -> a -> a
+ a
1
      g :: a -> b -> b
g = b -> a -> b -> b
f (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n)
      bnd :: (a, a)
bnd = (a
1,a
nforall a. Num a => a -> a -> a
-a
1)
  in  forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a b -> Array a (Complex b) -> b
wlp (forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (a, a)
bnd
              (forall a b. (a -> b) -> [a] -> [b]
map (\a
t -> a -> b -> b
g a
t (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
t)) (forall a. Ix a => (a, a) -> [a]
range (a, a)
bnd))) Array a (Complex b)
z