hsignal (empty) → 0.1.0.1
raw patch · 11 files changed
+790/−0 lines, 11 filesdep +basedep +haskell98dep +hmatrixbuild-type:Customsetup-changed
Dependencies added: base, haskell98, hmatrix, storable-complex
Files
- CHANGES +2/−0
- INSTALL +35/−0
- LICENSE +2/−0
- README +12/−0
- Setup.lhs +5/−0
- configure +3/−0
- configure.hs +132/−0
- hsignal.cabal +55/−0
- lib/Numeric/Signal.hs +200/−0
- lib/Numeric/Signal/Internal.hs +169/−0
- lib/Numeric/Signal/signal-aux.c +175/−0
+ CHANGES view
@@ -0,0 +1,2 @@+0.1.0.1: + first version
+ INSTALL view
@@ -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+
+ LICENSE view
@@ -0,0 +1,2 @@+Copyright Alberto Ruiz 2006-2007+GPL license
+ README view
@@ -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
+ Setup.lhs view
@@ -0,0 +1,5 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple++> main = defaultMainWithHooks autoconfUserHooks
+ configure view
@@ -0,0 +1,3 @@+#! /bin/sh++runhaskell configure.hs $*
+ configure.hs view
@@ -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 <gsl/gsl_sf_gamma.h>\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 <gsl/gsl_sf_gamma.h>\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 <gsl/gsl_sf_exp.h>\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"
+ hsignal.cabal view
@@ -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 <at> gmail <dot> 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+
+ lib/Numeric/Signal.hs view
@@ -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 <at> gmail <dot> 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++-----------------------------------------------------------------------------
+ lib/Numeric/Signal/Internal.hs view
@@ -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 <at> gmail <dot> 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++-----------------------------------------------------------------------------
+ lib/Numeric/Signal/signal-aux.c view
@@ -0,0 +1,175 @@+#include <gsl/gsl_complex.h>++#include <gsl/gsl_math.h>+#include <gsl/gsl_fft_real.h>+#include <gsl/gsl_fft_complex.h>+#include <gsl/gsl_vector.h>+#include <gsl/gsl_blas.h>+#include <gsl/gsl_poly.h>++#include <stdio.h>++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;+}