diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -92,3 +92,6 @@ 0.2.2.1: fixed unwrap + +0.2.3: + added cross covariance and friends diff --git a/hsignal.cabal b/hsignal.cabal --- a/hsignal.cabal +++ b/hsignal.cabal @@ -1,5 +1,5 @@ Name: hsignal -Version: 0.2.2.1 +Version: 0.2.3 License: BSD3 License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010, 2011 diff --git a/lib/Numeric/Signal.hs b/lib/Numeric/Signal.hs --- a/lib/Numeric/Signal.hs +++ b/lib/Numeric/Signal.hs @@ -26,6 +26,9 @@ -- * Analytic Signal analytic_signal,analytic_power,analytic_phase, unwrap, + -- * Statistics + cross_covariance,cross_correlation, + auto_covariance,auto_correlation, -- * Preprocessing detrend, resize, @@ -240,6 +243,45 @@ -- | resize the vector to length n by resampling resize :: S.Filterable a => Int -> Vector a -> Vector a resize n v = S.downsample_ (dim v `div` n) v + +----------------------------------------------------------------------------- + +-- | cross covariance of two signals +-- the cross correlation is computed by dividing the result +-- by the product of the two standard deviations +cross_covariance :: S.Filterable a => + Int -- ^ maximum delay + -> Vector a -- ^ time series + -> Vector a -- ^ time series + -> (a,a,Vector a) -- ^ (sd_x,sd_y,cov_xy) +cross_covariance = S.cross_covariance_ + +-- | cross correlation of two signals +cross_correlation :: S.Filterable a => + Int -- ^ maximum delay + -> Vector a -- ^ time series + -> Vector a -- ^ time series + -> Vector a -- ^ result +cross_correlation l x y = let (sx,sy,r) = S.cross_covariance_ l x y + in mapVector (/ (sx*sy)) r + +-- | auto covariance of two signals +-- the auto correlation is computed by dividing the result +-- by the variance +auto_covariance :: S.Filterable a => + Int -- ^ maximum delay + -> Vector a -- ^ time series + -> (a,Vector a) -- ^ (var,cov_xx) +auto_covariance l v = let (sd,_,r) = cross_covariance l v v + in (sd*sd,r) + +-- | auto correlation of two signals +auto_correlation :: S.Filterable a => + Int -- ^ maximum delay + -> Vector a -- ^ time series + -> Vector a -- ^ result +auto_correlation l v = let (var,r) = auto_covariance l v + in mapVector (/ var) r ----------------------------------------------------------------------------- diff --git a/lib/Numeric/Signal/Internal.hs b/lib/Numeric/Signal/Internal.hs --- a/lib/Numeric/Signal/Internal.hs +++ b/lib/Numeric/Signal/Internal.hs @@ -20,7 +20,7 @@ Filterable(..), freqz, pwelch, - hilbert, + hilbert ) where import Data.Packed.Development(createVector,vec,app1,app2,app3,app4) @@ -79,6 +79,11 @@ polyEval_ :: Vector a -- ^ the real coefficients -> Vector (Complex Double) -- ^ the points at which to be evaluated -> Vector (Complex Double) -- ^ the values + -- | the cross covariance of two signals + cross_covariance_ :: Int -- ^ maximum delay + -> Vector a -- ^ time series + -> Vector a -- ^ time series + -> (a,a,Vector a) -- ^ (sd_x,sd_y,cross_covariance) ----------------------------------------------------------------------------- @@ -131,7 +136,8 @@ deriv_ = derivD unwrap_ = unwrapD polyEval_ = polyEval - + cross_covariance_ = crossCovarianceD + instance Filterable Float where fromDouble = single filter_ = filterF @@ -141,6 +147,7 @@ deriv_ = derivF unwrap_ = unwrapF polyEval_ c = polyEval (double c) + cross_covariance_ = crossCovarianceF ----------------------------------------------------------------------------- @@ -340,6 +347,38 @@ return r foreign import ccall "signal-aux.h unwrap_float" signal_unwrap_float :: CInt -> PF -> CInt -> PF -> IO CInt + +----------------------------------------------------------------------------- + +-- | compute the cross covariance of two signals +crossCovarianceD :: Int -> Vector Double -> Vector Double -> (Double,Double,Vector Double) +crossCovarianceD l x y = unsafePerformIO $ do + r <- createVector (2*l) + alloca $ \sx -> + alloca $ \sy -> do + app3 (signal_cross_covariance_double (fromIntegral l) sx sy) vec x vec y vec r "cross_covariance" + sx' <- peek sx + sy' <- peek sy + return (sx',sy',r) + +foreign import ccall "signal-aux.h cross_covariance_double" + signal_cross_covariance_double :: CInt -> PD -> PD -> CInt -> PD -> CInt + -> PD -> CInt -> PD -> IO CInt + +-- | compute the cross covariance of two signals +crossCovarianceF :: Int -> Vector Float -> Vector Float -> (Float,Float,Vector Float) +crossCovarianceF l x y = unsafePerformIO $ do + r <- createVector (2*l) + alloca $ \sx -> + alloca $ \sy -> do + app3 (signal_cross_covariance_float (fromIntegral l) sx sy) vec x vec y vec r "cross_covariance" + sx' <- peek sx + sy' <- peek sy + return (sx',sy',r) + +foreign import ccall "signal-aux.h cross_covariance_float" + signal_cross_covariance_float :: CInt -> PF -> PF -> CInt -> PF -> CInt + -> PF -> CInt -> PF -> IO CInt ----------------------------------------------------------------------------- diff --git a/lib/Numeric/Signal/signal-aux.c b/lib/Numeric/Signal/signal-aux.c --- a/lib/Numeric/Signal/signal-aux.c +++ b/lib/Numeric/Signal/signal-aux.c @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include @@ -349,3 +351,78 @@ return 0; } +int cross_covariance_double(int max_lag, + double* sx, + double* sy, + int xs, const double* x, + int ys, const double* y, + int rs, double* r) +{ + if (xs != ys) return 2000; // BAD_SIZE + if (rs != 2*max_lag) return 2000; // BAD_SIZE + + double mx = gsl_stats_mean(x,1,xs); + double my = gsl_stats_mean(y,1,ys); + + *sx = gsl_stats_sd(x,1,xs); + *sy = gsl_stats_sd(y,1,ys); + + int delay; + int i, j; + + double sxy; + + for (delay=-max_lag; delay < max_lag; delay++) { + sxy = 0; + for (i=0;i= xs) sxy += (x[i]-mx)*(-my); + else sxy += (x[i]-mx)*(y[j]-my); + /* or should it be: + if (j < 0 || j >= xs) continue; + else sxy += (x[i]-mx)*(y[j]-my); + */ + } + r[delay+max_lag] = sxy/xs; + } + + return 0; +} + +int cross_covariance_float(int max_lag, + float* sx, + float* sy, + int xs, const float* x, + int ys, const float* y, + int rs, double* r) +{ + if (xs != ys) return 2000; // BAD_SIZE + if (rs != 2*max_lag) return 2000; // BAD_SIZE + + float mx = gsl_stats_float_mean(x,1,xs); + float my = gsl_stats_float_mean(y,1,ys); + + *sx = gsl_stats_float_sd(x,1,xs); + *sy = gsl_stats_float_sd(y,1,ys); + + int delay; + int i, j; + + float sxy; + + for (delay=-max_lag; delay < max_lag; delay++) { + sxy = 0; + for (i=0;i= xs) sxy += (x[i]-mx)*(-my); + else sxy += (x[i]-mx)*(y[j]-my); + /* or should it be: + if (j < 0 || j >= xs) continue; + else sxy += (x[i]-mx)*(y[j]-my); + */ + } + r[delay+max_lag] = sxy/xs; + } + + return 0; +}