hsignal 0.1.2.1 → 0.1.2.2
raw patch · 8 files changed
+184/−81 lines, 8 filesdep −haskell98dep ~hmatrix-gsl-statsdep ~hstatistics
Dependencies removed: haskell98
Dependency ranges changed: hmatrix-gsl-stats, hstatistics
Files
- CHANGES +10/−0
- INSTALL +9/−27
- LICENSE +1/−1
- hsignal.cabal +12/−8
- lib/Numeric/Signal.hs +35/−4
- lib/Numeric/Signal/Internal.hs +45/−5
- lib/Numeric/Signal/Multichannel.hs +26/−35
- lib/Numeric/Signal/signal-aux.c +46/−1
CHANGES view
@@ -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+
INSTALL view
@@ -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
LICENSE view
@@ -1,2 +1,2 @@-Copyright Alberto Ruiz 2006-2007+Copyright A.V.H. McPhail 2010 GPL license
hsignal.cabal view
@@ -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
lib/Numeric/Signal.hs view
@@ -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 -----------------------------------------------------------------------------
lib/Numeric/Signal/Internal.hs view
@@ -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++-----------------------------------------------------------------------------+++
lib/Numeric/Signal/Multichannel.hs view
@@ -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 -----------------------------------------------------------------------------
lib/Numeric/Signal/signal-aux.c view
@@ -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; }