diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -30,3 +30,13 @@ fixed signal-aux.c filter scaling error fixed hilbert changed Multichannel.fromList to createMultichannel to avoid name conflicts + +0.1.2.2: + changed Multichannel read/write type signature to reflect hmatrix change + removed haskell98 dependency + added resample/resize/deriv + removed attempt to compare against sampling of histPDF (requires random number calls -> seed/IO) + added Multichannel.toMatrix + fixed scale in pwelch + reflected changes in hmatrix exports + diff --git a/INSTALL b/INSTALL --- a/INSTALL +++ b/INSTALL @@ -1,35 +1,17 @@ ----------------------------------------------- - A simple signal processing library for Haskell +A signal processing library for Haskell ----------------------------------------------- INSTALLATION -Recommended method (ok in Ubuntu/Debian systems): - $ cabal install hsignal - -INSTALLATION ON WINDOWS ---------------------------------------- - -1) Install a recent ghc (e.g. ghc-6.10.3) - -2) Install cabal-install. A binary for windows can be obtained from: - - http://www.haskell.org/cabal/release/cabal-install-0.6.2/cabal.exe - - Put it somewhere in the path, for instance in c:\ghc\ghc-6.10.3\bin - -3) Download and uncompress hmatrix-x.y.z.tar.gz from Hackage: - - http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hmatrix - -4) Open a terminal, cd to the hmatrix folder, and run - - > cabal install - -5) Download and uncompress hsignal-x.y.z.tar.gz from Hackage: - - http://hackage.haskell.org/cgi-bin/hackage-scripts/package/hsignal +cabal install hsignal -6) Open a terminal, cd to the hsignal folder, and run +OR - > cabal install +tar xzf hsignal-x.y.z.tar.gz +cd hsignal +runhaskell Setup.lhs configure +runhaskell Setup.lhs build +runhaskell Setup.lhs hadock +runhaskell Setup.lhs install diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,2 +1,2 @@ -Copyright Alberto Ruiz 2006-2007 +Copyright A.V.H. McPhail 2010 GPL license diff --git a/hsignal.cabal b/hsignal.cabal --- a/hsignal.cabal +++ b/hsignal.cabal @@ -1,5 +1,5 @@ Name: hsignal -Version: 0.1.2.1 +Version: 0.1.2.2 License: GPL License-file: LICENSE Author: Vivian McPhail @@ -7,10 +7,14 @@ Stability: provisional 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, +Description: + Purely functional interface to signal processing based on hmatrix + . + Provides data types for manipulating EEG data, including reading from BDF data format files + . + When hmatrix is installed with -fvector, the vector type is Data.Vector.Storable + from the vector package. +Category: Math tested-with: GHC ==6.12.1 cabal-version: >=1.8 @@ -23,12 +27,12 @@ library Build-Depends: base >= 3 && < 5, - haskell98, mtl, + mtl, array, bytestring, storable-complex, ghc-binary, hmatrix >= 0.10.0, - hmatrix-gsl-stats >= 0.1.1.5, - hstatistics >= 0.2.0.5 + hmatrix-gsl-stats >= 0.1.2.1, + hstatistics >= 0.2.0.8 Extensions: ForeignFunctionInterface diff --git a/lib/Numeric/Signal.hs b/lib/Numeric/Signal.hs --- a/lib/Numeric/Signal.hs +++ b/lib/Numeric/Signal.hs @@ -14,13 +14,20 @@ ----------------------------------------------------------------------------- module Numeric.Signal ( + -- * Filtering hamming, pwelch, fir,standard_fir,broadband_fir, freqzF,freqzN, filter,broadband_filter, + -- * Analytic Signal analytic_signal,analytic_power,analytic_phase, - detrend + unwrap, + -- * Preprocessing + detrend, + downsample,resize, + -- * Utility functions + deriv ) where ----------------------------------------------------------------------------- @@ -38,7 +45,8 @@ import Numeric.GSL.Vector -import Numeric.LinearAlgebra.Algorithms +--import Numeric.LinearAlgebra.Algorithms +import Numeric.LinearAlgebra.Linear --import Numeric.LinearAlgebra.Interface import qualified Numeric.GSL.Fourier as F @@ -76,9 +84,9 @@ -> (Vector Double,Vector Double) -- ^ (frequency index,power density) pwelch s w v = let w' = max s w -- make window at least sampling rate r = S.pwelch w' v - sd = recip $ (fromIntegral s)/2 + sd = (fromIntegral s)/2 -- scale for sampling rate - r' = scale sd r + r' = scale (recip sd) r f = linspace ((w `div` 2) + 1) (0,sd) in (f,r') @@ -221,5 +229,28 @@ t = linspace ln (1.0,fromIntegral ln) (c0,c1,_,_,_,_) = linear t x in x - (scale c1 t + scalar c0) + +----------------------------------------------------------------------------- + +-- | take one sample from every n samples in the original +downsample :: Int -> Vector Double -> Vector Double +downsample = S.downsample + + +-- | resize the vector to length n by resampling +resize :: Int -> Vector Double -> Vector Double +resize n v = downsample (dim v `div` n) v + +----------------------------------------------------------------------------- + +-- | the difference between consecutive elements of a vector +deriv :: Vector Double -> Vector Double +deriv = S.deriv + +----------------------------------------------------------------------------- + +-- | unwrap the phase of signal (input expected to be within (-pi,pi) +unwrap :: Vector Double -> Vector Double +unwrap = S.unwrap ----------------------------------------------------------------------------- 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 @@ -18,20 +18,24 @@ hamming, filter, freqz, - hilbert, pwelch, - complex_power + hilbert, + complex_power, + downsample, + deriv, + unwrap ) where import Data.Packed.Development(createVector,vec,app1,app2,app3,app4) import Data.Packed.Vector import Data.Packed(Container(..)) -import Numeric.LinearAlgebra.Algorithms +--import Numeric.LinearAlgebra.Algorithms +import Numeric.LinearAlgebra.Linear import qualified Numeric.GSL.Fourier as F import Foreign -import Complex +import Data.Complex import Foreign.C.Types import Prelude hiding(filter) @@ -157,7 +161,7 @@ ----------------------------------------------------------------------------- --- | the complex power : real $ v * (conj c) +-- | the complex power : real $ v * (conj v) complex_power :: Vector (Complex Double) -- ^ input -> Vector Double -- ^ output complex_power v = unsafePerformIO $ do @@ -168,3 +172,39 @@ foreign import ccall "signal-aux.h complex_power" signal_complex_power :: CInt -> PC -> CInt -> PD -> IO CInt ----------------------------------------------------------------------------- + +-- | resample, take one sample every n samples in the original +downsample :: Int -> Vector Double -> Vector Double +downsample n v = unsafePerformIO $ do + r <- createVector (dim v `div` n) + app2 (signal_downsample $ fromIntegral n) vec v vec r "downsample" + return r + +foreign import ccall "signal-aux.h downsample" signal_downsample :: CInt -> CInt -> PD -> CInt -> PD -> IO CInt + +----------------------------------------------------------------------------- + +-- | the difference between consecutive elements of a vector +deriv :: Vector Double -> Vector Double +deriv v = unsafePerformIO $ do + r <- createVector (dim v - 1) + app2 (signal_deriv) vec v vec r "diff" + return r + +foreign import ccall "signal-aux.h vector_deriv" signal_deriv :: CInt -> PD -> CInt -> PD -> IO CInt + +----------------------------------------------------------------------------- + +-- | unwrap the phase of signal (input expected to be within (-pi,pi) +unwrap :: Vector Double -> Vector Double +unwrap v = unsafePerformIO $ do + r <- createVector $ dim v + app2 signal_unwrap vec v vec r "unwrap" + return r + +foreign import ccall "signal-aux.h unwrap" signal_unwrap :: CInt -> PD -> CInt -> PD -> IO CInt + +----------------------------------------------------------------------------- + + + 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 @@ -12,7 +12,7 @@ -- -- Signal processing functions, multichannel datatype -- --- link with '-threaded' and run with +RTS Nn, where n is the number of CPUs +-- link with '-threaded' and run with +RTS -Nn, where n is the number of CPUs -- ----------------------------------------------------------------------------- @@ -22,6 +22,7 @@ sampling_rate,precision,channels,samples, detrended,filtered, getChannel,getChannels, + toMatrix, mapConcurrently, detrend,filter, slice, @@ -53,6 +54,7 @@ import Foreign.Storable --import Numeric.GSL.Vector +import Numeric.LinearAlgebra.Linear --import Numeric.LinearAlgebra.Algorithms --import qualified Numeric.GSL.Fourier as F @@ -98,21 +100,9 @@ ----------------------------------------------------------------------------- -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 + Vectors Vector a) => Binary (Multichannel a) where put (MC s p c l de f d) = do put s put p @@ -121,7 +111,7 @@ put de put f put $! fmap convert d - where convert v = let (mi,ma) = (vMin v,vMax v) + where convert v = let (mi,ma) = (vectorMin v,vectorMax v) v' = mapVector (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word32))) v in (mi,ma,(v' :: Vector Word32)) @@ -142,20 +132,21 @@ readMultichannel :: (Binary a, Storable a, Ord a, RealFrac a, - Vectors a) => FilePath -> IO (Multichannel a) + Vectors Vector a) => FilePath -> IO (Multichannel a) readMultichannel = decodeFile writeMultichannel :: (Binary a, Storable a, Ord a, RealFrac a, - Vectors a) => FilePath -> Multichannel a -> IO () + Vectors Vector a) => FilePath -> Multichannel a -> IO () writeMultichannel = encodeFile ----------------------------------------------------------------------------- -- | create a multichannel data type -createMultichannel :: Storable a => Int -- ^ sampling rate +createMultichannel :: Storable a + => Int -- ^ sampling rate -> Int -- ^ bits of precision - -> [Vector a] -- ^ data + -> [Vector a] -- ^ data -> Multichannel a -- ^ datatype createMultichannel s p d = let c = length d in MC s p c (dim $ head d) False Nothing (I.listArray (1,c) d) @@ -184,6 +175,10 @@ getChannels :: Multichannel a -> I.Array Int (Vector a) getChannels d = _data d +-- | convert the data to a matrix with channels as rows +toMatrix :: Element a => Multichannel a -> Matrix a +toMatrix = fromRows . I.elems . _data + -- | was the data detrended? detrended :: Multichannel a -> Bool detrended = _detrended @@ -219,14 +214,16 @@ -} -- | 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 :: 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 +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 @@ -248,7 +245,8 @@ ----------------------------------------------------------------------------- -- | extract a slice of the data -slice :: Storable a => Int -- ^ starting sample number +slice :: Storable a + => Int -- ^ starting sample number -> Int -- ^ length -> Multichannel a -> Multichannel a @@ -257,25 +255,18 @@ ----------------------------------------------------------------------------- --- | calculate the mutual information of the phase between pairs of channels --- | the lower half of the matrix displays mutual information between pairs of channels --- | the upper half of the matrix displays mutual information between the data (row) and a sample --- | from the same distribution as the column +-- | calculate the mutual information of the phase between pairs of channels (fills upper half of matrix) mi_phase :: Multichannel Double -- ^ input data -> Matrix Double mi_phase m = let d = _data m - s = samples m histarray = mapArrayConcurrently (H.fromLimits 128 (-pi,pi)) d c = channels m pairs = I.array ((1,1),(c,c)) $ map (\(a,b) -> ((a,b),((a,b),d I.! a,d I.! b))) (range ((1,1),(c,c))) hist2array = mapArrayConcurrently (\(j,x,y) -> (j,H2.addVector (H2.emptyLimits 128 128 (-pi,pi) (-pi,pi)) x y)) pairs - mi = mapArrayConcurrently (doMI histarray d s) hist2array + mi = mapArrayConcurrently (doMI histarray d) hist2array in fromArray2D mi - where doMI histarray d s ((x,y),h2) + 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 = let hypdf = H.fromHistogram (histarray I.! y) - ys = replicate s (H.sample hypdf) - dy' = fromList ys - in SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,dy') + | otherwise = 0 ----------------------------------------------------------------------------- diff --git a/lib/Numeric/Signal/signal-aux.c b/lib/Numeric/Signal/signal-aux.c --- a/lib/Numeric/Signal/signal-aux.c +++ b/lib/Numeric/Signal/signal-aux.c @@ -137,7 +137,8 @@ else if (j == (rs-1)) r[j] = s[j]; else r[j] = s[j] + s[fs-j+1]; - r[j] /= (num_windows * w / 2); + r[j] /= num_windows; + r[j] = sqrt(r[j]); } gsl_fft_complex_wavetable_free (wavetable); gsl_fft_complex_workspace_free (workspace); @@ -174,6 +175,50 @@ for (i = 0; i < cs; i++) r[i] = c[i].dat[0]*c[i].dat[0] + c[i].dat[1]*c[i].dat[1]; + + return 0; +} + +int downsample(int n, int xs, const double* x, int rs, double* r) +{ + if (rs != xs/n) return 2000; // BAD_SIZE + + int i; + + for (i = 0; i < rs; i++) + r[i] = x[i*n]; + + return 0; +} + +int vector_deriv(int xs, const double* x, int rs, double* r) +{ + if (rs != xs - 1) return 2000; // BAD_SIZE + + int i; + + for (i = 0; i < rs; i++) + r[i] = x[i+1] - x[i]; + + return 0; +} + +int unwrap(int xs, const double* x, int rs, double* r) +{ + if (rs != xs) return 2000; // BAD_SIZE + + int i; + + r[0] = x[0]; + + int j = 0; + + for (i = 1; i < rs; i++) { + if (x[i] < x[i-1]) { + j += 1; + } + r[i] = x[i] + j*2*M_PI; + } return 0; }