hsignal-0.1.1.1: lib/Numeric/Signal/Multichannel.hs
{-# OPTIONS_GHC -fglasgow-exts #-}
{-# OPTIONS_GHC -XUndecidableInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module : Numeric.Signal.Multichannel
-- Copyright : (c) Alexander Vivian Hugh McPhail 2010
-- License : GPL-style
--
-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com
-- Stability : provisional
-- Portability : uses Concurrency
--
-- Signal processing functions, multichannel datatype
--
-- link with '-threaded' and run with +RTS Nn, where n is the number of CPUs
--
-----------------------------------------------------------------------------
module Numeric.Signal.Multichannel (
Multichannel,readMultichannel,writeMultichannel,
fromList,
sampling_rate,precision,channels,samples,
detrended,filtered,
getChannel,getChannels,
mapConcurrently,
detrend,filter,
slice
) where
-----------------------------------------------------------------------------
import qualified Numeric.Signal as S
--import Complex
import qualified Data.Array.IArray as I
import Control.Concurrent
--import Control.Concurrent.MVar
import System.IO.Unsafe(unsafePerformIO)
--import qualified Data.List as L
import Data.Packed.Vector hiding(fromList)
--import Data.Packed(Container(..))
import Data.Binary
import Foreign.Storable
--import Numeric.GSL.Vector
--import Numeric.LinearAlgebra.Algorithms
--import qualified Numeric.GSL.Fourier as F
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
, _detrended :: Bool -- ^ was the data detrended?
, _filtered :: Maybe (Int,Int) -- ^ if filtered the passband
, _data :: I.Array Int (Vector a) -- ^ data
}
-----------------------------------------------------------------------------
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
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 :: (Binary a, Storable a,
Ord a, RealFrac a,
Vectors a) => FilePath -> IO (Multichannel a)
readMultichannel = decodeFile
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 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
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 -> 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
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 ((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 }
-----------------------------------------------------------------------------