hsignal 0.1.0.5 → 0.1.1.1
raw patch · 4 files changed
+174/−55 lines, 4 filesdep +hmatrix-gsl-statsdep −hstatisticsdep ~hmatrix
Dependencies added: hmatrix-gsl-stats
Dependencies removed: hstatistics
Dependency ranges changed: hmatrix
Files
- CHANGES +9/−0
- hsignal.cabal +5/−5
- lib/Numeric/Signal/EEG/BDF.hs +25/−12
- lib/Numeric/Signal/Multichannel.hs +135/−38
CHANGES view
@@ -12,3 +12,12 @@ 0.1.0.5: fixed stack overflow in loadBDF by making [Double] -> Vector Double step strict++0.1.1.1:+ updated to reflect change from using hstatistics to hmatrix-gsl-stats+ added detrend/filter/slice to Multichannel+ changed mapConcurrently argument order+ changed getChannels to return array+ added mapArrayConcurrently+ improved performance of BDF readData+ changed BDF to store data as Floats, not Double
hsignal.cabal view
@@ -1,5 +1,5 @@ Name: hsignal-Version: 0.1.0.5+Version: 0.1.1.1 License: GPL License-file: LICENSE Author: Vivian McPhail@@ -23,11 +23,11 @@ library Build-Depends: base >= 3 && < 5,- storable-complex, haskell98, mtl, bytestring,- hmatrix >= 0.9.3,- hstatistics >= 0.1.0.5,+ haskell98, mtl, array,- ghc-binary+ bytestring, storable-complex, ghc-binary,+ hmatrix >= 0.10.0,+ hmatrix-gsl-stats >= 0.1.1.4 Extensions: ForeignFunctionInterface
lib/Numeric/Signal/EEG/BDF.hs view
@@ -129,7 +129,7 @@ , prefilter :: ![String] , samples :: ![Int] , reserved :: ![String] - , data_ :: ![Vector Double] + , data_ :: ![Vector Float] } --deriving(Show) getString :: Int -> BSM String @@ -164,7 +164,7 @@ --reverseBits w = foldl setBit 0 . foldl ((++) . \b -> if testBit w b then [9-b] else []) [] [1..8] reverseBits w = foldr (\b r -> if not $ testBit w b then setBit r (7-b) else r) 0 [0..7] -get24Bit :: (Int -> Double) -> BSM Double +get24Bit :: (Int -> Float) -> BSM Float get24Bit f = do b1 <- getN 1 b2 <- getN 1 @@ -173,18 +173,32 @@ return $ f $ (to32 b3) `shiftL` 16 .|. (to32 b2) `shiftL` 8 .|. (to32 b1) where to32 = fromIntegral {- . reverseBits -} . BS.head -readRecord :: (Int -> Double) -> Int -> BSM (Vector Double) +readRecord :: (Int -> Float) -> Int -> BSM (Vector Float) readRecord f s = do m <- replicateM s $ get24Bit f return $! fromList m -readRecordBlock :: (Int -> Int -> Double) -> [(Int,Int)] -> BSM [Vector Double] +readRecordBlock :: (Int -> Int -> Float) -> [(Int,Int)] -> BSM [Vector Float] readRecordBlock f ss = mapM (\(i,s) -> readRecord (f i) s) ss -convert :: [Int] -> [Int] -> [Int] -> [Int] -> Int -> Int -> Double +convert :: [Int] -> [Int] -> [Int] -> [Int] -> Int -> Int -> Float convert p_min p_max d_min d_max i = \x -> (fromIntegral x) - (fromIntegral (d_min !! i))*(fromIntegral ((p_max !! i) - (p_min !! i)))/(fromIntegral ((d_max !! i) - (d_min !! i))) -readData :: (Int -> Int -> Double) -> Int -> [(Int,Int)] -> BSM [Vector Double] +readData :: (Int -> Int -> Float) -> Int -> [(Int,Int)] -> BSM [Vector Float] +readData f rs ss = do + lift $ putStrLn $ "channels: " ++ (show $ length ss) ++ ", records: " ++ (show rs) ++ ", samples per record: " ++ (show $ snd $ head ss) + d <- mapM (\x -> do + lift $ putStrLn $ "Record: " ++ show x + readRecordBlock f ss) [1..rs] + -- let v = rotate d + -- lift $ putStrLn $ "vectors: " ++ (show $ length v) + -- lift $ putStrLn $ "slices: " ++ (show $ length $ head v) + return $! map join $! rotate d + where rotate [] = [] + rotate xs@((_:[]):_) = [concat xs] + rotate ((x:xs):xss) = (x : (map head xss)) : (rotate (xs : (map tail xss))) + +{- readData _ 0 _ = error "readData, zeroth record requested" readData f 1 ss = do lift $ putStrLn $ "data record: 1" @@ -195,7 +209,7 @@ ds <- readData f n ss let result = zipWith (\x y -> join [x,y]) bs ds return $ result - +-} readBDF :: BSM (Maybe BDF) readBDF = do @@ -259,14 +273,13 @@ , data_ = data_' } -loadBDF :: FilePath -> IO (Maybe (M.Multichannel Double)) +loadBDF :: FilePath -> IO (Maybe (M.Multichannel Float)) loadBDF fn = do bs <- BS.readFile fn (bdf,bs') <- runStateT readBDF bs m <- case bdf of (Just b) -> do - let s = (num_records b) * (head $ samples b) `div` (duration b) - return $ Just $ M.fromList s 24 (data_ b) + return $ Just $ M.fromList ((head $ samples b) * (num_records b)) 24 (data_ b) _ -> do putStrLn "File not read" return Nothing @@ -275,13 +288,13 @@ return m {- -getDPConversion :: BDF -> [Word32 -> Double] +getDPConversion :: BDF -> [Word32 -> Float] getDPConversion b = let chan = channels b tup4 = zip4 (phys_min b) (phys_max b) (dig_min b) (dig_max b) in map (\(a,b,c,d) -> \x -> ((fromIntegral x) - (fromIntegral c))*(fromIntegral (b-a))/(fromIntegral (d-c))) tup4 -writeLine :: [Double] -> String +writeLine :: [Float] -> String writeLine = ((++ "\n") . concat . intersperse " " . map show) main = do
lib/Numeric/Signal/Multichannel.hs view
@@ -1,4 +1,5 @@ {-# OPTIONS_GHC -fglasgow-exts #-}+{-# OPTIONS_GHC -XUndecidableInstances #-} ----------------------------------------------------------------------------- -- | -- Module : Numeric.Signal.Multichannel@@ -18,14 +19,17 @@ module Numeric.Signal.Multichannel ( Multichannel,readMultichannel,writeMultichannel, fromList,- sampling_rate,precision,channels,+ sampling_rate,precision,channels,samples,+ detrended,filtered, getChannel,getChannels,- mapConcurrently+ mapConcurrently,+ detrend,filter,+ slice ) where ----------------------------------------------------------------------------- ---import qualified Numeric.Signal as S+import qualified Numeric.Signal as S --import Complex @@ -50,54 +54,103 @@ --import qualified Numeric.GSL.Fourier as F ---import Prelude hiding(filter)+import Prelude hiding(filter) +--import Control.Monad+++{-+-------------------------------------------------------------------++instance (Binary a, Storable a) => Binary (Vector a) where+ put v = do+ let d = dim v+ put d+ mapM_ (\i -> put $ v @> i) [0..(d-1)]+ get = do+ d <- get+ xs <- replicateM d get+ return $ fromList xs++-------------------------------------------------------------------+-}+ ----------------------------------------------------------------------------- -- | data type with multiple channels data Multichannel a = MC {- _sampling_rate :: Int -- ^ sampling rate- , _precision :: Int -- ^ bits of precision- , _channels :: Int -- ^ number of channels- , _length :: Int -- ^ length in samples+ _sampling_rate :: Int -- ^ sampling rate+ , _precision :: Int -- ^ bits of precision+ , _channels :: Int -- ^ number of channels+ , _length :: Int -- ^ length in samples+ , _detrended :: Bool -- ^ was the data detrended?+ , _filtered :: Maybe (Int,Int) -- ^ if filtered the passband , _data :: I.Array Int (Vector a) -- ^ data } ----------------------------------------------------------------------------- -instance Binary (Multichannel Double) where- put (MC s p c l d) = do- put s- put p- put c- put l- put d+class Vectors a where+ vMax :: Vector a -> a+ vMin :: Vector a -> a++instance Vectors Float where+ vMax = vectorFMax+ vMin = vectorFMin++instance Vectors Double where+ vMax = vectorMax+ vMin = vectorMin++instance (Binary a, Storable a, + Ord a, RealFrac a,+ Vectors a) => Binary (Multichannel a) where+ put (MC s p c l de f d) = do+ put s+ put p+ put c+ put l+ put de+ put f+ put $! fmap convert d+ where convert v = let (mi,ma) = (vMin v,vMax v)+ v' = mapVector (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word32))) v+ in (mi,ma,(v' :: Vector Word32)) +++ get = do s <- get p <- get c <- get l <- get- d <- get- return (MC s p c l d)+ de <- get+ f <- get+ d <- (get :: Get (I.Array Int (a,a,Vector Word32)))+ return $! (MC s p c l de f (seq d (fmap convert) d))+ where convert (mi,ma,v) = mapVector (\x -> ((fromIntegral x) :: a) / (fromIntegral (maxBound :: Word32)) * (ma - mi) + mi) v ----------------------------------------------------------------------------- -readMultichannel :: FilePath -> IO (Multichannel Double)+readMultichannel :: (Binary a, Storable a, + Ord a, RealFrac a,+ Vectors a) => FilePath -> IO (Multichannel a) readMultichannel = decodeFile -writeMultichannel :: FilePath -> Multichannel Double -> IO ()+writeMultichannel :: (Binary a, Storable a, + Ord a, RealFrac a,+ Vectors a) => FilePath -> Multichannel a -> IO () writeMultichannel = encodeFile ----------------------------------------------------------------------------- -- | create a multichannel data type-fromList :: Storable a =>- Int -- ^ sampling rate- -> Int -- ^ bits of precision- -> [Vector a] -- ^ data- -> Multichannel a -- ^ datatype+fromList :: Storable a => Int -- ^ sampling rate+ -> Int -- ^ bits of precision+ -> [Vector a] -- ^ data+ -> Multichannel a -- ^ datatype fromList s p d = let c = length d- in MC s p c (dim $ head d) (I.listArray (1,c) d)+ in MC s p c (dim $ head d) False Nothing (I.listArray (1,c) d) -- | the sampling rate sampling_rate :: Multichannel a -> Int@@ -120,23 +173,67 @@ getChannel c d = (_data d) I.! c -- | extract all channels-getChannels :: Multichannel a -> [Vector a]-getChannels d = I.elems $ _data d+getChannels :: Multichannel a -> I.Array Int (Vector a)+getChannels d = _data d +-- | was the data detrended?+detrended :: Multichannel a -> Bool+detrended = _detrended++-- | was the data filtered?+filtered :: Multichannel a -> Maybe (Int,Int)+filtered = _filtered+ ----------------------------------------------------------------------------- -- | map a function executed concurrently-mapConcurrently :: Storable b => - Multichannel a -- ^ input data- -> (Vector a -> Vector b) -- ^ the function to be mapped- -> Multichannel b -- ^ output data-mapConcurrently (MC sr p c _ d) f = unsafePerformIO $ do- results <- newMVar []- mapM_ (forkIO . applyFunction results f) $ I.assocs d- vectors <- takeMVar results- return $ MC sr p c (dim $ snd $ head vectors) (I.array (1,c) vectors)- where applyFunction results f' (i,e) = do+mapArrayConcurrently :: (a -> b) -- ^ function to map+ -> I.Array Int a -- ^ input+ -> I.Array Int b -- ^ output+mapArrayConcurrently f d = unsafePerformIO $ do+ results <- newMVar []+ mapM_ (forkIO . applyFunction results f) $ I.assocs d+ vectors <- takeMVar results+ return $ I.array (I.bounds d) vectors+ where applyFunction results f' (j,e) = do let o = f' e- modifyMVar_ results (\x -> return ((i,o):x))+ modifyMVar_ results (\x -> return ((j,o):x))++-- | 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+mapConcurrently f (MC sr p c _ de fi d) = let d' = mapArrayConcurrently f d+ in MC sr p c (dim $ d' I.! 1) de fi d'++-- | map a function+mapMC :: Storable b => (Vector a -> Vector b) -- ^ the function to be mapped + -> Multichannel a -- ^ input data+ -> Multichannel b -- ^ output data+mapMC f (MC sr p c _ de fi d) = let d' = fmap f d+ in MC sr p c (dim $ d' I.! 1) de fi d' +-----------------------------------------------------------------------------++-- | detrend the data with a specified window size+detrend :: Int -> Multichannel Double -> Multichannel Double+detrend w m = let m' = mapConcurrently (S.detrend w) m+ in m' { _detrended = True }+++-- | filter the data with the given passband+filter :: (Int,Int) -> Multichannel Double -> Multichannel Double+filter pb m = let m' = mapConcurrently (S.broadband_filter (_sampling_rate m) pb) m+ in m' { _filtered = Just pb }++-----------------------------------------------------------------------------++-- | extract a slice of the data+slice :: Int -- ^ starting sample number+ -> Int -- ^ length+ -> Multichannel Double + -> Multichannel Double+slice j w m = let m' = mapConcurrently (subVector j w) m+ in m' { _length = w }+ -----------------------------------------------------------------------------