hsignal 0.1.3.6 → 0.2.1.0
raw patch · 6 files changed
+249/−26 lines, 6 filesdep ~hmatrixdep ~hmatrix-gsl-statsdep ~hstatistics
Dependency ranges changed: hmatrix, hmatrix-gsl-stats, hstatistics
Files
- CHANGES +13/−0
- hsignal.cabal +5/−4
- lib/Numeric/Signal.hs +3/−1
- lib/Numeric/Signal/Internal.hs +3/−0
- lib/Numeric/Signal/Multichannel.hs +107/−21
- lib/Numeric/Signal/Noise.hs +118/−0
CHANGES view
@@ -73,3 +73,16 @@ cleaned up mi_phase generate H(X;X) for reference value +0.2.0.0:+ added instances to store complex values+ Multichannel must now be [Complex] Float/Double+ Multichannel Double is now stored to greater accuracy+ breaks file compatability - can be loaded as Multichannel Float+ added histograms to Multichannel++0.2.0.1:+ added Filterable instance for Complex Float++0.2.1.0:+ added Noise.hs: pinkNoise, spatialNoise, powerNoise+ fixed bug in detrend that occurs if a zero window occurs at end
hsignal.cabal view
@@ -1,5 +1,5 @@ Name: hsignal-Version: 0.1.3.6+Version: 0.2.1.0 License: BSD3 License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010, 2011@@ -37,15 +37,16 @@ mtl > 2, array, bytestring, storable-complex, binary,- hmatrix >= 0.10.0.0,- hmatrix-gsl-stats >= 0.1.2.9,- hstatistics >= 0.2.2.7+ hmatrix >= 0.10.0.4,+ hmatrix-gsl-stats >= 0.1.2.10,+ hstatistics >= 0.2.2.10 Extensions: ForeignFunctionInterface hs-source-dirs: lib Exposed-modules: Numeric.Signal Numeric.Signal.Multichannel+ Numeric.Signal.Noise Numeric.Signal.EEG other-modules: Numeric.Signal.Internal Numeric.Signal.EEG.BDF
lib/Numeric/Signal.hs view
@@ -222,7 +222,9 @@ -> Vector Double -- ^ data to be detrended -> Vector Double -- ^ detrended data detrend w v = let windows = dim v `div` w- ws = takesV ((replicate windows w) ++ [dim v - (windows * w)]) v+ re = dim v - (windows * w)+ re' = if re == 0 then [] else [re]+ ws = takesV ((replicate windows w) ++ re') v ds = map detrend' ws windows' = (dim v - (w `div` 2)) `div` w ws' = takesV (((w `div` 2):(replicate windows' w)) ++ [dim v - (w `div` 2) - (windows' * w)]) v
lib/Numeric/Signal/Internal.hs view
@@ -117,6 +117,9 @@ foreign import ccall "signal-aux.h vector_complex_convolve" signal_vector_complex_convolve :: CInt -> PC -> CInt -> PC -> CInt -> PC -> IO CInt +instance Convolvable (Vector (Complex Float)) where+ convolve x y = single $ F.ifft $ (F.fft (double x) * F.fft (double y))+ ----------------------------------------------------------------------------- instance Filterable Double where
lib/Numeric/Signal/Multichannel.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE UndecidableInstances,+ FlexibleInstances, FlexibleContexts, TypeFamilies, ScopedTypeVariables #-}@@ -18,6 +19,9 @@ -- ----------------------------------------------------------------------------- +-- IncoherentInstances,++ module Numeric.Signal.Multichannel ( Multichannel,readMultichannel,writeMultichannel, createMultichannel,@@ -28,6 +32,7 @@ mapConcurrently, detrend,filter, slice,+ histograms, mi_phase ) where @@ -98,10 +103,7 @@ ----------------------------------------------------------------------------- -instance (Binary a, Storable a, - Ord a, RealFrac a,- Container Vector a,- Product a) => Binary (Multichannel a) where+instance Binary (Multichannel Double) where put (MC s p c l de f d) = do put s put p@@ -111,8 +113,8 @@ 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 :: Word32))) v- in (mi,ma,v' :: Vector Word32) + v' = mapVector (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word64))) v+ in (mi,ma,v' :: Vector Word64) get = do s <- get@@ -121,22 +123,92 @@ l <- get de <- get f <- get+ (d :: I.Array Int (a,a,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 (a,a,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 ((a,a,Vector Word64),(a,a,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 ((a,a,Vector Word32),(a,a,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 a, Storable a, - Ord a, RealFrac a,- Container Vector a,- Product a) => FilePath -> IO (Multichannel a)+readMultichannel :: (Binary (Multichannel a)) => FilePath -> IO (Multichannel a) readMultichannel = decodeFile -writeMultichannel :: (Binary a, Storable a, - Ord a, RealFrac a,- Container Vector a,- Product a) => FilePath -> Multichannel a -> IO ()+writeMultichannel :: (Binary (Multichannel a)) => FilePath -> Multichannel a -> IO () writeMultichannel = encodeFile -----------------------------------------------------------------------------@@ -255,17 +327,31 @@ ----------------------------------------------------------------------------- +-- | calculate histograms+histograms :: (S.Filterable a, Double ~ DoubleOf a) =>+ I.Array Int (Vector a)+ -> Int -> (Double,Double) + -> Int -> Int -> (Double,Double) -> (Double,Double) -- ^ bins and ranges+ -> (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)++-----------------------------------------------------------------------------+ -- | calculate the mutual information of the phase between pairs of channels (fills upper half of matrix) mi_phase :: (S.Filterable a, Double ~ DoubleOf a) => Multichannel a -- ^ input data -> Matrix Double-mi_phase m = let d = fmap double $ _data m- (bl,bu) = I.bounds d- br = ((bl,bl),(bu,bu))- histarray = mapArrayConcurrently (H.fromLimits 128 (-pi,pi)) d- pairs = I.array br $ map (\(a,b) -> ((a,b),((a,b),d I.! a,d I.! b))) (range br)- hist2array = mapArrayConcurrently (\(j,x,y) -> (j,H2.addVector (H2.emptyLimits 128 128 (-pi,pi) (-pi,pi)) x y)) pairs- mi = mapArrayConcurrently (doMI histarray d) hist2array+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)
+ lib/Numeric/Signal/Noise.hs view
@@ -0,0 +1,118 @@+--{-# LANGUAGE UndecidableInstances,+-- FlexibleInstances,+-- FlexibleContexts,+-- TypeFamilies,+-- ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Signal.Noise+-- Copyright : (c) Alexander Vivian Hugh McPhail 2010+-- License : GPL-style+--+-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com+-- Stability : provisional+-- Portability : uses Concurrency+--+-- Noise generation functions+--+-----------------------------------------------------------------------------++-- IncoherentInstances,+++module Numeric.Signal.Noise (+ pinkNoise+ , spatialNoise+ , powerNoise+ ) where++-----------------------------------------------------------------------------++--import qualified Numeric.Signal as S++--import Complex++--import qualified Data.Array.IArray as I+--import Data.Ix++--import Data.Word++--import System.IO.Unsafe(unsafePerformIO)++--import qualified Data.List as L++--import Data.Binary++--import Foreign.Storable+++import Numeric.Container+import Numeric.LinearAlgebra++import qualified Numeric.GSL.Fourier as F++--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)++-------------------------------------------------------------------++-- | The method is briefly descirbed in Lennon, J.L. "Red-shifts and red+-- herrings in geographical ecology", Ecography, Vol. 23, p101-113 (2000)+--+-- Matlab version Written by Jon Yearsley 1 May 2004+-- j.yearsley@macaulay.ac.uk+--+-- Creates 1/f scale invariant spatial noise+spatialNoise :: Double -- ^ β: spectral distribution+ -- 0: White noise+ -- -1: Pink noise+ -- -2: Brownian noise+ -> Int -> Int -- ^ matrix dimensions+ -> Int -- ^ random seed+ -> Matrix Double+spatialNoise b r' c' s = let c = fromIntegral c'+ r = fromIntegral r'+ pre_x = linspace c' (0::Double,c-1)+ post_x = linspace c' (c,1)+ freq_x = mapVector (/c) $ join [pre_x,post_x]+ u = fromRows (replicate (2*r') freq_x)+ pre_y = linspace r' (0::Double,r-1)+ post_y = linspace r' (r,1)+ freq_y = mapVector (/c) $ join [pre_y,post_y]+ v = fromColumns (replicate (2*c') freq_y)+ s_f = liftMatrix (mapVector (**(b/2))) ((u**2) + (v**2))+ s_f' = liftMatrix (mapVector (\x -> if isInfinite x then 0 else x)) s_f+ phi = reshape (2*c') (randomVector s Uniform (4*r'*c'))+ in subMatrix (1,1) (r',c') $ fst $ fromComplex $ fromRows $ map F.ifft $ toRows $ ((complex $ s_f'**0.5) * (toComplex (cos(2*pi*phi),sin(2*pi*phi))))++-- | 1/f scale invariant noise+pinkNoise :: + Double -- ^ β: spectral distribution+ -- 0: White noise+ -- -1: Pink noise+ -- -2: Brownian (red) noise+ -> Int -- ^ samples+ -> Int -- ^ random seed+ -> Vector Double +pinkNoise b s r = let pre = linspace s (0::Double,fromIntegral (s-1))+ post = linspace s (fromIntegral s,1)+ freq = join [pre/(fromIntegral s),post/(fromIntegral s)]+ s_f = mapVector (**(b/2)) (freq**2) + s_f' = mapVector (\x -> if isInfinite x then 0 else x) s_f+ phi = randomVector r Uniform (2*s)+ in subVector 0 s $ fst $ fromComplex $ F.ifft ((complex $ s_f'**0.5) * (toComplex (cos(2*pi*phi),sin(2*pi*phi))))++-- | generate noise from a power spectrum+powerNoise :: Vector Double -- ^ the power spectrum+ -> Int -- ^ random seed+ -> Vector Double+powerNoise psd r = let ln = dim psd+ freq = join [fromList [0],psd, (fromList . reverse . tail . toList) psd]+ phi = randomVector r Uniform (2*ln)+ in (fromIntegral ln) * (subVector 0 (ln-1) $ fst $ fromComplex $ F.ifft ((complex $ freq) * (toComplex (cos(2*pi*phi),sin(2*pi*phi)))))