diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -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 diff --git a/hsignal.cabal b/hsignal.cabal --- a/hsignal.cabal +++ b/hsignal.cabal @@ -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 diff --git a/lib/Numeric/Signal.hs b/lib/Numeric/Signal.hs --- a/lib/Numeric/Signal.hs +++ b/lib/Numeric/Signal.hs @@ -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 diff --git a/lib/Numeric/Signal/Internal.hs b/lib/Numeric/Signal/Internal.hs --- a/lib/Numeric/Signal/Internal.hs +++ b/lib/Numeric/Signal/Internal.hs @@ -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 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 @@ {-# 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) diff --git a/lib/Numeric/Signal/Noise.hs b/lib/Numeric/Signal/Noise.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Signal/Noise.hs @@ -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 gmail 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)))))