module Numeric.Signal.Internal (
Convolvable(..),
Filterable(..),
freqz,
pwelch,
hilbert
) where
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import qualified Numeric.GSL.Fourier as F
import Foreign
import Foreign.C.Types
import Prelude hiding(filter)
import System.IO.Unsafe(unsafePerformIO)
infixr 1 #
a # b = applyRaw a b
type PD = Ptr Double
type PC = Ptr (Complex Double)
type PF = Ptr Float
class Convolvable a where
convolve :: a -> a -> a
class (Storable a, Container Vector a, Num (Vector a)
, Convert a, Floating (Vector a), RealElement a
, Num a)
=> Filterable a where
fromDouble :: Vector Double -> Vector a
filter_ :: Vector a
-> Vector a
-> Vector a
-> Vector a
hamming_ :: Int
-> Vector a
complex_power_ :: Vector (Complex Double)
-> Vector a
downsample_ :: Int -> Vector a -> Vector a
deriv_ :: Vector a -> Vector a
unwrap_ :: Vector a -> Vector a
polyEval_ :: Vector a
-> Vector (Complex Double)
-> Vector (Complex Double)
cross_covariance_ :: Int
-> Vector a
-> Vector a
-> (a,a,Vector a)
cumulative_sum_ :: Vector a
-> Vector a
instance Convolvable (Vector Double) where
convolve x y = fst $ fromComplex $ F.ifft $ (F.fft (complex x) * F.fft (complex y))
convolve_vector_double c a = unsafePerformIO $ do
r <- createVector (size a)
(c # a # r # id) signal_vector_double_convolve #| "signalDoubleConvolve"
return r
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_vector_float c a = unsafePerformIO $ do
r <- createVector (size a)
(c # a # r # id ) signal_vector_float_convolve #| "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
convolve x y = F.ifft $ (F.fft x * F.fft y)
convolve_vector_complex c a = unsafePerformIO $ do
r <- createVector (size a)
(c # a # r # id) signal_vector_complex_convolve #| "signalComplexConvolve"
return r
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
fromDouble = id
filter_ = filterD
hamming_ = hammingD
complex_power_ = complex_powerD
downsample_ = downsampleD
deriv_ = derivD
unwrap_ = unwrapD
polyEval_ = polyEval
cross_covariance_ = crossCovarianceD
cumulative_sum_ = cumSumD
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)
cross_covariance_ = crossCovarianceF
cumulative_sum_ = cumSumF
filterD :: Vector Double
-> Vector Double
-> Vector Double
-> Vector Double
filterD l k v = unsafePerformIO $ do
r <- createVector (size v)
(l # k # v # r # id) signal_filter_double #| "signalFilter"
return r
foreign import ccall "signal-aux.h filter_double" signal_filter_double :: CInt -> PD -> CInt -> PD -> CInt -> PD -> CInt -> PD -> IO CInt
filterF :: Vector Float
-> Vector Float
-> Vector Float
-> Vector Float
filterF l k v = unsafePerformIO $ do
r <- createVector (size v)
(l # k # v # r # id) signal_filter_float #| "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 :: Vector Double -> Vector (Complex Double)
hilbert v = unsafePerformIO $ do
let r = complex v
(r # id) signal_hilbert #| "hilbert"
return r
foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> IO CInt
pwelch :: Int
-> Vector Double
-> Vector Double
pwelch w v = unsafePerformIO $ do
let r = konst 0.0 ((w `div` 2) + 1)
(complex v # r # id) (signal_pwelch $ fromIntegral w) #| "pwelch"
return r
foreign import ccall "signal-aux.h pwelch" signal_pwelch :: CInt -> CInt -> PC -> CInt -> PD -> IO CInt
hammingD :: Int
-> Vector Double
hammingD l
| l == 1 = konst 1.0 1
| otherwise = unsafePerformIO $ do
r <- createVector l
(r # id) signal_hamming_double #| "Hamming"
return r
foreign import ccall "signal-aux.h hamming_double" signal_hamming_double :: CInt -> PD -> IO CInt
hammingF :: Int
-> Vector Float
hammingF l
| l == 1 = konst 1.0 1
| otherwise = unsafePerformIO $ do
r <- createVector l
(r # id) signal_hamming_float #| "Hamming"
return r
foreign import ccall "signal-aux.h hamming_float" signal_hamming_float :: CInt -> PF -> IO CInt
freqz :: (Filterable a, Complex Double ~ ComplexOf (DoubleOf a)
,Filterable (DoubleOf a)) =>
Vector a
-> Vector a
-> Vector a
-> Vector a
freqz b a w = let k = max (size b) (size a)
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 = size v
in if d < n then vjoin [v,(konst 0.0 (nd))]
else v
polyEval :: Vector Double
-> Vector (Complex Double)
-> Vector (Complex Double)
polyEval c z = unsafePerformIO $ do
r <- createVector (size z)
(c # z # r # id) signal_real_poly_complex_eval #| "polyEval"
return r
foreign import ccall "signal-aux.h real_poly_complex_eval" signal_real_poly_complex_eval :: CInt -> PD -> CInt -> PC -> CInt -> PC -> IO CInt
complex_powerD :: Vector (Complex Double)
-> Vector Double
complex_powerD v = unsafePerformIO $ do
r <- createVector (size v)
(v # r # id) signal_complex_power_double #| "complex_power"
return r
foreign import ccall "signal-aux.h complex_power_double" signal_complex_power_double :: CInt -> PC -> CInt -> PD -> IO CInt
complex_powerF :: Vector (Complex Double)
-> Vector Float
complex_powerF v = unsafePerformIO $ do
r <- createVector (size v)
(v # r # id) signal_complex_power_float #| "complex_power"
return r
foreign import ccall "signal-aux.h complex_power_float" signal_complex_power_float :: CInt -> PC -> CInt -> PF -> IO CInt
downsampleD :: Int -> Vector Double -> Vector Double
downsampleD n v = unsafePerformIO $ do
r <- createVector (size v `div` n)
(v # r # id) (signal_downsample_double $ fromIntegral n) #| "downsample"
return r
foreign import ccall "signal-aux.h downsample_double" signal_downsample_double :: CInt -> CInt -> PD -> CInt -> PD -> IO CInt
downsampleF :: Int -> Vector Float -> Vector Float
downsampleF n v = unsafePerformIO $ do
r <- createVector (size v `div` n)
(v # r # id) (signal_downsample_float $ fromIntegral n) #| "downsample"
return r
foreign import ccall "signal-aux.h downsample_float" signal_downsample_float :: CInt -> CInt -> PF -> CInt -> PF -> IO CInt
derivD :: Vector Double -> Vector Double
derivD v = unsafePerformIO $ do
r <- createVector (size v 1)
(v # r # id) (signal_diff_double) #| "diff"
return r
foreign import ccall "signal-aux.h vector_diff_double" signal_diff_double :: CInt -> PD -> CInt -> PD -> IO CInt
derivF :: Vector Float -> Vector Float
derivF v = unsafePerformIO $ do
r <- createVector (size v 1)
(v # r # id) (signal_diff_float) #| "diff"
return r
foreign import ccall "signal-aux.h vector_diff_float" signal_diff_float :: CInt -> PF -> CInt -> PF -> IO CInt
unwrapD :: Vector Double -> Vector Double
unwrapD v = unsafePerformIO $ do
r <- createVector $ size v
(v # r # id) signal_unwrap_double #| "unwrap"
return r
foreign import ccall "signal-aux.h unwrap_double" signal_unwrap_double :: CInt -> PD -> CInt -> PD -> IO CInt
unwrapF :: Vector Float -> Vector Float
unwrapF v = unsafePerformIO $ do
r <- createVector $ size v
(v # r # id) signal_unwrap_float #| "unwrap"
return r
foreign import ccall "signal-aux.h unwrap_float" signal_unwrap_float :: CInt -> PF -> CInt -> PF -> IO CInt
crossCovarianceD :: Int -> Vector Double -> Vector Double -> (Double,Double,Vector Double)
crossCovarianceD l x y = unsafePerformIO $ do
r <- createVector (2*l)
alloca $ \sx ->
alloca $ \sy -> do
(x # y # r # id) (signal_cross_covariance_double (fromIntegral l) sx sy) #| "cross_covariance"
sx' <- peek sx
sy' <- peek sy
return (sx',sy',r)
foreign import ccall "signal-aux.h cross_covariance_double"
signal_cross_covariance_double :: CInt -> PD -> PD -> CInt -> PD -> CInt
-> PD -> CInt -> PD -> IO CInt
crossCovarianceF :: Int -> Vector Float -> Vector Float -> (Float,Float,Vector Float)
crossCovarianceF l x y = unsafePerformIO $ do
r <- createVector (2*l)
alloca $ \sx ->
alloca $ \sy -> do
(x # y # r # id) (signal_cross_covariance_float (fromIntegral l) sx sy) #| "cross_covariance"
sx' <- peek sx
sy' <- peek sy
return (sx',sy',r)
foreign import ccall "signal-aux.h cross_covariance_float"
signal_cross_covariance_float :: CInt -> PF -> PF -> CInt -> PF -> CInt
-> PF -> CInt -> PF -> IO CInt
cumSumD :: Vector Double -> Vector Double
cumSumD v = unsafePerformIO $ do
r <- createVector (size v)
(v # r # id) signal_cum_sum_double #| "cumSumD"
return r
cumSumF :: Vector Float -> Vector Float
cumSumF v = unsafePerformIO $ do
r <- createVector (size v)
(v # r # id) signal_cum_sum_float #| "cumSumF"
return r
foreign import ccall "signal-aux.h cum_sum_double"
signal_cum_sum_double :: CInt -> PD -> CInt -> PD -> IO CInt
foreign import ccall "signal-aux.h cum_sum_float"
signal_cum_sum_float :: CInt -> PF -> CInt -> PF -> IO CInt