dsp-0.1: DSP/Correlation.hs

-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Correlation
-- 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-correlation.
-- These formulas can be found in most DSP textbooks.
-- 
-- In the following routines, x and y are assumed to be of the same
-- length.
--
-----------------------------------------------------------------------------

module DSP.Correlation (rxy, rxy_b, rxy_u, rxx, rxx_b, rxx_u) where

import Data.Array
import Data.Complex

-- * Functions

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

-- | raw cross-correllation

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

rxy x y k | k >= 0 = sum [ x!(i+k) * (conjugate (y!i)) | i <- [0..(n-1-k)] ]
          | k < 0  = conjugate (rxy y x (-k))
    where n = snd (bounds x) + 1

-- | biased cross-correllation

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

rxy_b x y k = (rxy x y k) / (fromIntegral n)
    where n = snd (bounds x) + 1

-- | unbiased cross-correllation

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

rxy_u x y k = (rxy x y k) / (fromIntegral (n-(abs k)))
    where n = snd (bounds x) + 1

-- autocorrellation

-- We define autocorrelation in terms of the cross correlation routines.

-- | raw auto-correllation

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

rxx   x k = rxy   x x k

-- | biased auto-correllation

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

rxx_b x k = rxy_b x x k

-- | unbiased auto-correllation

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

rxx_u x k = rxy_u x x k

----------------------------------------------------------------------------
-- test routines
----------------------------------------------------------------------------

x = array (0,4) [ (0, 1 :+ 0), 
		  (1, 0 :+ 1), 
		  (2, (-1) :+  0), 
		  (3, 0 :+ (-1)), 
		  (4, 1 :+ 0) ]

y = array (0,4) [ (0, 1 :+ 0), 
		  (1, (-1) :+ 0), 
		  (2, 1 :+ 0), 
		  (3, (-1) :+ 0), 
		  (4, 1 :+ 0) ]

r = map (rxy_b x y) [ 0, 1, 2 ]

verify = r == [ (0.2 :+ 0.0), (0.0 :+ 0.0), (0.0 :+ 0.2) ]