hsignal 0.1.1.1 → 0.1.2.1
raw patch · 7 files changed
+111/−53 lines, 7 filesdep +hstatisticsdep ~hmatrix-gsl-stats
Dependencies added: hstatistics
Dependency ranges changed: hmatrix-gsl-stats
Files
- CHANGES +9/−0
- hsignal.cabal +4/−3
- lib/Numeric/Signal.hs +14/−11
- lib/Numeric/Signal/EEG/BDF.hs +1/−1
- lib/Numeric/Signal/Internal.hs +3/−4
- lib/Numeric/Signal/Multichannel.hs +63/−21
- lib/Numeric/Signal/signal-aux.c +17/−13
CHANGES view
@@ -21,3 +21,12 @@ added mapArrayConcurrently improved performance of BDF readData changed BDF to store data as Floats, not Double++0.1.2.1:+ fixed mapConcurrently to wait for all threads+ fixed Signal.detrend window sizes + fixed Signal.detrend mathematics+ fixed broadband_filter to call the correct filter (Numeric.Signal)+ fixed signal-aux.c filter scaling error+ fixed hilbert+ changed Multichannel.fromList to createMultichannel to avoid name conflicts
hsignal.cabal view
@@ -1,5 +1,5 @@ Name: hsignal-Version: 0.1.1.1+Version: 0.1.2.1 License: GPL License-file: LICENSE Author: Vivian McPhail@@ -10,7 +10,7 @@ Description: Purely functional interface to signal processing based on hmatrix provides data types for manipulating EEG data, including reading from BDF-Category: Math+Category: Math, tested-with: GHC ==6.12.1 cabal-version: >=1.8@@ -27,7 +27,8 @@ array, bytestring, storable-complex, ghc-binary, hmatrix >= 0.10.0,- hmatrix-gsl-stats >= 0.1.1.4+ hmatrix-gsl-stats >= 0.1.1.5,+ hstatistics >= 0.2.0.5 Extensions: ForeignFunctionInterface
lib/Numeric/Signal.hs view
@@ -29,7 +29,7 @@ import Numeric.GSL.Fitting.Linear -import Complex+import Data.Complex import qualified Data.List as L @@ -53,9 +53,12 @@ -> Int -- ^ sampling rate -> Vector Double -- ^ input signal -> Vector Double -- ^ output signal-filter b a s v = let sd = (fromIntegral s) * 2.0- sc = recip $ (sd / (fromIntegral $ max (dim b) (dim a))) * sd- in scale sc $ S.filter b a v+filter b a s v = let len = dim v+ w = min s len+ start = (negate . fromList . reverse . toList . subVector 0 w) v+ finish = (negate . fromList . reverse . toList . subVector (len-w) w) v+ v' = join [start,v,finish]+ in subVector s len $ S.filter b a v' ----------------------------------------------------------------------------- @@ -100,7 +103,7 @@ -> Vector Double -- ^ input signal -> Vector Double -- ^ output signal broadband_filter s f v = let b = broadband_fir s f- in S.filter b (constant 1.0 1) v+ in filter b (scalar 1.0) s v ----------------------------------------------------------------------------- @@ -208,15 +211,15 @@ -> Vector Double -- ^ data to be detrended -> Vector Double -- ^ detrended data detrend w v = let windows = dim v `div` w- ws = takesV ((replicate windows w) ++ [dim v - (windows * w) + 1]) v+ ws = takesV ((replicate windows w) ++ [dim v - (windows * w)]) v ds = map detrend' ws windows' = (dim v - (w `div` 2)) `div` w- ws' = takesV (((w `div` 2):(replicate windows w)) ++ [dim v - (w `div` 2) - (windows' * w) + 1]) v+ ws' = takesV (((w `div` 2):(replicate windows' w)) ++ [dim v - (w `div` 2) - (windows' * w)]) v ds' = map detrend' ws'- in (join ds + join ds') / 2+ in (join ds + join ds') / 2 where detrend' x = let ln = dim x- y = linspace ln (1.0,fromIntegral ln)- (c0,c1,_,_,_,_) = linear x y- in x - (fromList [c0]) - (scale c1 x)+ t = linspace ln (1.0,fromIntegral ln)+ (c0,c1,_,_,_,_) = linear t x+ in x - (scale c1 t + scalar c0) -----------------------------------------------------------------------------
lib/Numeric/Signal/EEG/BDF.hs view
@@ -279,7 +279,7 @@ (bdf,bs') <- runStateT readBDF bs m <- case bdf of (Just b) -> do - return $ Just $ M.fromList ((head $ samples b) * (num_records b)) 24 (data_ b) + return $ Just $ M.createMultichannel (head $ samples b) 24 (data_ b) _ -> do putStrLn "File not read" return Nothing
lib/Numeric/Signal/Internal.hs view
@@ -92,13 +92,12 @@ -- | Hilbert transform with original vector as real value, transformed as imaginary hilbert :: Vector Double -> Vector (Complex Double) hilbert v = unsafePerformIO $ do- r <- createVector (dim v)- let v' = complex v+ let r = complex v -- could use (comp v) to make a complex vector in haskell rather than C- app2 signal_hilbert vec v' vec r "hilbert"+ app1 signal_hilbert vec r "hilbert" return r -foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> CInt -> PC -> IO CInt+foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> IO CInt -----------------------------------------------------------------------------
lib/Numeric/Signal/Multichannel.hs view
@@ -18,13 +18,14 @@ module Numeric.Signal.Multichannel ( Multichannel,readMultichannel,writeMultichannel,- fromList,+ createMultichannel, sampling_rate,precision,channels,samples, detrended,filtered, getChannel,getChannels, mapConcurrently, detrend,filter,- slice+ slice,+ mi_phase ) where -----------------------------------------------------------------------------@@ -34,6 +35,7 @@ --import Complex import qualified Data.Array.IArray as I+import Data.Ix import Control.Concurrent --import Control.Concurrent.MVar@@ -42,7 +44,8 @@ --import qualified Data.List as L -import Data.Packed.Vector hiding(fromList)+import Data.Packed.Vector+import Data.Packed.Matrix --import Data.Packed(Container(..)) import Data.Binary@@ -54,9 +57,14 @@ --import qualified Numeric.GSL.Fourier as F +import qualified Numeric.GSL.Histogram as H+import qualified Numeric.GSL.Histogram2D as H2++import qualified Numeric.Statistics.Information as SI+ import Prelude hiding(filter) ---import Control.Monad+import Control.Monad(replicateM) {-@@ -145,11 +153,11 @@ ----------------------------------------------------------------------------- -- | create a multichannel data type-fromList :: Storable a => Int -- ^ sampling rate- -> Int -- ^ bits of precision- -> [Vector a] -- ^ data- -> Multichannel a -- ^ datatype-fromList s p d = let c = length d+createMultichannel :: Storable a => Int -- ^ sampling rate+ -> Int -- ^ bits of precision+ -> [Vector a] -- ^ data+ -> Multichannel a -- ^ datatype+createMultichannel s p d = let c = length d in MC s p c (dim $ head d) False Nothing (I.listArray (1,c) d) -- | the sampling rate@@ -187,19 +195,30 @@ ----------------------------------------------------------------------------- -- | map a function executed concurrently-mapArrayConcurrently :: (a -> b) -- ^ function to map- -> I.Array Int a -- ^ input- -> I.Array Int b -- ^ output+mapArrayConcurrently :: Ix i => (a -> b) -- ^ function to map+ -> I.Array i a -- ^ input+ -> I.Array i b -- ^ output mapArrayConcurrently f d = unsafePerformIO $ do- results <- newMVar []- mapM_ (forkIO . applyFunction results f) $ I.assocs d- vectors <- takeMVar results+ let b = I.bounds d+ results <- replicateM (rangeSize b) newEmptyMVar+ mapM_ (forkIO . applyFunction f) $ zip results (I.assocs d)+ vectors <- mapM takeMVar results return $ I.array (I.bounds d) vectors- where applyFunction results f' (j,e) = do- let o = f' e- modifyMVar_ results (\x -> return ((j,o):x))+ where applyFunction f' (m,(j,e)) = putMVar m (j,f' e) +{- -- | map a function executed concurrently+mapListConcurrently :: (a -> b) -- ^ function to map+ -> [a] -- ^ input+ -> [b] -- ^ output+mapListConcurrently f d = unsafePerformIO $ do+ results <- replicateM (length d) newEmptyMVar+ mapM_ (forkIO . applyFunction f) zip results d+ mapM takeMVar results+ where applyFunction f' (m,e) = putMVar m (f' e)+-}++-- | map a function executed concurrently mapConcurrently :: Storable b => (Vector a -> Vector b) -- ^ the function to be mapped -> Multichannel a -- ^ input data -> Multichannel b -- ^ output data@@ -229,11 +248,34 @@ ----------------------------------------------------------------------------- -- | extract a slice of the data-slice :: Int -- ^ starting sample number+slice :: Storable a => Int -- ^ starting sample number -> Int -- ^ length- -> Multichannel Double - -> Multichannel Double+ -> Multichannel a + -> Multichannel a slice j w m = let m' = mapConcurrently (subVector j w) m in m' { _length = w }++-----------------------------------------------------------------------------++-- | calculate the mutual information of the phase between pairs of channels+-- | the lower half of the matrix displays mutual information between pairs of channels+-- | the upper half of the matrix displays mutual information between the data (row) and a sample+-- | from the same distribution as the column+mi_phase :: Multichannel Double -- ^ input data+ -> Matrix Double+mi_phase m = let d = _data m+ s = samples m+ histarray = mapArrayConcurrently (H.fromLimits 128 (-pi,pi)) d+ c = channels m+ pairs = I.array ((1,1),(c,c)) $ map (\(a,b) -> ((a,b),((a,b),d I.! a,d I.! b))) (range ((1,1),(c,c)))+ hist2array = mapArrayConcurrently (\(j,x,y) -> (j,H2.addVector (H2.emptyLimits 128 128 (-pi,pi) (-pi,pi)) x y)) pairs+ mi = mapArrayConcurrently (doMI histarray d s) hist2array+ in fromArray2D mi+ where doMI histarray d s ((x,y),h2) + | x < y = SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,d I.! y)+ | otherwise = let hypdf = H.fromHistogram (histarray I.! y)+ ys = replicate s (H.sample hypdf) + dy' = fromList ys+ in SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,dy') -----------------------------------------------------------------------------
lib/Numeric/Signal/signal-aux.c view
@@ -2,6 +2,7 @@ #include <gsl/gsl_math.h> #include <gsl/gsl_fft_real.h>+#include <gsl/gsl_fft_halfcomplex.h> #include <gsl/gsl_fft_complex.h> #include <gsl/gsl_vector.h> #include <gsl/gsl_blas.h>@@ -59,21 +60,19 @@ for (i = 0; i < vs; i++) { r[i] = 0;- for (j = 1; j < N; j++) {- if (i - j > 0) r[i] -= (l[j+1]/L)*v[i-j];+ for (j = 0; j < N; j++) {+ if (i - j > 0) r[i] -= (l[j+1])*v[i-j]; } for (j = 0; j < M; j++) {- if (i - j > 0) r[i] += (k[j+1]/K)*r[i-j];+ if (i - j > 0) r[i] += (k[j+1])*r[i-j]; } } return 0; } -int hilbert(int vs, const gsl_complex* v, int rs, gsl_complex* r)+int hilbert(int rs, gsl_complex* r) {- if (vs != rs) return 2000; // BAD_SIZE-- int s = vs;+ int s = rs; gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s); gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);@@ -81,17 +80,22 @@ // forward fourier transform gsl_fft_complex_forward ((double*)r, 1, s, wavetable, workspace); // zero negative coefficients and double positive+ int i;- for (i = 0; i < s; i++) {- if (i <= s/2) {- r[i].dat[0] *= sqrt(2.0);- r[i].dat[1] *= sqrt(2.0);+ int m = s/2;+ for (i = 1; i < s; i++) {+ if (i <= m) {+ r[i].dat[0] *= 2;+ r[i].dat[1] *= 2; }+ else if (s % 2 == 0 && i == m+1) {+ } else {- r[i].dat[0] = 0.0;- r[i].dat[1] = 0.0;+ r[i].dat[0] = 0;+ r[i].dat[1] = 0; } }+ // inverse fourier transform gsl_fft_complex_inverse ((double*)r, 1, s, wavetable, workspace);