module Numeric.Signal.Multichannel (
Multichannel,readMultichannel,writeMultichannel,
createMultichannel,
sampling_rate,precision,channels,samples,
detrended,filtered,
getChannel,getChannels,
toMatrix,
mapConcurrently,
detrend,filter,
slice,
histograms,
entropy_delta_phase,mi_phase
) where
import qualified Numeric.Signal as S
import qualified Data.Array.IArray as I
import Data.Ix
import Control.Concurrent
import System.IO.Unsafe(unsafePerformIO)
import Data.Binary
import Data.Maybe
import Foreign.Storable
import Numeric.LinearAlgebra
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(replicateM)
data Multichannel a = MC {
_sampling_rate :: Int
, _precision :: Int
, _channels :: Int
, _length :: Int
, _detrended :: Bool
, _filtered :: Maybe (Int,Int)
, _data :: I.Array Int (Vector a)
}
instance Binary (Multichannel Double) 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) = (minElement v,maxElement v)
v' = mapVector (\x -> round $ (x mi)/(ma mi) * (fromIntegral (maxBound :: Word64))) v
in (mi,ma,v' :: Vector Word64)
get = do
s <- get
p <- get
c <- get
l <- get
de <- get
f <- get
(d :: I.Array Int (Double,Double,Vector Word64)) <- get
return $! (MC s p c l de f (seq d (fmap convert) d))
where convert (mi,ma,v) = mapVector (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word64)) * (ma mi) + mi) v
instance Binary (Multichannel Float) 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) = (minElement v,maxElement v)
v' = mapVector (\x -> round $ (x mi)/(ma mi) * (fromIntegral (maxBound :: Word64))) v
in (mi,ma,v' :: Vector Word64)
get = do
s <- get
p <- get
c <- get
l <- get
de <- get
f <- get
(d :: I.Array Int (Float,Float,Vector Word32)) <- get
return $! (MC s p c l de f (seq d (fmap convert) d))
where convert (mi,ma,v) = mapVector (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word32)) * (ma mi) + mi) v
instance Binary (Multichannel (Complex Double)) 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 ((\(r,i) -> (convert r, convert i)) . fromComplex) d
where convert v = let (mi,ma) = (minElement v,maxElement v)
v' = mapVector (\x -> round $ (x mi)/(ma mi) * (fromIntegral (maxBound :: Word64))) v
in (mi,ma,v' :: Vector Word64)
get = do
s <- get
p <- get
c <- get
l <- get
de <- get
f <- get
(d :: I.Array Int ((Double,Double,Vector Word64),(Double,Double,Vector Word64))) <- get
return $! (MC s p c l de f (seq d (fmap (\(r,i) -> toComplex (convert r,convert i)) d)))
where convert (mi,ma,v) = mapVector (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word64)) * (ma mi) + mi) v
instance Binary (Multichannel (Complex Float)) 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 ((\(r,i) -> (convert r, convert i)) . fromComplex) d
where convert v = let (mi,ma) = (minElement v,maxElement 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 :: I.Array Int ((Float,Float,Vector Word32),(Float,Float,Vector Word32))) <- get
return $! (MC s p c l de f (seq d (fmap (\(r,i) -> toComplex (convert r,convert i)) d)))
where convert (mi,ma,v) = mapVector (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word32)) * (ma mi) + mi) v
readMultichannel :: (Binary (Multichannel a)) => FilePath -> IO (Multichannel a)
readMultichannel = decodeFile
writeMultichannel :: (Binary (Multichannel a)) => FilePath -> Multichannel a -> IO ()
writeMultichannel = encodeFile
createMultichannel :: Storable a
=> Int
-> Int
-> [Vector a]
-> Multichannel a
createMultichannel s p d = let c = length d
in MC s p c (dim $ head d) False Nothing (I.listArray (1,c) d)
sampling_rate :: Multichannel a -> Int
sampling_rate = _sampling_rate
precision :: Multichannel a -> Int
precision = _precision
channels :: Multichannel a -> Int
channels = _channels
samples :: Multichannel a -> Int
samples = _length
getChannel :: Int -> Multichannel a -> Vector a
getChannel c d = (_data d) I.! c
getChannels :: Multichannel a -> I.Array Int (Vector a)
getChannels d = _data d
toMatrix :: Element a => Multichannel a -> Matrix a
toMatrix = fromRows . I.elems . _data
detrended :: Multichannel a -> Bool
detrended = _detrended
filtered :: Multichannel a -> Maybe (Int,Int)
filtered = _filtered
mapArrayConcurrently :: Ix i => (a -> b)
-> I.Array i a
-> (I.Array i b)
mapArrayConcurrently f d = unsafePerformIO $ do
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 b vectors
where applyFunction f' (m,(j,e)) = putMVar m (j,f' e)
mapConcurrently :: Storable b
=> (Vector a -> Vector b)
-> Multichannel a
-> Multichannel b
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'
mapMC :: Storable b
=> (Vector a -> Vector b)
-> Multichannel a
-> Multichannel b
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 :: Int -> Multichannel Double -> Multichannel Double
detrend w m = let m' = mapConcurrently (S.detrend w) m
in m' { _detrended = True }
filter :: (S.Filterable a, Double ~ DoubleOf a, Container Vector (Complex a), Convert (Complex a)) =>
(Int,Int) -> Multichannel a -> Multichannel a
filter pb m = let m' = mapConcurrently (S.broadband_filter (_sampling_rate m) pb) m
in m' { _filtered = Just pb }
slice :: Storable a
=> Int
-> Int
-> Multichannel a
-> Multichannel a
slice j w m = let m' = mapConcurrently (subVector j w) m
in m' { _length = w }
histograms :: (S.Filterable a, Double ~ DoubleOf a) =>
I.Array Int (Vector a)
-> Int -> (Double,Double)
-> Int -> Int -> (Double,Double) -> (Double,Double)
-> (I.Array Int H.Histogram,I.Array (Int,Int) H2.Histogram2D)
histograms d' b (l,u) bx by (lx,ux) (ly,uy)
= let d = fmap double d'
(bl,bu) = I.bounds d
br = ((bl,bl),(bu,bu))
histarray = mapArrayConcurrently (H.fromLimits b (l,u)) d
pairs = I.array br $ map (\(m,n) -> ((m,n),(d I.! m,d I.! n))) (range br)
hist2array = mapArrayConcurrently (\(x,y) -> (H2.addVector (H2.emptyLimits bx by (lx,ux) (ly,uy)) x y)) pairs
in (histarray,hist2array)
entropy_delta_phase :: (S.Filterable a, Double ~ DoubleOf a) =>
Multichannel a
-> Matrix Double
entropy_delta_phase m = let d = _data m
c = _channels m
b = ((1,1),(c,c))
r = I.range b
diff = I.listArray b (map (\i@(x,y) -> (i,if x <= y then Just (double $ (d I.! y)(d I.! x)) else Nothing)) r) :: I.Array (Int,Int) ((Int,Int),Maybe (Vector Double))
h = mapArrayConcurrently (maybe Nothing (\di -> Just $ H.fromLimits 128 ((2)*pi,2*pi) di)) (fmap snd diff)
ent = mapArrayConcurrently (\(i,difvec) -> case difvec of
Nothing -> 0 :: Double
Just da -> SI.entropy (fromJust (h I.! i)) da) diff
in fromArray2D ent
mi_phase :: (S.Filterable a, Double ~ DoubleOf a) =>
Multichannel a
-> Matrix Double
mi_phase m = let d = _data m
(histarray,hist2array) = histograms d 128 (pi,pi) 128 128 (pi,pi) (pi,pi)
indhist = I.listArray (I.bounds hist2array) (I.assocs hist2array)
mi = mapArrayConcurrently (doMI histarray (fmap double d)) indhist
in fromArray2D mi
where doMI histarray d ((x,y),h2)
| x <= y = SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,d I.! y)
| otherwise = 0