hsignal 0.1.0.3 → 0.1.0.4
raw patch · 7 files changed
+377/−10 lines, 7 filesdep +bytestringdep +mtldep ~hstatistics
Dependencies added: bytestring, mtl
Dependency ranges changed: hstatistics
Files
- CHANGES +3/−0
- hsignal.cabal +7/−4
- lib/Numeric/Signal.hs +2/−1
- lib/Numeric/Signal/EEG.hs +21/−0
- lib/Numeric/Signal/EEG/BDF.hs +297/−0
- lib/Numeric/Signal/Internal.hs +5/−2
- lib/Numeric/Signal/Multichannel.hs +42/−3
CHANGES view
@@ -6,3 +6,6 @@ 0.1.0.3: added Numeric.Signal.Multichannel++0.1.0.4:+ added Numeric.Signal.EEG
hsignal.cabal view
@@ -1,5 +1,5 @@ Name: hsignal-Version: 0.1.0.3+Version: 0.1.0.4 License: GPL License-file: LICENSE Author: Alexander Vivian Hugh McPhail@@ -8,6 +8,8 @@ Homepage: http://code.haskell.org/hsignal Synopsis: Signal processing Description: Purely functional interface to signal processing based on hmatrix+ provides data types for manipulating EEG data,+ including reading from BDF Category: Math tested-with: GHC ==6.12.1 @@ -21,9 +23,9 @@ library Build-Depends: base >= 3 && < 5,- storable-complex, haskell98,+ storable-complex, haskell98, mtl, bytestring, hmatrix >= 0.9.3,- hstatistics,+ hstatistics >= 0.1.0.5, array Extensions: ForeignFunctionInterface@@ -31,7 +33,9 @@ hs-source-dirs: lib Exposed-modules: Numeric.Signal Numeric.Signal.Multichannel+ Numeric.Signal.EEG other-modules: Numeric.Signal.Internal+ Numeric.Signal.EEG.BDF C-sources: lib/Numeric/Signal/signal-aux.c ghc-prof-options: -auto@@ -39,7 +43,6 @@ ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans -fno-warn-unused-binds- -threaded if os(OSX) extra-lib-dirs: /opt/local/lib/
lib/Numeric/Signal.hs view
@@ -39,6 +39,7 @@ import Numeric.GSL.Vector import Numeric.LinearAlgebra.Algorithms+--import Numeric.LinearAlgebra.Interface import qualified Numeric.GSL.Fourier as F @@ -216,6 +217,6 @@ where detrend' x = let ln = dim x y = linspace ln (1.0,fromIntegral ln) (c0,c1,_,_,_,_) = linear x y- in x - c0 - c1 * x+ in x - (fromList [c0]) - (scale c1 x) -----------------------------------------------------------------------------
+ lib/Numeric/Signal/EEG.hs view
@@ -0,0 +1,21 @@+{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Signal.EEG +-- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com +-- Stability : provisional +-- Portability : +-- +-- EEG functions +-- +----------------------------------------------------------------------------- + +module Numeric.Signal.EEG ( + loadBDF + ) where + + +import Numeric.Signal.EEG.BDF
+ lib/Numeric/Signal/EEG/BDF.hs view
@@ -0,0 +1,297 @@+{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Signal.EEG.BDF +-- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com +-- Stability : provisional +-- Portability : +-- +-- EEG BDF data file functions +-- +----------------------------------------------------------------------------- + +module Numeric.Signal.EEG.BDF ( + loadBDF + ) where + +import qualified Data.ByteString as BS + +--import Data.Char +--import Data.List +import Data.Word +import Data.Bits +--import Data.Array.Storable +import Data.Packed.Vector + +import qualified Numeric.Signal.Multichannel as M + +import Control.Monad hiding(join) +import Control.Monad.State hiding(join) + +--import System.IO + +--import Debug.Trace + +{- BDF Header format: http://www.biosemi.com/faq/file_format.htm +Length in bytes BDF Header: EDF Header: Description +8 bytes Byte 1: "255" (non ascii) Byte 1: "0" (ASCII) Identification code +Bytes 2-8 : "BIOSEMI" (ASCII) Bytes 2-8 : " "(ASCII) +80 bytes +User text input (ASCII) + Local subject identification +80 bytes +User text input (ASCII) + Local recording identification +8 bytes +dd.mm.yy (ASCII) + Startdate of recording +8 bytes +hh.mm.ss (ASCII) + Starttime of recording +8 bytes +(ASCII) + Number of bytes in header record +44 bytes +"24BIT" (ASCII) + +"BIOSEMI" (ASCII) + Version of data format. +8 bytes +(ASCII) + Number of data records "-1" if unknown +8 bytes +e.g.: "1" (ASCII) + Duration of a data record, in seconds +4 bytes +e.g.: "257" or "128" (ASCII) + Number of channels (N) in data record +N x 16 bytes +e.g.: "Fp1", "Fpz", "Fp2", etc (ASCII) + Labels of the channels +N x 80 bytes +e.g.: "active electrode", "respiration belt" (ASCII) + Transducer type +N x 8 bytes +e.g.: "uV", "Ohm" (ASCII) + Physical dimension of channels +N x 8 bytes e.g.: "-262144" (ASCII) e.g.: "-32768" (ASCII) Physical minimum in units of physical dimension +N x 8 bytes e.g.: "262143" (ASCII) e.g.: "32767" (ASCII) Physical maximum in units of physical dimension +N x 8 bytes e.g.: "-8388608" (ASCII) e.g.: "-32768" (ASCII) Digital minimum +N x 8 bytes e.g.: "8388607" (ASCII) e.g.: "32767" (ASCII) Digital maximum +N x 80 bytes e.g.: "HP:DC; LP:410" e.g.: "HP:0,16; LP:500" Prefiltering +N x 8 bytes +For example: "2048" (ASCII) + +Number of samples in each data record +(Sample-rate if Duration of data record = "1") +N x 32 bytes +(ASCII) + Reserved +-} + +type BSM a = StateT BS.ByteString IO a + +getN :: Int -> BSM BS.ByteString +getN n = do + bs <- get + let (v,bs') = BS.splitAt n bs + put bs' + return v + +data Date = Date { day :: Int, month :: Int, year :: Int } + deriving(Eq,Ord,Show) + +data Time = Time { hour :: Int, minute :: Int, second :: Int } + deriving(Eq,Ord,Show) + +data BDF = BDF { + id_ :: Word8 + , type_ :: String + , subject :: String + , recording :: String + , date :: Date + , time :: Time + , head_bytes :: Int + , data_version :: String + , num_records :: Int + , duration :: Int + , channels :: Int + , chan_labels :: [String] + , tran_type :: [String] + , dimensions :: [String] + , phys_min :: [Int] + , phys_max :: [Int] + , dig_min :: [Int] + , dig_max :: [Int] + , prefilter :: [String] + , samples :: [Int] + , reserved :: [String] + , data_ :: [Vector Double] + } --deriving(Show) + +getString :: Int -> BSM String +getString n = do + bs <- getN n + return $ (reverse . dropWhile (== ' ') . tail . reverse . tail . show) bs + +getInt :: Int -> BSM Int +getInt n = do + s <- getString n + return $ read s + +getDate = do + s <- getString 8 + return $ strToDate' s + +strToDate' :: String -> Date +strToDate' (d1:d2:'.':m1:m2:'.':y1:y2:[]) = Date (read [d1,d2]) (read [m1,m2]) (read [y1,y2]) +strToDate' _ = error "strToDate" + +getTime = do + s <- getString 8 + return $ strToTime' s + +strToTime' :: String -> Time +strToTime' (h1:h2:'.':m1:m2:'.':s1:s2:[]) = Time (read [h1,h2]) (read [m1,m2]) (read [s1,s2]) +strToTime' _ = error "strToTime" + +reverseBits :: Word8 -> Word8 +--reverseBits w = foldl setBit 0 $ snd . unzip . filter ((== False) . fst) $ zip (map (testBit w) [1..8]) $ reverse [1..8] +--reverseBits = foldl setBit 0 . snd . unzip . filter ((False ==) . fst) . ($ reverse [0..7]) . zip . flip map [0..7] . testBit +--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 :: BSM Word32 +get24Bit f = do + b1 <- getN 1 + b2 <- getN 1 + b3 <- getN 1 +-- return $ (fromIntegral $ BS.head b1) `shiftL` 16 .|. (fromIntegral $ BS.head b2) `shiftL` 8 .|. (fromIntegral $ BS.head b3) + return $ f $ (to32 b3) `shiftL` 16 .|. (to32 b2) `shiftL` 8 .|. (to32 b1) + where to32 = fromIntegral {- . reverseBits -} . BS.head + +--readRecord :: Int -> BSM [Word32] +readRecord f s = do + m <- replicateM s $ get24Bit f + return $ fromList m +--readRecordBlock :: [(Int,Int)] -> BSM [[Word32]] +readRecordBlock f ss = mapM (\(i,s) -> readRecord (f i) s) ss + +convert :: [Int] -> [Int] -> [Int] -> [Int] -> Int -> Int -> Double +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] -> BSM [[Word32]] +readData _ 0 _ = error "readData, zeroth record requested" +readData f 1 ss = readRecordBlock f ss +readData f (n+1) ss = do + bs <- readRecordBlock f ss + ds <- readData f n ss + return $ zipWith (\x y -> join [x,y]) bs ds + +readBDF :: BSM (Maybe BDF) +readBDF = do + id_' <- getN 1 + type_' <- getString 7 + if (not $ BS.head id_' == 255 && type_' == "BIOSEMI") + then do + lift $ putStrLn "Error: File not BDF Format" + return Nothing + else do + subject' <- getString 80 + recording' <- getString 80 + date' <- getDate + time' <- getTime + head_bytes' <- getInt 8 + data_version' <- getString 44 + num_records' <- getInt 8 + if (num_records' == -1) + then do + lift $ putStrLn "This file is probably a valid BDF file..." + lift $ putStrLn " but this program cannot grok an unspecified" + lift $ putStrLn " number of records" + lift $ putStrLn "So complain to the software author" + return Nothing + else do + duration' <- getInt 8 + channels' <- getInt 4 + chan_labels' <- replicateM channels' $ getString 16 + tran_type' <- replicateM channels' $ getString 80 + dimensions' <- replicateM channels' $ getString 8 + phys_min' <- replicateM channels' $ getInt 8 + phys_max' <- replicateM channels' $ getInt 8 + dig_min' <- replicateM channels' $ getInt 8 + dig_max' <- replicateM channels' $ getInt 8 + prefilter' <- replicateM channels' $ getString 80 + samples' <- replicateM channels' $ getInt 8 + reserved' <- replicateM channels' $ getString 32 + data_' <- readData (convert phys_min' phys_max' dig_min' dig_max') num_records' (zip [1..] samples') + return $ Just $ BDF { + id_ = BS.head id_' + , type_ = type_' + , subject = subject' + , recording = recording' + , date = date' + , time = time' + , head_bytes = head_bytes' + , data_version = data_version' + , num_records = num_records' + , duration = duration' + , channels = channels' + , chan_labels = chan_labels' + , tran_type = tran_type' + , dimensions = dimensions' + , phys_min = phys_min' + , phys_max = phys_max' + , dig_min = dig_min' + , dig_max = dig_max' + , prefilter = prefilter' + , samples = samples' + , reserved = reserved' + , data_ = data_' + } + +loadBDF :: FilePath -> IO (Maybe (M.Multichannel Double)) +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) + _ -> do + putStrLn "File not read" + return Nothing + when (not (BS.null bs')) $ do + putStrLn "data remaining..." + return m + +{- +getDPConversion :: BDF -> [Word32 -> Double] +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 = ((++ "\n") . concat . intersperse " " . map show) + +main = do + bs <- BS.readFile "VivianMeditation.bdf" + (bdf,bs') <- runStateT readBDF bs + case bdf of + (Just b) -> do + let conv = getDPConversion b + let b' = map (zipWith ($) conv) (data_ b) + putStrLn $ show $ length $ head b' + writeFile "VivianMeditation.txt" $ concat $ map writeLine b' + _ -> putStrLn "File not read" + when (not (BS.null bs')) $ + do + putStrLn "data remaining..." + return () + return () +-} +
lib/Numeric/Signal/Internal.hs view
@@ -29,6 +29,7 @@ import Numeric.LinearAlgebra.Algorithms +import qualified Numeric.GSL.Fourier as F import Foreign import Complex import Foreign.C.Types@@ -49,7 +50,8 @@ ----------------------------------------------------------------------------- instance Convolvable (Vector Double) where- convolve = convolve_vector_double+ convolve x y = fst $ fromComplex $ F.ifft $ (F.fft (complex x) * F.fft (complex y))+-- convolve = convolve_vector_double convolve_vector_double c a = unsafePerformIO $ do r <- createVector (dim a)@@ -61,7 +63,8 @@ ----------------------------------------------------------------------------- instance Convolvable (Vector (Complex Double)) where- convolve = convolve_vector_complex+ convolve x y = F.ifft $ (F.fft x * F.fft y)+-- convolve = convolve_vector_complex convolve_vector_complex c a = unsafePerformIO $ do r <- createVector (dim a)
lib/Numeric/Signal/Multichannel.hs view
@@ -17,6 +17,9 @@ module Numeric.Signal.Multichannel ( Multichannel,+ fromList,+ sampling_rate,precision,channels,+ getChannel,getChannels, mapConcurrently ) where @@ -29,13 +32,13 @@ import qualified Data.Array.IArray as I import Control.Concurrent-import Control.Concurrent.MVar+--import Control.Concurrent.MVar import System.IO.Unsafe(unsafePerformIO) --import qualified Data.List as L -import Data.Packed.Vector+import Data.Packed.Vector hiding(fromList) --import Data.Packed(Container(..)) import Foreign.Storable@@ -49,6 +52,7 @@ ----------------------------------------------------------------------------- +-- | data type with multiple channels data Multichannel a = MC { _sampling_rate :: Int -- ^ sampling rate , _precision :: Int -- ^ bits of precision@@ -59,9 +63,44 @@ ----------------------------------------------------------------------------- +-- | 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+ in MC s p c (dim $ head d) (I.listArray (1,c) d) +-- | the sampling rate+sampling_rate :: Multichannel a -> Int+sampling_rate = _sampling_rate++-- | the bits of precision+precision :: Multichannel a -> Int+precision = _precision++-- | the number of channels+channels :: Multichannel a -> Int+channels = _channels++-- | the length, in samples+samples :: Multichannel a -> Int+samples = _length++-- | extract one channel+getChannel :: Int -> Multichannel a -> Vector a+getChannel c d = (_data d) I.! c++-- | extract all channels+getChannels :: Multichannel a -> [Vector a]+getChannels d = I.elems $ _data d++-----------------------------------------------------------------------------+ -- | map a function executed concurrently-mapConcurrently :: Storable b => Multichannel a -- ^ input data+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