module Numeric.Signal (
hamming,
pwelch,
fir,standard_fir,broadband_fir,
freqzF,freqzN,
filter,broadband_filter,
analytic_signal,analytic_power,analytic_phase
) where
import qualified Numeric.Signal.Internal as S
import Complex
import qualified Data.List as L
import Data.Packed.Vector
import Data.Packed(Container(..))
import Numeric.GSL.Vector
import Numeric.LinearAlgebra.Instances()
import Numeric.LinearAlgebra.Linear(Linear(..))
import qualified Numeric.GSL.Fourier as F
import Prelude hiding(filter)
filter :: Vector Double
-> Vector Double
-> Int
-> Vector Double
-> Vector Double
filter b a s v = let sd = (fromIntegral s) * 2.0
sc = recip $ (sd / (fromIntegral $ max (dim b) (dim a))) * sd
in scale sc $ S.filter b a v
hamming :: Int
-> Vector Double
hamming = S.hamming
pwelch :: Int
-> Int
-> Vector Double
-> (Vector Double,Vector Double)
pwelch s w v = let w' = max s w
r = S.pwelch w' v
sd = recip $ (fromIntegral s)/2
r' = scale sd r
f = linspace ((w `div` 2) + 1) (0,sd)
in (f,r')
broadband_fir :: Int
-> (Int,Int)
-> Vector Double
broadband_fir s (l,h) = let o = 501
ny = (fromIntegral s) / 2.0
fl = (fromIntegral l) / ny
fh = (fromIntegral h) / ny
f = [0, fl*0.95, fl, fh, fh*1.05, 1]
m = [0,0,1,1,0,0]
be = zip f m
in standard_fir o be
broadband_filter :: Int
-> (Int,Int)
-> Vector Double
-> Vector Double
broadband_filter s f v = let b = broadband_fir s f
in S.filter b (constant 1.0 1) v
standard_fir :: Int -> [(Double,Double)] -> Vector Double
standard_fir o be = let grid = calc_grid o
trans = grid `div` 16
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
in floor $ 2.0 ** ((fromIntegral next_power) :: Double)
fir :: Int
-> [(Double,Double)]
-> Int
-> Int
-> Vector Double
-> Vector Double
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
m' = interpolate f m f'
grid = interpolate f' m' $ map (\x -> (fromIntegral x)/(fromIntegral gn)) [0..(gn1)]
grid' = map (\x -> x :+ 0) grid
(b,_) = fromComplex $ F.ifft $ fromList $ grid' ++ (reverse (drop 1 grid'))
b' = join [subVector ((dim b)mid1) (mid+1) b, subVector 1 (mid+1) b]
in b' * w
floor_zero x
| x < 0.0 = 0.0
| otherwise = x
ceil_one x
| x > 1.0 = 1.0
| otherwise = x
diff :: Double -> [Double] -> [Double]
diff _ [] = []
diff _ [x] = [x]
diff inc (x1:x2:xs)
| x1 == x2 = (floor_zero $ x1inc):x1:(ceil_one $ x1+inc):(diff inc (L.filter (/= x2) xs))
| otherwise = x1:(diff inc (x2:xs))
interpolate :: [Double] -> [Double] -> [Double] -> [Double]
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' x y xp = let Just i = L.findIndex (> xp) x
in (interpolate'' i x y xp)
interpolate'' :: Int -> [Double] -> [Double] -> Double -> Double
interpolate'' i x y xp = let x0 = x !! (i1)
y0 = y !! (i1)
x1 = x !! i
y1 = y !! i
in y0 + (xp x0) * ((y1 y0)/(x1x0))
freqzF :: Vector Double
-> Vector Double
-> Int
-> Vector Double
-> Vector Double
freqzF b a s f = S.freqz b a ((2*pi/(fromIntegral s)) * f)
freqzN :: Vector Double
-> Vector Double
-> Int
-> Int
-> (Vector Double,Vector Double)
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)
analytic_signal :: Vector Double -> Vector (Complex Double)
analytic_signal = S.hilbert
analytic_power :: Vector (Complex Double) -> Vector Double
analytic_power = S.complex_power
analytic_phase :: Vector (Complex Double) -> Vector Double
analytic_phase v = let (r,c) = fromComplex v
in vectorZipR ATan2 c r