-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Covariance
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- This module contains routines to perform cross- and auto-covariance
-- These formulas can be found in most DSP textbooks.
--
-- In the following routines, x and y are assumed to be of the same
-- length.
--
-----------------------------------------------------------------------------


-- TODO: fix these routines to handle the case were x and y are different
-- lengths.

-- TODO: Cxx(X) = Var(X), but I'm not sure how the lag works into that

module DSP.Covariance (cxy, cxy_b, cxy_u, cxx, cxx_b, cxx_u) where

import Data.Array
import Data.Complex

import DSP.Correlation
import Numeric.Statistics.Moment

-- | raw cross-covariance
--
-- We define covariance in terms of correlation.
--
-- Cxy(X,Y) = E[(X - E[X])(Y - E[Y])]
--          = E[XY] - E[X]E[Y]
--          = Rxy(X,Y) - E[X]E[Y]

-- cxy x y k | k >= 0 = sum [ (x!(i+k) - xm) * ((conjugate (y!i)) - ym) | i <- [0..(n-1-k)] ]

cxy :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x
                                 -> Array a (Complex b) -- ^ y
                                 -> a                   -- ^ k
                                 -> Complex b           -- ^ C_xy[k]

cxy :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
cxy Array a (Complex b)
x Array a (Complex b)
y a
k =
   if a
k forall a. Ord a => a -> a -> Bool
>= a
0
     then forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
rxy Array a (Complex b)
x Array a (Complex b)
y a
k forall a. Num a => a -> a -> a
- Complex b
xm forall a. Num a => a -> a -> a
* Complex b
ym
     else forall a. Num a => Complex a -> Complex a
conjugate (forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
cxy Array a (Complex b)
y Array a (Complex b)
x (-a
k))
    where xm :: Complex b
xm = forall a. Fractional a => [a] -> a
mean (forall i e. Array i e -> [e]
elems Array a (Complex b)
x)
          ym :: Complex b
ym = forall a. Fractional a => [a] -> a
mean (forall a b. (a -> b) -> [a] -> [b]
map forall a. Num a => Complex a -> Complex a
conjugate (forall i e. Array i e -> [e]
elems Array a (Complex b)
y))

-- | raw auto-covariance
--
-- Cxx(X,X) = E[(X - E[X])(X - E[X])]
--          = E[XX] - E[X]E[X]
--          = Rxy(X,X) - E[X]^2

-- We define this explicitly to prevent the mean from being calculated
-- twice.

-- cxx x k | k >= 0 = sum [ (x!(i+k) - xm) * (conjugate (x!i - xm)) | i <- [0..(n-1-k)] ]

cxx :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x
                                 -> a                   -- ^ k
                                 -> Complex b           -- ^ C_xx[k]

cxx :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
cxx Array a (Complex b)
x a
k =
   if a
k forall a. Ord a => a -> a -> Bool
>= a
0
     then forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
rxx Array a (Complex b)
x a
k forall a. Num a => a -> a -> a
- forall a. Fractional a => [a] -> a
mean (forall i e. Array i e -> [e]
elems Array a (Complex b)
x) forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2::Int)
     else forall a. Num a => Complex a -> Complex a
conjugate (forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
cxx Array a (Complex b)
x (-a
k))

-- Define the biased and unbiased versions in terms of the above.

-- | biased cross-covariance

cxy_b :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x
                                 -> Array a (Complex b) -- ^ y
                                 -> a                   -- ^ k
                                 -> Complex b           -- ^ C_xy[k] \/ N

cxy_b :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
cxy_b Array a (Complex b)
x Array a (Complex b)
y a
k = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
cxy Array a (Complex b)
x Array a (Complex b)
y a
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n
    where n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
x) forall a. Num a => a -> a -> a
+ a
1

-- | unbiased cross-covariance

cxy_u :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x
                                 -> Array a (Complex b) -- ^ y
                                 -> a                   -- ^ k
                                 -> Complex b           -- ^ C_xy[k] \/ (N-k)

cxy_u :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
cxy_u Array a (Complex b)
x Array a (Complex b)
y a
k = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> Array a (Complex b) -> a -> Complex b
cxy Array a (Complex b)
x Array a (Complex b)
y a
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
n forall a. Num a => a -> a -> a
- forall a. Num a => a -> a
abs a
k)
    where n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
x) forall a. Num a => a -> a -> a
+ a
1

-- | biased auto-covariance

cxx_b :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x
                                 -> a                   -- ^ k
                                 -> Complex b           -- ^ C_xx[k] \/ N

cxx_b :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
cxx_b Array a (Complex b)
x a
k = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
cxx Array a (Complex b)
x a
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n
    where n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
x) forall a. Num a => a -> a -> a
+ a
1

-- | unbiased auto-covariance

cxx_u :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x
                                 -> a                   -- ^ k
                                 -> Complex b           -- ^ C_xx[k] \/ (N-k)

cxx_u :: forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
cxx_u Array a (Complex b)
x a
k = forall a b.
(Ix a, Integral a, RealFloat b) =>
Array a (Complex b) -> a -> Complex b
cxx Array a (Complex b)
x a
k forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
n forall a. Num a => a -> a -> a
- forall a. Num a => a -> a
abs a
k)
    where n :: a
n = forall a b. (a, b) -> b
snd (forall i e. Array i e -> (i, i)
bounds Array a (Complex b)
x) forall a. Num a => a -> a -> a
+ a
1