hsignal 0.2.2.1 → 0.2.3
raw patch · 5 files changed
+164/−3 lines, 5 files
Files
- CHANGES +3/−0
- hsignal.cabal +1/−1
- lib/Numeric/Signal.hs +42/−0
- lib/Numeric/Signal/Internal.hs +41/−2
- lib/Numeric/Signal/signal-aux.c +77/−0
CHANGES view
@@ -92,3 +92,6 @@ 0.2.2.1: fixed unwrap++0.2.3:+ added cross covariance and friends
hsignal.cabal view
@@ -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
lib/Numeric/Signal.hs view
@@ -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 -----------------------------------------------------------------------------
lib/Numeric/Signal/Internal.hs view
@@ -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 -----------------------------------------------------------------------------
lib/Numeric/Signal/signal-aux.c view
@@ -7,6 +7,8 @@ #include <gsl/gsl_vector.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_poly.h>+#include <gsl/gsl_statistics_double.h>+#include <gsl/gsl_statistics_float.h> #include <stdio.h> @@ -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;i++) {+ j = i+delay;+ if (j < 0 || j >= 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;i++) {+ j = i+delay;+ if (j < 0 || j >= 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;+}