diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -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 diff --git a/hsignal.cabal b/hsignal.cabal --- a/hsignal.cabal +++ b/hsignal.cabal @@ -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 diff --git a/lib/Numeric/Signal/EEG/BDF.hs b/lib/Numeric/Signal/EEG/BDF.hs --- a/lib/Numeric/Signal/EEG/BDF.hs +++ b/lib/Numeric/Signal/EEG/BDF.hs @@ -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 diff --git a/lib/Numeric/Signal/Multichannel.hs b/lib/Numeric/Signal/Multichannel.hs --- a/lib/Numeric/Signal/Multichannel.hs +++ b/lib/Numeric/Signal/Multichannel.hs @@ -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 } + -----------------------------------------------------------------------------