diff --git a/CHANGES b/CHANGES new file mode 100644 --- /dev/null +++ b/CHANGES @@ -0,0 +1,2 @@ +0.1.0.1: + first version diff --git a/INSTALL b/INSTALL new file mode 100644 --- /dev/null +++ b/INSTALL @@ -0,0 +1,35 @@ +----------------------------------------------- + A simple 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 + +6) Open a terminal, cd to the hsignal folder, and run + + > cabal install + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +Copyright Alberto Ruiz 2006-2007 +GPL license diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,12 @@ +----------------------------------------- + A simple scientific library for Haskell +----------------------------------------- + +INSTALLATION + +See the INSTALL file. + +ACKNOWLEDGEMENTS ----------------------------------------------------- + +I thank ALberto Ruiz for hmatrix, especially for setup files which I +have shamelessly copied and modified diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,5 @@ +#! /usr/bin/env runhaskell + +> import Distribution.Simple + +> main = defaultMainWithHooks autoconfUserHooks diff --git a/configure b/configure new file mode 100644 --- /dev/null +++ b/configure @@ -0,0 +1,3 @@ +#! /bin/sh + +runhaskell configure.hs $* diff --git a/configure.hs b/configure.hs new file mode 100644 --- /dev/null +++ b/configure.hs @@ -0,0 +1,132 @@ +#! /usr/bin/env runhaskell +{- configure.hs for hsignal, copied from hmatrix + ------------------------ + + GSL and LAPACK may require auxiliary libraries which depend on OS, + distribution, and implementation. This script tries to to find out + the correct link command for your system. + Suggestions and contributions are welcome. + + By default we try to link -lgsl -llapack. This works in ubuntu/debian, + both with and without ATLAS. + If this fails we try different sets of additional libraries which are + known to work in some systems. + + The desired libraries can also be explicitly given by the user using cabal + flags (e.g., -fmkl, -faccelerate) or --configure-option=link:lib1,lib2,lib3,... + +-} + +import System +import Data.List(isPrefixOf, intercalate) +import Distribution.Simple.LocalBuildInfo +import Distribution.Simple.Configure +import Distribution.PackageDescription + +-- possible additional dependencies for the desired libs (by default gsl lapack) + +opts = [ "" -- Ubuntu/Debian + , "blas" + , "blas cblas" + , "cblas" + , "gslcblas" + , "blas gslcblas" + , "f77blas" + , "f77blas cblas atlas gcc_s" -- Arch Linux (older version of atlas-lapack) + , "blas gslcblas gfortran" -- Arch Linux with normal blas and lapack + ] + +-- compile a simple program with symbols from GSL and LAPACK with the given libs +testprog buildInfo libs fmks = + "echo \"#include \nint main(){zgesvd_(); gsl_sf_gamma(5);}\"" + ++" > /tmp/dummy.c; gcc " + ++ (join $ ccOptions buildInfo) ++ " " + ++ (join $ cppOptions buildInfo) ++ " " + ++ (join $ map ("-I"++) $ includeDirs buildInfo) + ++" /tmp/dummy.c -o /tmp/dummy " + ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " " + ++ (prepend "-l" $ libs) ++ " " + ++ (prepend "-framework " fmks) ++ " > /dev/null 2> /dev/null" + +join = intercalate " " +prepend x = unwords . map (x++) . words + +check buildInfo libs fmks = (ExitSuccess ==) `fmap` system (testprog buildInfo libs fmks) + +-- simple test for GSL +gsl buildInfo = "echo \"#include \nint main(){gsl_sf_gamma(5);}\"" + ++" > /tmp/dummy.c; gcc " + ++ (join $ ccOptions buildInfo) ++ " " + ++ (join $ cppOptions buildInfo) ++ " " + ++ (join $ map ("-I"++) $ includeDirs buildInfo) + ++ " /tmp/dummy.c -o /tmp/dummy " + ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" + ++ " > /dev/null 2> /dev/null" + +-- test for gsl >= 1.12 +gsl112 buildInfo = + "echo \"#include \nint main(){gsl_sf_exprel_n_CF_e(1,1,0);}\"" + ++" > /tmp/dummy.c; gcc /tmp/dummy.c " + ++ (join $ ccOptions buildInfo) ++ " " + ++ (join $ cppOptions buildInfo) ++ " " + ++ (join $ map ("-I"++) $ includeDirs buildInfo) + ++" -o /tmp/dummy " + ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" + ++ " > /dev/null 2> /dev/null" + + +checkCommand c = (ExitSuccess ==) `fmap` system c + +-- test different configurations until the first one works +try _ _ _ [] = return Nothing +try i b f (opt:rest) = do + ok <- check i (b ++ " " ++ opt) f + if ok then return (Just opt) + else try i b f rest + +-- read --configure-option=link:lib1,lib2,lib3,etc +linkop = "link:" +getUserLink = concatMap (g . drop (length linkop)) . filter (isPrefixOf linkop) + where g = map cs + cs ',' = ' ' + cs x = x + +main = do + putStr "Checking foreign libraries..." + + args <- getArgs + Just bInfo <- maybeGetPersistBuildConfig "dist" + + let Just lib = library . localPkgDescr $ bInfo + buildInfo = libBuildInfo lib + base = unwords . extraLibs $ buildInfo + fwks = unwords . frameworks $ buildInfo + auxpref = getUserLink args + + -- We extract the desired libs from hsignal.cabal (using a cabal flags) + -- and from a posible --configure-option=link:lib1,lib2,lib3 + -- by default the desired libs are gsl lapack. + + let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref + fullOpts = map ((pref++" ")++) opts + + r <- try buildInfo base fwks fullOpts + case r of + Nothing -> do + putStrLn " FAIL" + g <- checkCommand $ gsl buildInfo + if g + then putStrLn " *** Sorry, I can't link LAPACK." + else putStrLn " *** Sorry, I can't link GSL." + putStrLn " *** Please make sure that the appropriate -dev packages are installed." + putStrLn " *** You can also specify the required libraries using" + putStrLn " *** cabal install hsignal --configure-option=link:lib1,lib2,lib3,etc." + writeFile "hsignal.buildinfo" ("buildable: False\n") + Just ops -> do + putStrLn " OK" + g <- checkCommand $ gsl112 buildInfo + writeFile "hsignal.buildinfo" $ "extra-libraries: " ++ + ops ++ "\n" ++ + if g + then "" + else "cc-options: -DGSL110\n" diff --git a/hsignal.cabal b/hsignal.cabal new file mode 100644 --- /dev/null +++ b/hsignal.cabal @@ -0,0 +1,55 @@ +Name: hsignal +Version: 0.1.0.1 +License: GPL +License-file: LICENSE +Author: Alexander Vivian Hugh McPhail +Maintainer: haskell.vivian.mcphail gmail com +Stability: provisional +Homepage: http://code.haskell.org/hsignal +Synopsis: Signal processing +Description: Purely functional interface to signal processing based on hmatrix +Category: Math +tested-with: GHC ==6.10.4 + +cabal-version: >=1.2 + +build-type: Custom + +extra-source-files: configure configure.hs README INSTALL CHANGES +extra-tmp-files: hmatrix.buildinfo + +library + + Build-Depends: base >= 3 && < 5, + storable-complex, haskell98, + hmatrix + + Extensions: ForeignFunctionInterface + + hs-source-dirs: lib + Exposed-modules: Numeric.Signal + other-modules: Numeric.Signal.Internal + C-sources: lib/Numeric/Signal/signal-aux.c + + ghc-prof-options: -auto + + ghc-options: -Wall -fno-warn-missing-signatures + -fno-warn-orphans + -fno-warn-unused-binds + + if os(OSX) + extra-lib-dirs: /opt/local/lib/ + include-dirs: /opt/local/include/ + extra-libraries: gsl + frameworks: Accelerate + +-- The extra-libraries required for GSL +-- should now be automatically detected by configure(.hs) + + extra-libraries: + extra-lib-dirs: + + source-repository head + type: darcs + location: http://code.haskell.org/hsignal + diff --git a/lib/Numeric/Signal.hs b/lib/Numeric/Signal.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Signal.hs @@ -0,0 +1,200 @@ +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Signal +-- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : haskell.vivian.mcphail gmail com +-- Stability : provisional +-- Portability : uses FFI +-- +-- Signal processing functions +-- +----------------------------------------------------------------------------- + +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) + +----------------------------------------------------------------------------- + +-- | filters the signal +filter :: Vector Double -- ^ zero coefficients + -> Vector Double -- ^ pole coefficients + -> Int -- ^ sampling rate + -> Vector Double -- ^ input signal + -> Vector Double -- ^ output signal +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 + +----------------------------------------------------------------------------- + +-- | 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 + -> Vector Double -- ^ input signal + -> (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 + -- scale for sampling rate + r' = scale sd r + f = linspace ((w `div` 2) + 1) (0,sd) + in (f,r') + +----------------------------------------------------------------------------- + +-- | a broadband FIR +broadband_fir :: Int -- ^ sampling rate + -> (Int,Int) -- ^ (lower,upper) frequency cutoff + -> Vector Double -- ^ filter coefficients +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 + +-- | 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 + in S.filter b (constant 1.0 1) 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 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) + + +-- | 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 + -> Int -- ^ grid spacing + -> Int -- ^ transition width + -> Vector Double -- ^ smoothing window (size is order + 1) + -> Vector Double -- ^ 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 + 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' = join [subVector ((dim b)-mid-1) (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 $ 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 _ _ [] = [] +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 !! (i-1) + y0 = y !! (i-1) + x1 = x !! i + y1 = y !! i + in y0 + (xp - x0) * ((y1 - y0)/(x1-x0)) + +----------------------------------------------------------------------------- + +-- | 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 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 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) + +----------------------------------------------------------------------------- + +-- | an analytic signal is the original signal with Hilbert-transformed signal as imaginary component +analytic_signal :: Vector Double -> Vector (Complex Double) +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 + +-- | the phase of an analytic signal +analytic_phase :: Vector (Complex Double) -> Vector Double +analytic_phase v = let (r,c) = fromComplex v + in vectorZipR ATan2 c r + +----------------------------------------------------------------------------- diff --git a/lib/Numeric/Signal/Internal.hs b/lib/Numeric/Signal/Internal.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Signal/Internal.hs @@ -0,0 +1,169 @@ +{-# OPTIONS_GHC -fglasgow-exts #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Signal.Internal +-- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : haskell.vivian.mcphail gmail com +-- Stability : provisional +-- Portability : uses FFI +-- +-- low-level interface +-- +----------------------------------------------------------------------------- + +module Numeric.Signal.Internal ( + Convolvable(..), + hamming, + filter, + freqz, + hilbert, + pwelch, + complex_power + ) where + +import Data.Packed.Development(createVector,vec,app1,app2,app3,app4) +import Data.Packed.Vector +import Data.Packed(Container(..)) + +import Numeric.LinearAlgebra.Instances() +import Numeric.LinearAlgebra.Linear(Linear(..)) + +import Foreign +import Complex +import Foreign.C.Types + +import Prelude hiding(filter) + +----------------------------------------------------------------------------- + +type PD = Ptr Double +type PC = Ptr (Complex Double) + +----------------------------------------------------------------------------- + +class Convolvable a where + -- | convolve two containers, output is the size of the second argument, no zero padding + convolve :: a -> a -> a + +----------------------------------------------------------------------------- + +instance Convolvable (Vector Double) where + convolve = convolve_vector_double + +convolve_vector_double c a = unsafePerformIO $ do + r <- createVector (dim a) + app3 signal_vector_double_convolve vec c vec a vec r "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 (Complex Double)) where + convolve = convolve_vector_complex + +convolve_vector_complex c a = unsafePerformIO $ do + r <- createVector (dim a) + app3 signal_vector_complex_convolve vec c vec a vec r "signalComplexConvolve" + return r + +foreign import ccall "signal-aux.h vector_complex_convolve" signal_vector_complex_convolve :: CInt -> PC -> CInt -> PC -> CInt -> PC -> IO CInt + +----------------------------------------------------------------------------- + +-- | filters the signal +filter :: Vector Double -- ^ zero coefficients + -> Vector Double -- ^ pole coefficients + -> Vector Double -- ^ input signal + -> Vector Double -- ^ output signal +filter l k v = unsafePerformIO $ do + r <- createVector (dim v) + app4 signal_filter 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 + +----------------------------------------------------------------------------- + +-- | Hilbert transform with original vector as real value, transformed as imaginary +hilbert :: Vector Double -> Vector (Complex Double) +hilbert v = unsafePerformIO $ do + r <- createVector (dim v) + let v' = complex v + -- could use (comp v) to make a complex vector in haskell rather than C + app2 signal_hilbert vec v' vec r "hilbert" + return r + +foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> CInt -> PC -> IO CInt + +----------------------------------------------------------------------------- + +-- | Welch (1967) power spectrum density using periodogram/FFT method +pwelch :: Int -- ^ window size (multiple of 2) + -> Vector Double -- ^ input signal + -> Vector Double -- ^ power density +pwelch w v = unsafePerformIO $ do + let r = constant 0.0 ((w `div` 2) + 1) + app2 (signal_pwelch $ fromIntegral w) vec (complex v) vec r "pwelch" + return r + +foreign import ccall "signal-aux.h pwelch" signal_pwelch :: CInt -> CInt -> PC -> CInt -> PD -> IO CInt + +----------------------------------------------------------------------------- + +-- | coefficients of a Hamming window +hamming :: Int -- ^ length + -> Vector Double -- ^ the Hamming coeffficents +hamming l + | l == 1 = constant 1.0 1 + | otherwise = unsafePerformIO $ do + r <- createVector l + app1 signal_hamming vec r "Hamming" + return r + +foreign import ccall "signal-aux.h hamming" signal_hamming :: CInt -> PD -> 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 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) + +postpad v n = let d = dim v + in if d < n then join [v,(constant 0.0 (n-d))] + else v + +----------------------------------------------------------------------------- + +-- | evaluate a real coefficient polynomial for complex arguments +polyEval :: Vector Double -- ^ the real coefficients + -> Vector (Complex Double) -- ^ the points at which to be evaluated + -> Vector (Complex Double) -- ^ the values +polyEval c z = unsafePerformIO $ do + r <- createVector (dim z) + app3 signal_real_poly_complex_eval vec c vec z vec r "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 + +----------------------------------------------------------------------------- + +-- | the complex power : real $ v * (conj c) +complex_power :: Vector (Complex Double) -- ^ input + -> Vector Double -- ^ output +complex_power v = unsafePerformIO $ do + r <- createVector (dim v) + app2 signal_complex_power 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 + +----------------------------------------------------------------------------- diff --git a/lib/Numeric/Signal/signal-aux.c b/lib/Numeric/Signal/signal-aux.c new file mode 100644 --- /dev/null +++ b/lib/Numeric/Signal/signal-aux.c @@ -0,0 +1,175 @@ +#include + +#include +#include +#include +#include +#include +#include + +#include + +int vector_double_convolve(int cs, const double* c, int as, const double* a, int rs, double* 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_complex_convolve(int cs, const gsl_complex* c, int as, const gsl_complex* a, int rs, gsl_complex* 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].dat[0] = 0; + r[i].dat[1] = 0; + for (j = (li >= 0 ? li : 0); j < (ri < as ? ri : (as - 1)); j++) { + r[i].dat[0] += a[j].dat[0]*c[j+h+1].dat[0]-a[j].dat[1]*c[j+h+1].dat[1]; + r[i].dat[1] += a[j].dat[0]*c[j+h+1].dat[1]+a[j].dat[1]*c[j+h+1].dat[0]; + } + } + return 0; +} + +int filter(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 + + int i,j; + + double L = l[0]; + double K = k[0]; + + int N = ls - 1; + int M = ks - 1; + + for (i = 0; i < vs; i++) { + r[i] = 0; + for (j = 1; j < N; j++) { + if (i - j > 0) r[i] -= (l[j+1]/L)*v[i-j]; + } + for (j = 0; j < M; j++) { + if (i - j > 0) r[i] += (k[j+1]/K)*r[i-j]; + } + } + return 0; +} + +int hilbert(int vs, const gsl_complex* v, int rs, gsl_complex* r) +{ + if (vs != rs) return 2000; // BAD_SIZE + + int s = vs; + + gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s); + gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s); + + // forward fourier transform + gsl_fft_complex_forward ((double*)r, 1, s, wavetable, workspace); + // zero negative coefficients and double positive + int i; + for (i = 0; i < s; i++) { + if (i <= s/2) { + r[i].dat[0] *= sqrt(2.0); + r[i].dat[1] *= sqrt(2.0); + } + else { + r[i].dat[0] = 0.0; + r[i].dat[1] = 0.0; + } + } + // inverse fourier transform + gsl_fft_complex_inverse ((double*)r, 1, s, wavetable, workspace); + + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); + + return 0; +} + +int pwelch(int w, int vs, const gsl_complex* v, int rs, double* r) +{ + if (w > vs) return 2000; // BAD_SIZE + + int i,j; + + int fs = w; + + int num_windows = vs / fs; // ignore end + + double s[fs]; + for (i = 0; i < fs; i++) s[i] = 0; + + gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (fs); + gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (fs); + + gsl_complex* f = malloc(sizeof(gsl_complex)*fs); + gsl_vector_view F = gsl_vector_view_array((double*)f, 2*fs); + + gsl_vector_view X; + + for (i = 0; i < num_windows; i++) { + X = gsl_vector_view_array((double*)(&v[i*fs]), 2*fs); // v is gsl_complex* + gsl_blas_dcopy(&X.vector,&F.vector); + gsl_fft_complex_forward ((double*)f, 1, fs, wavetable, workspace); + for (j = 0; j < fs; j++) s[j] += f[j].dat[0]*f[j].dat[0] + f[j].dat[1]*f[j].dat[1]; + } + for (j = 0; j < rs; j++) { + if (j == 0) r[j] = s[j]; + else if (j == (rs-1)) r[j] = s[j]; + else r[j] = s[j] + s[fs-j+1]; + + r[j] /= (num_windows * w / 2); + } + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); + + free(f); + + return 0; +} + +int hamming(int rs, double* 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; + + for (i = 0; i < zs; i++) + r[i] = gsl_poly_complex_eval(c,cs,z[i]); + + return 0; +} + +int complex_power(int cs, const gsl_complex* c, int rs, double* 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; +}