hsignal 0.1.2.5 → 0.1.3

raw patch · 6 files changed

+340/−113 lines, 6 filesdep +base-unicode-symbolsPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: base-unicode-symbols

API changes (from Hackage documentation)

Files

• CHANGES +3/−0
• hsignal.cabal +5/−4
• lib/Numeric/Signal.hs +77/−64
• lib/Numeric/Signal/Internal.hs +172/−38
• lib/Numeric/Signal/Multichannel.hs +5/−3
• lib/Numeric/Signal/signal-aux.c +78/−4

@@ -48,3 +48,6 @@  0.1.2.5: 		hmatrix 0.10.0.0 release++0.1.3:+		Add Float versions of some functions

@@ -1,5 +1,5 @@ Name:               hsignal-Version:            0.1.2.5+Version:            0.1.3 License:            BSD3 License-file:       LICENSE Author:             Vivian McPhail@@ -36,6 +36,7 @@                         mtl,                          array,                         bytestring, storable-complex, binary,+                        base-unicode-symbols,                         hmatrix >= 0.10.0.0,                         hmatrix-gsl-stats >= 0.1.2.6,                         hstatistics >= 0.2.2.5@@ -78,7 +79,7 @@     extra-libraries:     extra-lib-dirs: -    source-repository head-        type:     darcs-        location: darcs get http://code.haskell.org/hsignal+source-repository head+    type:     darcs+    location: http://code.haskell.org/hsignal 

@@ -1,4 +1,6 @@-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UnicodeSyntax #-} ----------------------------------------------------------------------------- -- | -- Module      :  Numeric.Signal@@ -14,6 +16,8 @@ -----------------------------------------------------------------------------  module Numeric.Signal (+                       S.Convolvable(..),+                       S.Filterable(),                        -- * Filtering                        hamming,                        pwelch,@@ -25,9 +29,9 @@                        unwrap,                        -- * Preprocessing                        detrend,-                       downsample,resize,-                       -- * Utility functions-                       deriv+                       resize,+                       downsample,+                       deriv,                 ) where  -----------------------------------------------------------------------------@@ -37,7 +41,10 @@ import Numeric.GSL.Fitting.Linear   import Data.Complex+import Foreign.Storable() +--import Data.Function.Unicode+ import qualified Data.List as L  --import Data.Packed.Vector@@ -51,27 +58,21 @@ -----------------------------------------------------------------------------  -- | filters the signal-filter :: Vector Double -- ^ zero coefficients-       -> Vector Double -- ^ pole coefficients-       -> Int           -- ^ sampling rate-       -> Vector Double -- ^ input signal-       -> Vector Double -- ^ output signal+filter :: (S.Filterable a) +         ⇒ Vector a   -- ^ zero coefficients+       -> Vector a     -- ^ pole coefficients+       -> Int   -- ^ sampling rate+       -> Vector a     -- ^ input signal+       -> Vector a     -- ^ output signal filter b a s v = let len = dim v                      w = min s len                      start = (negate . fromList . reverse . toList . subVector 0 w) v                      finish = (negate . fromList . reverse . toList . subVector (len-w) w) v                      v' = join [start,v,finish]-                 in subVector s len $ S.filter b a v'+                 in subVector s len $ S.filter_ b a v'  -----------------------------------------------------------------------------                      --- | coefficients of a Hamming window-hamming :: Int           -- ^ length-        -> Vector Double -- ^ the Hamming coeffficents-hamming = S.hamming-------------------------------------------------------------------------------- -- | Welch (1967) power spectrum density using periodogram/FFT method pwelch :: Int            -- ^ sampling rate        -> Int            -- ^ window size@@ -88,9 +89,10 @@ -----------------------------------------------------------------------------  -- | a broadband FIR-broadband_fir :: Int           -- ^ sampling rate+broadband_fir :: (S.Filterable a, Double ~ DoubleOf a, Container Vector (Complex a), Convert (Complex a))  ⇒+                Int           -- ^ sampling rate               -> (Int,Int)     -- ^ (lower,upper) frequency cutoff-              -> Vector Double -- ^ filter coefficients   +              -> Vector a -- ^ filter coefficients    broadband_fir s (l,h) = let o = 501                             ny = (fromIntegral s) / 2.0                             fl = (fromIntegral l) / ny@@ -101,21 +103,23 @@                         in standard_fir o be  -- | a broadband filter-broadband_filter :: Int           -- ^ sampling rate-                 -> (Int,Int)     -- ^ (lower,upper) frequency cutoff-                 -> Vector Double -- ^ input signal-                 -> Vector Double -- ^ output signal-broadband_filter s f v = let b = broadband_fir s f+broadband_filter :: (S.Filterable a, Double ~ DoubleOf a, Container Vector (Complex a), Convert (Complex a)) +                   ⇒ Int        -- ^ sampling rate+                 -> (Int,Int)    -- ^ (lower,upper) frequency cutoff+                 -> Vector a            -- ^ input signal+                 -> Vector a            -- ^ output signal+broadband_filter s f v = let b = S.fromDouble $ broadband_fir s f                          in filter b (scalar 1.0) s v                                  -----------------------------------------------------------------------------  -- | standard FIR filter -- |   FIR filter with grid a power of 2 greater than the order, ramp = grid/16, hamming window-standard_fir :: Int -> [(Double,Double)] -> Vector Double+standard_fir :: (S.Filterable a, Double ~ DoubleOf a, Container Vector (Complex a), Convert (Complex a)) ⇒ +               Int -> [(a,a)] -> Vector a standard_fir o be = let grid  = calc_grid o                         trans = grid `div` 16-                    in fir o be grid trans $ S.hamming (o+1)+                    in fir o be grid trans $ S.hamming_ (o+1)  calc_grid :: Int -> Int calc_grid o = let next_power = ceiling (((log $ fromIntegral o) :: Double) / (log 2.0)) :: Int@@ -123,20 +127,22 @@   -- | produce an FIR filter-fir :: Int               -- ^ order (one less than the length of the filter)-    -> [(Double,Double)] -- ^ band edge frequency, nondecreasing, [0, f1, ..., f(n-1), 1]-                         -- ^ band edge magnitude+fir :: (S.Filterable a+      , Container Vector (Complex a), Convert (Complex a), Double ~ DoubleOf a) ⇒+      Int               -- ^ order (one less than the length of the filter)+    -> [(a,a)] -- ^ band edge frequency, nondecreasing, [0, f1, ..., f(n-1), 1]+                        -- ^ band edge magnitude     -> Int               -- ^ grid spacing     -> Int               -- ^ transition width-    -> Vector Double     -- ^ smoothing window (size is order + 1)-    -> Vector Double     -- ^ the filter coefficients+    -> Vector a     -- ^ smoothing window (size is order + 1)+    -> Vector a     -- ^ the filter coefficients fir o be gn tn w = let mid = o `div` 2                        (f,m) = unzip be-                       f' = diff (((fromIntegral gn) :: Double)/((fromIntegral tn) :: Double)/2.0) f+                       f' = diff (((fromIntegral gn))/((fromIntegral tn))/2.0) f                        m' = interpolate f m f'                        grid = interpolate f' m' $ map (\x -> (fromIntegral x)/(fromIntegral gn)) [0..(gn-1)]                        grid' = map (\x -> x :+ 0) grid-                       (b,_) = fromComplex $ F.ifft $ fromList $ grid' ++ (reverse (drop 1 grid'))+                       b = S.fromDouble $ fst $ fromComplex $ F.ifft $ double $ fromList $ grid' ++ (reverse (drop 1 grid'))                        b' = join [subVector ((dim b)-mid-1) (mid+1) b, subVector 1 (mid+1) b]                     in b' * w @@ -148,24 +154,24 @@     | x > 1.0   = 1.0     | otherwise = x -diff :: Double -> [Double] -> [Double]+diff :: S.Filterable a ⇒ a -> [a] -> [a] diff _ []  = [] diff _ [x] = [x] diff inc (x1:x2:xs)      | x1 == x2     = (floor_zero $ x1-inc):x1:(ceil_one $ x1+inc):(diff inc (L.filter (/= x2) xs))      | otherwise    = x1:(diff inc (x2:xs)) -interpolate :: [Double] -> [Double] -> [Double] -> [Double]+interpolate :: S.Filterable a ⇒ [a] -> [a] -> [a] -> [a] interpolate _ _ []      = [] interpolate x y (xp:xs) = if xp == 1.0                               then ((interpolate'' ((length x)-1) x y xp):(interpolate x y xs))                              else ((interpolate' x y xp):(interpolate x y xs)) -interpolate' :: [Double] -> [Double] -> Double -> Double+interpolate' :: S.Filterable a ⇒ [a] -> [a] -> a -> a interpolate' x y xp = let Just j = L.findIndex (> xp) x                       in (interpolate'' j x y xp) -interpolate'' :: Int -> [Double] -> [Double] -> Double -> Double+interpolate'' :: S.Filterable a ⇒ Int -> [a] -> [a] -> a -> a interpolate'' j x y xp = let x0 = x !! (j-1)                              y0 = y !! (j-1)                              x1 = x !! j@@ -175,19 +181,21 @@ -----------------------------------------------------------------------------  -- | determine the frequency response of a filter, given a vector of frequencies-freqzF :: Vector Double     -- ^ zero coefficients-       -> Vector Double     -- ^ pole coefficients-       -> Int               -- ^ sampling rate   -       -> Vector Double     -- ^ frequencies-       -> Vector Double     -- ^ frequency response+freqzF :: (S.Filterable a, Double ~ DoubleOf a, S.Filterable (DoubleOf a)) ⇒ +         Vector a     -- ^ zero coefficients+       -> Vector a       -- ^ pole coefficients+       -> Int     -- ^ sampling rate   +       -> Vector a       -- ^ frequencies+       -> Vector a       -- ^ frequency response freqzF b a s f = S.freqz b a ((2*pi/(fromIntegral s)) * f)  -- | determine the frequency response of a filter, given a number of points and sampling rate-freqzN :: Vector Double     -- ^ zero coefficients-       -> Vector Double     -- ^ pole coefficients-       -> Int               -- ^ sampling rate-       -> Int               -- ^ number of points-       -> (Vector Double,Vector Double)     -- ^ (frequencies,response)+freqzN :: (S.Filterable a, Enum a, Double ~ DoubleOf a) ⇒ +         Vector a     -- ^ zero coefficients+       -> Vector a       -- ^ pole coefficients+       -> Int     -- ^ sampling rate+       -> Int     -- ^ number of points+       -> (Vector a,Vector a)   -- ^ (frequencies,response) freqzN b a s n = let w' = linspace n (0,((fromIntegral n)-1)/(fromIntegral (2*n)))                      r = S.freqz b a ((2*pi)*w')                      in ((fromIntegral s)*w',r)@@ -199,12 +207,14 @@ analytic_signal = S.hilbert  -- | the power (amplitude^2 = v * (conj c)) of an analytic signal-analytic_power :: Vector (Complex Double) -> Vector Double-analytic_power = S.complex_power+analytic_power :: S.Filterable a ⇒ Vector (Complex Double) -> Vector a+analytic_power = S.complex_power_  -- | the phase of an analytic signal-analytic_phase :: Vector (Complex Double) -> Vector Double-analytic_phase = uncurry arctan2 . fromComplex+analytic_phase :: (S.Filterable a, Container Vector a+                 ,Double ~ DoubleOf a) ⇒ +                 Vector (Complex a) -> Vector a+analytic_phase = (uncurry arctan2) . fromComplex  ----------------------------------------------------------------------------- @@ -226,25 +236,28 @@  ----------------------------------------------------------------------------- --- | 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+resize :: S.Filterable a ⇒ Int -> Vector a -> Vector a+resize n v = S.downsample_ (dim v `div` n) v  ----------------------------------------------------------------------------- --- | the difference between consecutive elements of a vector-deriv :: Vector Double -> Vector Double-deriv = S.deriv+-- | coefficients of a Hamming window+hamming :: S.Filterable a ⇒+          Int           -- ^ length+        -> Vector a -- ^ the Hamming coeffficents+hamming = S.hamming_ ------------------------------------------------------------------------------+-- | resample, take one sample every n samples in the original+downsample :: S.Filterable a ⇒ Int -> Vector a -> Vector a+downsample = S.downsample_ +-- | the difference between consecutive elements of a vector+deriv :: S.Filterable a ⇒ Vector a -> Vector a+deriv = S.deriv_+ -- | unwrap the phase of signal (input expected to be within (-pi,pi)-unwrap :: Vector Double -> Vector Double-unwrap = S.unwrap+unwrap :: S.Filterable a ⇒ Vector a -> Vector a+unwrap = S.unwrap_  -----------------------------------------------------------------------------

@@ -1,4 +1,6 @@-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-} ----------------------------------------------------------------------------- -- | -- Module      :  Numeric.Signal.Internal@@ -15,15 +17,10 @@  module Numeric.Signal.Internal (                 Convolvable(..),-                hamming,-                filter,+                Filterable(..),                 freqz,                 pwelch,                 hilbert,-                complex_power,-                downsample,-                deriv,-                unwrap                 ) where  import Data.Packed.Development(createVector,vec,app1,app2,app3,app4)@@ -44,6 +41,7 @@  type PD = Ptr Double                             type PC = Ptr (Complex Double)                  +type PF = Ptr Float  ----------------------------------------------------------------------------- @@ -53,6 +51,37 @@  ----------------------------------------------------------------------------- +class (Storable a, Container Vector a, Num (Vector a)+      , Convert a, Floating (Vector a), RealElement a+      , Num a)+       => Filterable a where+    -- | conver from Vector Double+    fromDouble :: Vector Double -> Vector a+--       b ~ ComplexOf a, Container Vector b, Convert b) => Filterable a where+    -- | filter a signal+    filter_ :: Vector a -- ^ zero coefficients+            -> Vector a -- ^ pole coefficients+            -> Vector a -- ^ input signal+            -> Vector a -- ^ output signal+    -- | coefficients of a Hamming window+    hamming_ :: Int           -- ^ length+           -> Vector a -- ^ the Hamming coeffficents+    -- | the complex power : real $ v * (conj v)+    complex_power_ :: Vector (Complex Double) -- ^ input+                  -> Vector a                       -- ^ output+     -- | resample, take one sample every n samples in the original+    downsample_ :: Int -> Vector a -> Vector a+    -- | the difference between consecutive elements of a vector+    deriv_ :: Vector a -> Vector a+    -- | unwrap the phase of signal (input expected to be within (-pi,pi)+    unwrap_ :: Vector a -> Vector a+    -- | evaluate a real coefficient polynomial for complex arguments+    polyEval_ :: Vector a           -- ^ the real coefficients+         -> Vector (Complex Double) -- ^ the points at which to be evaluated+         -> Vector (Complex Double) -- ^ the values++-----------------------------------------------------------------------------+ instance Convolvable (Vector Double) where     convolve x y = fst $ fromComplex $ F.ifft $ (F.fft (complex x) * F.fft (complex y)) --    convolve = convolve_vector_double@@ -64,6 +93,17 @@  foreign import ccall "signal-aux.h vector_double_convolve" signal_vector_double_convolve :: CInt -> PD -> CInt -> PD -> CInt -> PD -> IO CInt +instance Convolvable (Vector Float) where+    convolve x y = single $ fst $ fromComplex $ F.ifft $ (F.fft (complex $ double x) * F.fft (complex $ double y))+--    convolve = convolve_vector_double++convolve_vector_float c a = unsafePerformIO $ do+                             r <- createVector (dim a)+                             app3 signal_vector_float_convolve vec c vec a vec r "signalFloatConvolve"+                             return r++foreign import ccall "signal-aux.h vector_float_convolve" signal_vector_float_convolve :: CInt -> PF -> CInt -> PF -> CInt -> PF -> IO CInt+ -----------------------------------------------------------------------------  instance Convolvable (Vector (Complex Double)) where@@ -79,18 +119,52 @@  ----------------------------------------------------------------------------- +instance Filterable Double where+    fromDouble = id+    filter_ = filterD+    hamming_ = hammingD+    complex_power_ = complex_powerD+    downsample_ = downsampleD+    deriv_ = derivD+    unwrap_ = unwrapD+    polyEval_ = polyEval++instance Filterable Float where+    fromDouble = single+    filter_ = filterF+    hamming_ = hammingF+    complex_power_ = complex_powerF+    downsample_ = downsampleF+    deriv_ = derivF+    unwrap_ = unwrapF+    polyEval_ c = polyEval (double c)++-----------------------------------------------------------------------------+ -- | filters the signal-filter :: Vector Double -- ^ zero coefficients+filterD :: Vector Double -- ^ zero coefficients        -> Vector Double -- ^ pole coefficients        -> Vector Double -- ^ input signal        -> Vector Double -- ^ output signal-filter l k v = unsafePerformIO $ do+filterD l k v = unsafePerformIO $ do                r <- createVector (dim v)-               app4 signal_filter vec l vec k vec v vec r "signalFilter"+               app4 signal_filter_double vec l vec k vec v vec r "signalFilter"                return r -foreign import ccall "signal-aux.h filter" signal_filter :: CInt -> PD -> CInt -> PD -> CInt -> PD -> CInt -> PD -> IO CInt+foreign import ccall "signal-aux.h filter_double" signal_filter_double :: CInt -> PD -> CInt -> PD -> CInt -> PD -> CInt -> PD -> IO CInt +-- | filters the signal+filterF :: Vector Float -- ^ zero coefficients+       -> Vector Float -- ^ pole coefficients+       -> Vector Float -- ^ input signal+       -> Vector Float -- ^ output signal+filterF l k v = unsafePerformIO $ do+               r <- createVector (dim v)+               app4 signal_filter_float vec l vec k vec v vec r "signalFilter"+               return r++foreign import ccall "signal-aux.h filter_float" signal_filter_float :: CInt -> PF -> CInt -> PF -> CInt -> PF -> CInt -> PF -> IO CInt+ -----------------------------------------------------------------------------  -- | Hilbert transform with original vector as real value, transformed as imaginary@@ -119,28 +193,51 @@ -----------------------------------------------------------------------------  -- | coefficients of a Hamming window-hamming :: Int           -- ^ length+hammingD :: Int           -- ^ length         -> Vector Double -- ^ the Hamming coeffficents-hamming l +hammingD l      | l == 1          = constant 1.0 1     | otherwise       = unsafePerformIO $ do                         r <- createVector l-                        app1 signal_hamming vec r "Hamming"+                        app1 signal_hamming_double vec r "Hamming"                         return r -foreign import ccall "signal-aux.h hamming" signal_hamming :: CInt -> PD -> IO CInt+foreign import ccall "signal-aux.h hamming_double" signal_hamming_double :: CInt -> PD -> IO CInt +-- | coefficients of a Hamming window+hammingF :: Int           -- ^ length+        -> Vector Float -- ^ the Hamming coeffficents+hammingF l +    | l == 1          = constant 1.0 1+    | otherwise       = unsafePerformIO $ do+                        r <- createVector l+                        app1 signal_hamming_float vec r "Hamming"+                        return r++foreign import ccall "signal-aux.h hamming_float" signal_hamming_float :: CInt -> PF -> IO CInt+ -----------------------------------------------------------------------------  -- | determine the frequency response of a filter-freqz :: Vector Double     -- ^ zero coefficients-      -> Vector Double     -- ^ pole coefficients-      -> Vector Double     -- ^ points (between 0 and 2*pi)-      -> Vector Double     -- ^ response+{-freqz :: (Filterable a, Storable a, Container Vector a, Convert a, RealElement a,+         DoubleOf a ~ DoubleOf (RealOf b), RealElement c, c ~ DoubleOf a, c ~ DoubleOf (RealOf b),+         b ~ Complex a, b ~ ComplexOf a, Convert b, Container Vector b,+         Container Vector c, Convert c, b ~ ComplexOf a, b ~ ComplexOf c)+        ⇒ Vector a     -- ^ zero coefficients+      -> Vector a       -- ^ pole coefficients+      -> Vector a       -- ^ points (between 0 and 2*pi)+      -> Vector a       -- ^ response+-}+freqz :: (Filterable a, Complex Double ~ ComplexOf (DoubleOf a)+        ,Filterable (DoubleOf a)) => +        Vector a       -- ^ zero coefficients+      -> Vector a       -- ^ pole coefficients+      -> Vector a       -- ^ points (between 0 and 2*pi)+      -> Vector a       -- ^ response freqz b a w = let k = max (dim b) (dim a)-                  hb = polyEval (postpad b k) (exp (scale (0 :+ 1) ((complex w) :: Vector (Complex Double))))-                  ha = polyEval (postpad a k) (exp (scale (0 :+ 1) ((complex w) :: Vector (Complex Double))))-              in complex_power (hb / ha)+                  hb = polyEval_ (postpad b k) (exp (scale (0 :+ 1) ((complex $ double w))))+                  ha = polyEval_ (postpad a k) (exp (scale (0 :+ 1) ((complex $ double w))))+              in complex_power_ (hb / ha)  postpad v n = let d = dim v               in if d < n then join [v,(constant 0.0 (n-d))]@@ -162,47 +259,84 @@ -----------------------------------------------------------------------------  -- | the complex power : real $ v * (conj v)-complex_power :: Vector (Complex Double) -- ^ input+complex_powerD :: Vector (Complex Double) -- ^ input               -> Vector Double           -- ^ output-complex_power v = unsafePerformIO $ do+complex_powerD v = unsafePerformIO $ do                   r <- createVector (dim v)-                  app2 signal_complex_power vec v vec r "complex_power"+                  app2 signal_complex_power_double vec v vec r "complex_power"                   return r -foreign import ccall "signal-aux.h complex_power" signal_complex_power :: CInt -> PC -> CInt -> PD -> IO CInt+foreign import ccall "signal-aux.h complex_power_double" signal_complex_power_double :: CInt -> PC -> CInt -> PD -> IO CInt +-- | the complex power : real $ v * (conj v)+complex_powerF :: Vector (Complex Double) -- ^ input+              -> Vector Float             -- ^ output+complex_powerF v = unsafePerformIO $ do+                  r <- createVector (dim v)+                  app2 signal_complex_power_float vec v vec r "complex_power"+                  return r++foreign import ccall "signal-aux.h complex_power_float" signal_complex_power_float :: CInt -> PC -> CInt -> PF -> IO CInt+ -----------------------------------------------------------------------------  -- | resample, take one sample every n samples in the original-downsample :: Int -> Vector Double -> Vector Double-downsample n v = unsafePerformIO $ do+downsampleD :: Int -> Vector Double -> Vector Double+downsampleD n v = unsafePerformIO $ do                r <- createVector (dim v `div` n)-               app2 (signal_downsample $ fromIntegral n) vec v vec r "downsample"+               app2 (signal_downsample_double $ 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+foreign import ccall "signal-aux.h downsample_double" signal_downsample_double :: CInt -> CInt -> PD -> CInt -> PD -> IO CInt +-- | resample, take one sample every n samples in the original+downsampleF :: Int -> Vector Float -> Vector Float+downsampleF n v = unsafePerformIO $ do+               r <- createVector (dim v `div` n)+               app2 (signal_downsample_float $ fromIntegral n) vec v vec r "downsample"+               return r++foreign import ccall "signal-aux.h downsample_float" signal_downsample_float :: CInt -> CInt -> PF -> CInt -> PF -> IO CInt+ -----------------------------------------------------------------------------  -- | the difference between consecutive elements of a vector-deriv :: Vector Double -> Vector Double-deriv v = unsafePerformIO $ do+derivD :: Vector Double -> Vector Double+derivD v = unsafePerformIO $ do           r <- createVector (dim v - 1)-          app2 (signal_deriv) vec v vec r "diff"+          app2 (signal_deriv_double) vec v vec r "diff"           return r -foreign import ccall "signal-aux.h vector_deriv" signal_deriv :: CInt -> PD -> CInt -> PD -> IO CInt+foreign import ccall "signal-aux.h vector_deriv_double" signal_deriv_double :: CInt -> PD -> CInt -> PD -> IO CInt +-- | the difference between consecutive elements of a vector+derivF :: Vector Float -> Vector Float+derivF v = unsafePerformIO $ do+          r <- createVector (dim v - 1)+          app2 (signal_deriv_float) vec v vec r "diff"+          return r++foreign import ccall "signal-aux.h vector_deriv_float" signal_deriv_float :: CInt -> PF -> CInt -> PF -> IO CInt+ -----------------------------------------------------------------------------  -- | unwrap the phase of signal (input expected to be within (-pi,pi)-unwrap :: Vector Double -> Vector Double-unwrap v = unsafePerformIO $ do+unwrapD :: Vector Double -> Vector Double+unwrapD v = unsafePerformIO $ do            r <- createVector $ dim v-           app2 signal_unwrap vec v vec r "unwrap"+           app2 signal_unwrap_double vec v vec r "unwrap"            return r -foreign import ccall "signal-aux.h unwrap" signal_unwrap :: CInt -> PD -> CInt -> PD -> IO CInt+foreign import ccall "signal-aux.h unwrap_double" signal_unwrap_double :: CInt -> PD -> CInt -> PD -> IO CInt++-- | unwrap the phase of signal (input expected to be within (-pi,pi)+unwrapF :: Vector Float -> Vector Float+unwrapF v = unsafePerformIO $ do+           r <- createVector $ dim v+           app2 signal_unwrap_float vec v vec r "unwrap"+           return r++foreign import ccall "signal-aux.h unwrap_float" signal_unwrap_float :: CInt -> PF -> CInt -> PF -> IO CInt  ----------------------------------------------------------------------------- 

@@ -235,7 +235,8 @@   -- | filter the data with the given passband-filter :: (Int,Int) -> Multichannel Double -> Multichannel Double+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 } @@ -253,9 +254,10 @@ -----------------------------------------------------------------------------  -- | calculate the mutual information of the phase between pairs of channels (fills upper half of matrix)-mi_phase :: Multichannel Double      -- ^ input data+mi_phase :: (S.Filterable a, Double ~ DoubleOf a) ⇒+           Multichannel a      -- ^ input data          -> Matrix Double-mi_phase m = let d = _data m+mi_phase m = let d = fmap double $ _data 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)))

@@ -10,6 +10,23 @@  #include <stdio.h> +int vector_float_convolve(int cs, const float* c, int as, const float* a, int rs, float* r)+{+  int h = cs / 2;+  int li, ri;+  int i,j;++  for (i = 0; i < cs; i++) {+    li = i - h;+    ri = i + h;+    r[i] = 0;+    for (j = (li >= 0 ? li : 0); j < (ri < as ? ri : (as - 1)); j++) {+      r[i] += a[j]*c[j+h+1];+    }+  }+  return 0;+}+ int vector_double_convolve(int cs, const double* c, int as, const double* a, int rs, double* r) {   int h = cs / 2;@@ -46,7 +63,7 @@   return 0; } -int filter(int ls, const double* l, int ks, const double* k, int vs, const double* v, int rs, double* r)+int filter_double(int ls, const double* l, int ks, const double* k, int vs, const double* v, int rs, double* r) {   if (ls > vs || ks > vs) return 2000; // BAD_SIZE @@ -70,6 +87,30 @@   return 0; } +int filter_float(int ls, const float* l, int ks, const float* k, int vs, const float* v, int rs, float* r)+{+  if (ls > vs || ks > vs) return 2000; // BAD_SIZE++  int i,j;++  float L = l[0];+  float K = k[0];+  +  int N = ls - 1;+  int M = ks - 1;++  for (i = 0; i < vs; i++) {+    r[i] = 0;+    for (j = 0; j < N; j++) {+      if (i - j > 0) r[i] -= (l[j+1])*v[i-j];+    }+    for (j = 0; j < M; j++) {+      if (i - j > 0) r[i] += (k[j+1])*r[i-j];+    }+  }+  return 0;+}+ int hilbert(int rs, gsl_complex* r) {   int s = rs;@@ -148,7 +189,7 @@   return 0; } -int hamming(int rs, double* r)+int hamming_double(int rs, double* r) {   int i; @@ -157,6 +198,15 @@   return 0; } +int hamming_float(int rs, float* r)+{+  int i;++  for (i = 0; i < rs; i++) r[i] = 0.54 - 0.46 * cos(2*M_PI*i/rs);++  return 0;+}+ int real_poly_complex_eval(int cs, const double* c, int zs, const gsl_complex* z, int rs, gsl_complex* r) {   int i;@@ -167,7 +217,7 @@   return 0; }  -int complex_power(int cs, const gsl_complex* c, int rs, double* r)+int complex_power_double(int cs, const gsl_complex* c, int rs, double* r) {   if (rs != cs) return 2000; // BAD_SIZE @@ -179,7 +229,31 @@   return 0; } -int downsample(int n, int xs, const double* x, int rs, double* r)+int complex_power_float(int cs, const gsl_complex* c, int rs, float* r)+{+  if (rs != cs) return 2000; // BAD_SIZE++  int i;++  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_double(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 downsample_float(int n, int xs, const float* x, int rs, float* r) {   if (rs != xs/n) return 2000; // BAD_SIZE