diff --git a/DSP/Basic.hs b/DSP/Basic.hs --- a/DSP/Basic.hs +++ b/DSP/Basic.hs @@ -20,6 +20,28 @@ -- * Functions +-- | 'linspace' generates a list of values linearly spaced between specified +-- start and end values (array will include both start and end values). +-- +-- @linspace 0.0 1.0 5 == [ 0.0, 0.25, 0.5, 0.75 1.0 ]@ + +linspace :: Double -> Double -> Int -> [Double] +linspace a b n = + map (\ i -> a + (fromIntegral i) * inc) [(0::Int) .. (n - 1)] + where + inc = (b - a) / (fromIntegral (n - 1)) + + +-- | 'logspace' generates a list of values logarithmically spaced between the +-- values 10 ** start and 10 ** end (array will include both start and end values). +-- +-- @logspace 0.0 1.0 4 == [ 1.0, 2.1544, 4.6416, 10.0 ]@ + +logspace :: Double -> Double -> Int -> [Double] +logspace a b n = + map (\ x -> 10.0 ** x) $ linspace a b n + + -- | 'delay' is the unit delay function, eg, -- -- @delay1 [ 1, 2, 3 ] == [ 0, 1, 2, 3 ]@ diff --git a/DSP/Filter/FIR/FIR.hs b/DSP/Filter/FIR/FIR.hs --- a/DSP/Filter/FIR/FIR.hs +++ b/DSP/Filter/FIR/FIR.hs @@ -28,7 +28,7 @@ {-# specialize fir :: Array Int Float -> [Float] -> [Float] #-} {-# specialize fir :: Array Int Double -> [Double] -> [Double] #-} -fir :: Num a => Array Int a -- ^ h[n] +fir :: (Num a, Eq a) => Array Int a -- ^ h[n] -> [a] -- ^ x[n] -> [a] -- ^ y[n] @@ -41,7 +41,7 @@ where w = listArray (0,m) $ x : replicate m 0 m = snd $ bounds h --- This is for testing the symetric helpers. +-- This is for testing the symmetric helpers. fir0 :: Num a => Array Int a -> [a] -> [a] fir0 _ [] = [] @@ -49,7 +49,7 @@ where w = listArray (0,m) $ x : replicate m 0 m = snd $ bounds h --- Asymetric FIR +-- Asymmetric FIR {-# specialize fir'0 :: Array Int Float -> Array Int Float -> [Float] -> [Float] #-} {-# specialize fir'0 :: Array Int Double -> Array Int Double -> [Double] -> [Double] #-} @@ -63,7 +63,7 @@ w' = listArray (0,m) $ x : elems w m = snd $ bounds h --- Type 1: symetric FIR, even order / odd length +-- Type 1: symmetric FIR, even order / odd length {-# specialize fir'1 :: Array Int Float -> Array Int Float -> [Float] -> [Float] #-} {-# specialize fir'1 :: Array Int Double -> Array Int Double -> [Double] -> [Double] #-} @@ -79,7 +79,7 @@ m = snd $ bounds h m2 = m `div` 2 --- Type 2: symetric FIR, odd order / even length +-- Type 2: symmetric FIR, odd order / even length {-# specialize fir'2 :: Array Int Float -> Array Int Float -> [Float] -> [Float] #-} {-# specialize fir'2 :: Array Int Double -> Array Int Double -> [Double] -> [Double] #-} @@ -95,7 +95,7 @@ m = snd $ bounds h m2 = m `div` 2 --- Type 3: anti-symetric FIR, even order / odd length +-- Type 3: anti-symmetric FIR, even order / odd length {-# specialize fir'3 :: Array Int Float -> Array Int Float -> [Float] -> [Float] #-} {-# specialize fir'3 :: Array Int Double -> Array Int Double -> [Double] -> [Double] #-} @@ -111,7 +111,7 @@ m = snd $ bounds h m2 = m `div` 2 --- Type 4: anti-symetric FIR, off order / even length +-- Type 4: anti-symmetric FIR, off order / even length {-# specialize fir'4 :: Array Int Float -> Array Int Float -> [Float] -> [Float] #-} {-# specialize fir'4 :: Array Int Double -> Array Int Double -> [Double] -> [Double] #-} @@ -133,7 +133,7 @@ {-# specialize isFIRType1 :: Array Int Float -> Bool #-} {-# specialize isFIRType1 :: Array Int Double -> Bool #-} -isFIRType1 :: Num a => Array Int a -> Bool +isFIRType1 :: (Num a, Eq a) => Array Int a -> Bool isFIRType1 h = even m && (h' == (reverse h')) where m = snd $ bounds h h' = elems h @@ -141,7 +141,7 @@ {-# specialize isFIRType2 :: Array Int Float -> Bool #-} {-# specialize isFIRType2 :: Array Int Double -> Bool #-} -isFIRType2 :: Num a => Array Int a -> Bool +isFIRType2 :: (Num a, Eq a) => Array Int a -> Bool isFIRType2 h = odd m && (h' == (reverse h')) where m = snd $ bounds h h' = elems h @@ -149,7 +149,7 @@ {-# specialize isFIRType3 :: Array Int Float -> Bool #-} {-# specialize isFIRType3 :: Array Int Double -> Bool #-} -isFIRType3 :: Num a => Array Int a -> Bool +isFIRType3 :: (Num a, Eq a) => Array Int a -> Bool isFIRType3 h = even m && ha == reverse hb where m = snd $ bounds h h' = elems h @@ -160,7 +160,7 @@ {-# specialize isFIRType4 :: Array Int Float -> Bool #-} {-# specialize isFIRType4 :: Array Int Double -> Bool #-} -isFIRType4 :: Num a => Array Int a -> Bool +isFIRType4 :: (Num a, Eq a) => Array Int a -> Bool isFIRType4 h = odd m && ha == reverse hb where m = snd $ bounds h ha = elems h @@ -182,7 +182,7 @@ yt' :: [Double] yt' = fir ht xt --- This checks the symetric routines against fir'0 +-- This checks the symmetric routines against fir'0 h1 :: Array Int Double h1 = listArray (0,4) [ 1, 2, 3, 2, 1 ] diff --git a/DSP/Filter/FIR/Kaiser.hs b/DSP/Filter/FIR/Kaiser.hs --- a/DSP/Filter/FIR/Kaiser.hs +++ b/DSP/Filter/FIR/Kaiser.hs @@ -18,17 +18,17 @@ -- @Book{dsp, -- author = "Alan V. Oppenheim and Ronald W. Schafer", -- title = "Discrete-Time Signal Processing", --- publisher = "Pretice-Hall", +-- publisher = "Prentice-Hall", -- year = 1989, -- address = "Englewood Cliffs", --- series = {Pretice-Hall Signal Processing Series} +-- series = {Prentice-Hall Signal Processing Series} -- } module DSP.Filter.FIR.Kaiser (kaiser_lpf, kaiser_hpf) where import Data.Array -import DSP.Filter.FIR.Window +import DSP.Window import DSP.Filter.FIR.Taps -- Set the cutoff frequency to the middle of the transition band. This diff --git a/DSP/Filter/FIR/Sharpen.hs b/DSP/Filter/FIR/Sharpen.hs --- a/DSP/Filter/FIR/Sharpen.hs +++ b/DSP/Filter/FIR/Sharpen.hs @@ -36,7 +36,7 @@ -- | Filter shaprening routine -sharpen :: (Num a) => Array Int a -- ^ h[n] +sharpen :: (Num a, Eq a) => Array Int a -- ^ h[n] -> ([a] -> [a]) -- ^ function that implements the sharpened filter sharpen h x = step4 diff --git a/DSP/Filter/FIR/Smooth.hs b/DSP/Filter/FIR/Smooth.hs --- a/DSP/Filter/FIR/Smooth.hs +++ b/DSP/Filter/FIR/Smooth.hs @@ -42,7 +42,7 @@ where m1 (y:ys) = y : map (0.5*) ys p1 (_:ys) = map (0.5*) ys ++ [ 0, 0 ] --- Reflect makes the filter symetric (not sure where this is stated) +-- Reflect makes the filter symmetric (not sure where this is stated) reflect :: Fractional a => [a] -> [a] reflect (x:xs) = (map (0.5*) $ reverse xs) ++ x : (map (0.5*) xs) diff --git a/DSP/Filter/FIR/Taps.hs b/DSP/Filter/FIR/Taps.hs --- a/DSP/Filter/FIR/Taps.hs +++ b/DSP/Filter/FIR/Taps.hs @@ -18,10 +18,10 @@ @Book{dsp, author = "Alan V. Oppenheim and Ronald W. Schafer", title = "Discrete-Time Signal Processing", - publisher = "Pretice-Hall", + publisher = "Prentice-Hall", year = 1989, address = "Englewood Cliffs", - series = {Pretice-Hall Signal Processing Series} + series = {Prentice-Hall Signal Processing Series} } -} @@ -43,21 +43,21 @@ -- Lowpass tap function -lpf_tap :: (Integral a, Floating b) => b -> a -> b -> b +lpf_tap :: (Integral a, Floating b, Eq b) => b -> a -> b -> b lpf_tap wc m n | n-a == 0 = wc / pi | otherwise = sin (wc * (n-a)) / (pi * (n-a)) where a = (fromIntegral m) / 2 -- Highpass tap function -hpf_tap :: (Floating a1, Integral a) => a1 -> a -> a1 -> a1 +hpf_tap :: (Integral a, Floating b, Eq b) => b -> a -> b -> b hpf_tap wc m n | n-a == 0 = 1 - wc / pi | otherwise = sin (pi * (n-a)) / (pi * (n-a)) - lpf_tap wc m n where a = (fromIntegral m) / 2 -- Multiband tap function -mbf_tap :: (Floating b, Integral a) => [b] -> [b] -> a -> b -> b +mbf_tap :: (Integral a, Floating b, Eq b) => [b] -> [b] -> a -> b -> b mbf_tap (g:[]) (w:[]) m n = g * lpf_tap w m n mbf_tap (g1:g2:gs) (w:ws) m n = (g1-g2) * lpf_tap w m n + mbf_tap (g2:gs) ws m n mbf_tap _ _ _ _ = error "mbf_tap: bands out of sync" @@ -67,7 +67,7 @@ -- ws = symbol rate in normalized radians -- b = filter beta -rc_tap :: (Integral a, Floating a1) => a1 -> a1 -> a -> a1 -> a1 +rc_tap :: (Integral a, Floating b, Eq b) => b -> b -> a -> b -> b rc_tap ws b m n | n-a == 0 = 1 | den == 0 = 0 | otherwise = sin sarg / sarg * cos carg / den @@ -81,7 +81,7 @@ -- | Lowpass filter -lpf :: (Ix a, Integral a, Enum b, Floating b) => b -- ^ wc +lpf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -- ^ wc -> a -- ^ M -> Array a b -- ^ h[n] @@ -89,7 +89,7 @@ -- | Highpass filter -hpf :: (Ix a, Integral a, Enum b, Floating b) => b -- ^ wc +hpf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -- ^ wc -> a -- ^ M -> Array a b -- ^ h[n] @@ -97,7 +97,7 @@ -- | Bandpass filter -bpf :: (Ix a, Integral a, Enum b, Floating b) => b -- ^ wl +bpf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -- ^ wl -> b -- ^ wu -> a -- ^ M -> Array a b -- ^ h[n] @@ -106,7 +106,7 @@ -- | Bandstop filter -bsf :: (Ix a, Integral a, Enum b, Floating b) => b -- ^ wl +bsf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -- ^ wl -> b -- ^ wu -> a -- ^ M -> Array a b -- ^ h[n] @@ -115,7 +115,7 @@ -- | Multiband filter -mbf :: (Ix a, Integral a, Enum b, Floating b) => [b] -- ^ [mags] +mbf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => [b] -- ^ [mags] -> [b] -- ^ [w] -> a -- ^ M -> Array a b -- ^ h[n] @@ -124,7 +124,7 @@ -- | Raised-cosine filter -rc :: (Ix a, Integral a, Enum b, Floating b) => b -- ^ ws +rc :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -- ^ ws -> b -- ^ beta -> a -- ^ M -> Array a b -- ^ h[n] diff --git a/DSP/Filter/FIR/Window.hs b/DSP/Filter/FIR/Window.hs --- a/DSP/Filter/FIR/Window.hs +++ b/DSP/Filter/FIR/Window.hs @@ -1,172 +1,6 @@ ------------------------------------------------------------------------------ --- | --- Module : DSP.Filter.FIR.Window --- Copyright : (c) Matthew Donadio 1998 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- Commonly used window functions. Except for the Parzen window, the --- results of all of these /look/ right, but I have to check them against --- either Matlab or my C code. --- --- More windowing functions exist, but I have to dig through my papers to --- find the equations. --- ------------------------------------------------------------------------------ - --- TODO: These functions should probably be reworked to use list --- comprehensions... - -{- - -Reference: - -@Book{dsp, - author = "Alan V. Oppenheim and Ronald W. Schafer", - title = "Discrete-Time Signal Processing", - publisher = "Pretice-Hall", - year = 1989, - address = "Englewood Cliffs", - series = {Pretice-Hall Signal Processing Series} -} - -@Book{kay, - author = "Steven M. Kay", - title = "Modern Spectral Estimation: Theory \& Application", - publisher = "Prentice Hall", - year = 1988, - address = "Englewood Cliffs", - series = {Pretice-Hall Signal Processing Series} -} - --} - -module DSP.Filter.FIR.Window - (window, rectangular, bartlett, hanning, hamming, blackman, - kaiser, gen_hamming, parzen) where - -import DSP.Basic ((^!)) -import Data.Array - - --- | Applys a window, @w@, to a sequence @x@ - -window :: Array Int Double -- ^ w[n] - -> Array Int Double -- ^ x[n] - -> Array Int Double -- ^ w[n] * x[n] - -window w x = listArray (0,m) [ w!i * x!i | i <- [0..m] ] - where m = snd $ bounds w - --- | rectangular window - -rectangular :: Int -- ^ M - -> Array Int Double -- ^ w[n] - -rectangular m = listArray (0,m) $ repeat 1 - --- | Bartlett window - -bartlett :: Int -- ^ M - -> Array Int Double -- ^ w[n] - -bartlett = makeArray bartlett' - -bartlett' :: Double -> Double -> Double -bartlett' m n = - if n <= m / 2 - then 2 * n / m - else 2 - 2 * n / m - --- | Hanning window - -hanning :: Int -- ^ M - -> Array Int Double -- ^ w[n] - -hanning = makeArray hanning' - -hanning' :: Double -> Double -> Double -hanning' m n = 0.5 - 0.5 * cos(2 * pi * n / m) - --- | Hamming window - -hamming :: Int -- ^ M - -> Array Int Double -- ^ w[n] - -hamming = makeArray hamming' - -hamming' :: Double -> Double -> Double -hamming' m n = 0.54 - 0.46 * cos(2 * pi * n / m) - - --- | Blackman window - -blackman :: Int -- ^ M - -> Array Int Double -- ^ w[n] - -blackman = makeArray blackman' - -blackman' :: Double -> Double -> Double -blackman' m n = - 0.42 - 0.5 * cos(2 * pi * n / m) + - 0.08 * cos (4 * pi * n / m) - --- | Generalized Hamming window - -gen_hamming :: Double -- ^ alpha - -> Int -- ^ M - -> Array Int Double -- ^ w[n] - -gen_hamming = makeArray . gen_hamming' - -gen_hamming' :: Double -> Double -> Double -> Double -gen_hamming' a m n = a - (1 - a) * cos(2 * pi * n / m) - --- | rectangular window - -kaiser :: Double -- ^ beta - -> Int -- ^ M - -> Array Int Double -- ^ w[n] - -kaiser = makeArray . kaiser' - -kaiser' :: Double -> Double -> Double -> Double -kaiser' b m n = - let a = m / 2 - in i0 (b * sqrt (1 -((n-a)/a)^!2)) / i0 b - - --- Recursive computation of I0, the zeroth-order modified Bessel function --- of the first kind. - -i0 :: Double -> Double -i0 x = i0' x 2 1 - -i0' :: Double -> Double -> Double -> Double -i0' x d ds | ds < 1.0e-30 = 1 - | otherwise = ds * x^!2 / d^!2 + (i0' x (d+2) (ds * x^!2 / d^!2)) - --- I don't think this one is correct. Kay's book uses different variable --- conventions and I haven't deciphered them yet... - --- | rectangular window - -parzen :: Int -- ^ M - -> Array Int Double -- ^ w[n] - -parzen = makeArray parzen' - -parzen' :: Double -> Double -> Double -parzen' m n = - if n <= m / 2 - then 2 * (1-n/m) ^! 3 - (1-2*n/m) ^! 3 - else 2 * (1-n/m) ^! 3 +module DSP.Filter.FIR.Window {-# DEPRECATED "Use DSP.Window instead" #-} ( + module DSP.Window + ) where +import DSP.Window -makeArray :: (Double -> Double -> Double) -> Int -> Array Int Double -makeArray win m = - let md = fromIntegral m - in listArray (0,m) $ map (win md . fromIntegral) [(0::Int) ..] diff --git a/DSP/Filter/IIR/IIR.hs b/DSP/Filter/IIR/IIR.hs --- a/DSP/Filter/IIR/IIR.hs +++ b/DSP/Filter/IIR/IIR.hs @@ -35,10 +35,10 @@ @Book{dsp, author = "Alan V. Oppenheim and Ronald W. Schafer", title = "Discrete-Time Signal Processing", - publisher = "Pretice-Hall", + publisher = "Prentice-Hall", year = 1989, address = "Englewood Cliffs", - series = {Pretice-Hall Signal Processing Series} + series = {Prentice-Hall Signal Processing Series} } However, we differ in the convention of the sign of the poles, as @@ -243,7 +243,7 @@ {- specialize iir_df1 :: (Array Int Float, Array Int Float) -> [Float] -> [Float] -} {- specialize iir_df1 :: (Array Int Double, Array Int Double) -> [Double] -> [Double] -} -iir_df1 :: (Num a) => (Array Int a, Array Int a) -- ^ (b,a) +iir_df1 :: (Num a, Eq a) => (Array Int a, Array Int a) -- ^ (b,a) -> [a] -- ^ x[n] -> [a] -- ^ y[n] diff --git a/DSP/Multirate/Halfband.hs b/DSP/Multirate/Halfband.hs --- a/DSP/Multirate/Halfband.hs +++ b/DSP/Multirate/Halfband.hs @@ -27,7 +27,7 @@ -- | Halfband interpolator -hb_interp :: (Num a) => Array Int a -- ^ h[n] +hb_interp :: (Num a, Eq a) => Array Int a -- ^ h[n] -> [a] -- ^ x[n] -> [a] -- ^ y[n] @@ -39,7 +39,7 @@ -- | Halfband decimator -hb_decim :: (Num a) => Array Int a -- ^ h[n] +hb_decim :: (Num a, Eq a) => Array Int a -- ^ h[n] -> [a] -- ^ x[n] -> [a] -- ^ y[n] diff --git a/DSP/Multirate/Polyphase.hs b/DSP/Multirate/Polyphase.hs --- a/DSP/Multirate/Polyphase.hs +++ b/DSP/Multirate/Polyphase.hs @@ -29,7 +29,7 @@ -- | Polyphase interpolator -poly_interp :: Num a => Int -- ^ L +poly_interp :: (Num a, Eq a) => Int -- ^ L -> Array Int a -- ^ h[n] -> [a] -- ^ x[n] -> [a] -- ^ y[n] diff --git a/DSP/Window.hs b/DSP/Window.hs new file mode 100644 --- /dev/null +++ b/DSP/Window.hs @@ -0,0 +1,172 @@ +----------------------------------------------------------------------------- +-- | +-- Module : DSP.Window +-- Copyright : (c) Matthew Donadio 1998 +-- License : GPL +-- +-- Maintainer : m.p.donadio@ieee.org +-- Stability : experimental +-- Portability : portable +-- +-- Commonly used window functions. Except for the Parzen window, the +-- results of all of these /look/ right, but I have to check them against +-- either Matlab or my C code. +-- +-- More windowing functions exist, but I have to dig through my papers to +-- find the equations. +-- +----------------------------------------------------------------------------- + +-- TODO: These functions should probably be reworked to use list +-- comprehensions... + +{- + +Reference: + +@Book{dsp, + author = "Alan V. Oppenheim and Ronald W. Schafer", + title = "Discrete-Time Signal Processing", + publisher = "Prentice-Hall", + year = 1989, + address = "Englewood Cliffs", + series = {Prentice-Hall Signal Processing Series} +} + +@Book{kay, + author = "Steven M. Kay", + title = "Modern Spectral Estimation: Theory \& Application", + publisher = "Prentice Hall", + year = 1988, + address = "Englewood Cliffs", + series = {Prentice-Hall Signal Processing Series} +} + +-} + +module DSP.Window + (window, rectangular, bartlett, hanning, hamming, blackman, + kaiser, gen_hamming, parzen) where + +import DSP.Basic ((^!)) +import Data.Array + + +-- | Applys a window, @w@, to a sequence @x@ + +window :: Array Int Double -- ^ w[n] + -> Array Int Double -- ^ x[n] + -> Array Int Double -- ^ w[n] * x[n] + +window w x = listArray (0,m) [ w!i * x!i | i <- [0..m] ] + where m = snd $ bounds w + +-- | rectangular window + +rectangular :: Int -- ^ M + -> Array Int Double -- ^ w[n] + +rectangular m = listArray (0,m) $ repeat 1 + +-- | Bartlett window + +bartlett :: Int -- ^ M + -> Array Int Double -- ^ w[n] + +bartlett = makeArray bartlett' + +bartlett' :: Double -> Double -> Double +bartlett' m n = + if n <= m / 2 + then 2 * n / m + else 2 - 2 * n / m + +-- | Hanning window + +hanning :: Int -- ^ M + -> Array Int Double -- ^ w[n] + +hanning = makeArray hanning' + +hanning' :: Double -> Double -> Double +hanning' m n = 0.5 - 0.5 * cos(2 * pi * n / m) + +-- | Hamming window + +hamming :: Int -- ^ M + -> Array Int Double -- ^ w[n] + +hamming = makeArray hamming' + +hamming' :: Double -> Double -> Double +hamming' m n = 0.54 - 0.46 * cos(2 * pi * n / m) + + +-- | Blackman window + +blackman :: Int -- ^ M + -> Array Int Double -- ^ w[n] + +blackman = makeArray blackman' + +blackman' :: Double -> Double -> Double +blackman' m n = + 0.42 - 0.5 * cos(2 * pi * n / m) + + 0.08 * cos (4 * pi * n / m) + +-- | Generalized Hamming window + +gen_hamming :: Double -- ^ alpha + -> Int -- ^ M + -> Array Int Double -- ^ w[n] + +gen_hamming = makeArray . gen_hamming' + +gen_hamming' :: Double -> Double -> Double -> Double +gen_hamming' a m n = a - (1 - a) * cos(2 * pi * n / m) + +-- | rectangular window + +kaiser :: Double -- ^ beta + -> Int -- ^ M + -> Array Int Double -- ^ w[n] + +kaiser = makeArray . kaiser' + +kaiser' :: Double -> Double -> Double -> Double +kaiser' b m n = + let a = m / 2 + in i0 (b * sqrt (1 -((n-a)/a)^!2)) / i0 b + + +-- Recursive computation of I0, the zeroth-order modified Bessel function +-- of the first kind. + +i0 :: Double -> Double +i0 x = i0' x 2 1 + +i0' :: Double -> Double -> Double -> Double +i0' x d ds | ds < 1.0e-30 = 1 + | otherwise = ds * x^!2 / d^!2 + (i0' x (d+2) (ds * x^!2 / d^!2)) + +-- I don't think this one is correct. Kay's book uses different variable +-- conventions and I haven't deciphered them yet... + +-- | rectangular window + +parzen :: Int -- ^ M + -> Array Int Double -- ^ w[n] + +parzen = makeArray parzen' + +parzen' :: Double -> Double -> Double +parzen' m n = + if n <= m / 2 + then 2 * (1-n/m) ^! 3 - (1-2*n/m) ^! 3 + else 2 * (1-n/m) ^! 3 + + +makeArray :: (Double -> Double -> Double) -> Int -> Array Int Double +makeArray win m = + let md = fromIntegral m + in listArray (0,m) $ map (win md . fromIntegral) [(0::Int) ..] diff --git a/Makefile b/Makefile deleted file mode 100644 --- a/Makefile +++ /dev/null @@ -1,117 +0,0 @@ -LIBS= Matrix -LIBS+= Polynomial -LIBS+= Numeric -LIBS+= DSP - -APPS= FFTBench -APPS+= FFTTest -APPS+= Article -APPS+= IIRDemo -APPS+= FreqDemo -APPS+= NoiseDemo - -#OPT= -O -funbox-strict-fields -#OPT+= -fvia-C -O2-for-C - -#PROF= -prof -auto-all - -#WARN= -W -#WARN= -Wall - -#HEAP= -H128m - -HC= ghc-5.04.3 -#HC= ghc-6.0 - -GC= green-card - -GC_PATH= /usr/local/lib/green-card -GC_LIBS= $(GC_PATH)/StdDIS.o - -GSL_PATH= /usr/local/lib -GSL_INC= /usr/local/include -#GSL_LIBS= -lgsl -lgslcblas -GSL_LIBS= -lgsl -lcblas -latlas - -TARGET= ffi - -GCSRCS= Numeric/Special/Airy.gc -GCSRCS+= Numeric/Special/Bessel.gc -GCSRCS+= Numeric/Special/Clausen.gc -GCSRCS+= Numeric/Special/Ellint.gc -GCSRCS+= Numeric/Special/Elljac.gc -GCSRCS+= Numeric/Special/Erf.gc - -GCOBJS= $(GCSRCS:.gc=.o) - -HSFLAGS= -fno-glasgow-exts ${OPT} ${PROF} ${WARN} ${HEAP} -HSFLAGS+= -cpp -fffi -package lang -I$(GSL_INC) -i$(GC_PATH) - -INSTALLDIR= /usr/home/donadio/lib/hs - -URL= http://haskelldsp.sourceforge.net/ - -.SUFFIXES: -.SUFFIXES: .o .hs .gc - -.gc.o: - $(GC) -t $(TARGET) -i $(GC_PATH) $< - $(HC) $(HSFLAGS) -I. -I$(GSL_INC) -i$(GC_PATH) -package-name=Numeric -package lang -c $*.hs -o $*_hs.o - $(HC) $(HSFLAGS) -I. -I$(GSL_INC) -i$(GC_PATH) -package-name=Numeric -package lang -c $*_stub_$(TARGET).c -o $*_stub_$(TARGET).o - $(LD) -r -o $@ $*_hs.o $*_stub_$(TARGET).o -# $(RM) $*.hs - $(RM) $*_stub_$(TARGET).c $*_stub_$(TARGET).h - $(RM) $*_hs.o $*_stub_$(TARGET).o - -all: foreign libs - -apps: - for f in ${APPS}; do \ - $(HC) --make ${HSFLAGS} -o $$f demo/$$f.hs ;\ - done - -foreign: $(GCOBJS) - -libs: - for f in ${LIBS}; do \ - find $$f \( -name "*.lhs" -or -name "*.hs" \) -print | xargs $(HC) --make -package-name=$$f ${HSFLAGS} ;\ - done - -docs: - find ${LIBS} -name "*.hs" -print | xargs haddock -h -o doc -t "Haskell DSP Library" -s "${URL}" - -install: #libs - for f in ${LIBS}; do \ - rm -f libHS$$f.a HS$$f.o ;\ - find $$f -name "*.o" -print | xargs ar cqs libHS$$f.a ;\ - ld -r --whole-archive -o HS$$f.o libHS$$f.a ;\ - for i in `find $$f -name "*.hi" -print`; do \ - mkdir -p `dirname ${INSTALLDIR}/imports/$$i` ;\ - install -m 644 $$i ${INSTALLDIR}/imports/$$i ;\ - done ;\ - install -m 644 libHS$$f.a ${INSTALLDIR} ;\ - install -m 644 HS$$f.o ${INSTALLDIR} ;\ - installdir=${INSTALLDIR} ghc-pkg -f mpd.conf -u < pkg/$$f.pkg ;\ - done - -test: test.hs - $(HC) -cpp -package lang -i$(GC_PATH) -L$(GSL_PATH) $(GSL_LIBS) --make -o test test.hs - -snapshot: - tar cfz haskelldsp-snapshot.tar.gz ${LIBS} demo doc \ - COPYING README TODO Makefile - -clean: - find . -name "*.o" -print | xargs rm -f - find . -name "*.hi" -print | xargs rm -f - find . -name "*~" -print | xargs rm -f - rm -f $(GCSRCS:.gc=.hs) - rm -f $(GCSRCS:.gc=_stub_$(TARGET).c) - rm -f $(GCSRCS:.gc=_stub_$(TARGET).h) - rm -f *.a - rm -f *.prof - -realclean: clean - rm -f ${APPS} - rm -f *.core - rm -f haskelldsp-snapshot.tar.gz diff --git a/Numeric/Special/Airy.gc b/Numeric/Special/Airy.gc deleted file mode 100644 --- a/Numeric/Special/Airy.gc +++ /dev/null @@ -1,276 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Airy --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Airy functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Airy (airy_Ai, airy_Ai_e, - airy_Ai_scaled, airy_Ai_scaled_e, - airy_Ai_deriv, airy_Ai_deriv_e, - airy_Ai_deriv_scaled, airy_Ai_deriv_scaled_e, - airy_zero_Ai, airy_zero_Ai_e, - airy_zero_Ai_deriv, airy_zero_Ai_deriv_e, - airy_Bi, airy_Bi_e, - airy_Bi_scaled, airy_Bi_scaled_e, - airy_Bi_deriv, airy_Bi_deriv_e, - airy_Bi_deriv_scaled, airy_Bi_deriv_scaled_e, - airy_zero_Bi, airy_zero_Bi_e, - airy_zero_Bi_deriv, airy_zero_Bi_deriv_e - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun airy_Ai :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Ai(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Ai_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Ai_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Ai_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Ai_scaled(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Ai_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Ai_scaled_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Ai_deriv :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Ai_deriv(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Ai_deriv_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Ai_deriv_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Ai_deriv_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Ai_deriv_scaled(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Ai_deriv_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Ai_deriv_scaled_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_zero_Ai :: Int -> Double -%call (int s) -%code double z; -% z = gsl_sf_airy_zero_Ai(s); -%result (double z) - -%fun airy_zero_Ai_e :: Int -> (Double, Double) -%call (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_zero_Ai_e(s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_zero_Ai_deriv :: Int -> Double -%call (int s) -%code double z; -% z = gsl_sf_airy_zero_Ai_deriv(s); -%result (double z) - -%fun airy_zero_Ai_deriv_e :: Int -> (Double, Double) -%call (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_zero_Ai_deriv_e(s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Bi :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Bi(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Bi_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Bi_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Bi_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Bi_scaled(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Bi_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Bi_scaled_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Bi_deriv :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Bi_deriv(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Bi_deriv_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Bi_deriv_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_Bi_deriv_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_airy_Bi_deriv_scaled(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun airy_Bi_deriv_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_Bi_deriv_scaled_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_zero_Bi :: Int -> Double -%call (int s) -%code double z; -% z = gsl_sf_airy_zero_Bi(s); -%result (double z) - -%fun airy_zero_Bi_e :: Int -> (Double, Double) -%call (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_zero_Bi_e(s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun airy_zero_Bi_deriv :: Int -> Double -%call (int s) -%code double z; -% z = gsl_sf_airy_zero_Bi_deriv(s); -%result (double z) - -%fun airy_zero_Bi_deriv_e :: Int -> (Double, Double) -%call (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_airy_zero_Bi_deriv_e(s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) diff --git a/Numeric/Special/Bessel.gc b/Numeric/Special/Bessel.gc deleted file mode 100644 --- a/Numeric/Special/Bessel.gc +++ /dev/null @@ -1,874 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Bessel --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Bessel functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Bessel (bessel_J0, bessel_J0_e, - bessel_J1, bessel_J1_e, - bessel_Jn, bessel_Jn_e, - - bessel_Y0, bessel_Y0_e, - bessel_Y1, bessel_Y1_e, - bessel_Yn, bessel_Yn_e, - - bessel_I0, bessel_I0_e, - bessel_I1, bessel_I1_e, - bessel_In, bessel_In_e, - - bessel_I0_scaled, bessel_I0_scaled_e, - bessel_I1_scaled, bessel_I1_scaled_e, - bessel_In_scaled, bessel_In_scaled_e, - - bessel_K0, bessel_K0_e, - bessel_K1, bessel_K1_e, - bessel_Kn, bessel_Kn_e, - - bessel_K0_scaled, bessel_K0_scaled_e, - bessel_K1_scaled, bessel_K1_scaled_e, - bessel_Kn_scaled, bessel_Kn_scaled_e, - - bessel_j0, bessel_j0_e, - bessel_j1, bessel_j1_e, - bessel_jl, bessel_jl_e, - - bessel_y0, bessel_y0_e, - bessel_y1, bessel_y1_e, - bessel_yl, bessel_yl_e, - - bessel_i0_scaled, bessel_i0_scaled_e, - bessel_i1_scaled, bessel_i1_scaled_e, - bessel_il_scaled, bessel_il_scaled_e, - - bessel_k0_scaled, bessel_k0_scaled_e, - bessel_k1_scaled, bessel_k1_scaled_e, - bessel_kl_scaled, bessel_kl_scaled_e, - - bessel_Jnu, bessel_Jnu_e, - bessel_Ynu, bessel_Ynu_e, - bessel_Inu, bessel_Inu_e, - bessel_Inu_scaled, bessel_Inu_scaled_e, - bessel_Knu, bessel_Knu_e, - bessel_Knu_scaled, bessel_Knu_scaled_e, - bessel_lnKnu, bessel_lnKnu_e, - - bessel_zero_J0, bessel_zero_J0_e, - bessel_zero_J1, bessel_zero_J1_e, - bessel_zero_Jnu, bessel_zero_Jnu_e - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun bessel_J0 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_J0(x); -%result (double y) - -%fun bessel_J0_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_J0_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_J1 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_J1(x); -%result (double y) - -%fun bessel_J1_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_J1_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Jn :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_Jn(n, x); -%result (double y) - -%fun bessel_Jn_e :: Int -> Double -> (Double, Double) -%call (int n) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Jn_e(n, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Y0 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_Y0(x); -%result (double y) - -%fun bessel_Y0_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Y0_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Y1 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_Y1(x); -%result (double y) - -%fun bessel_Y1_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Y1_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Yn :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_Yn(n, x); -%result (double y) - -%fun bessel_Yn_e :: Int -> Double -> (Double, Double) -%call (int n) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Yn_e(n, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_I0 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_I0(x); -%result (double y) - -%fun bessel_I0_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_I0_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_I1 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_I1(x); -%result (double y) - -%fun bessel_I1_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_I1_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_In :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_In(n, x); -%result (double y) - -%fun bessel_In_e :: Int -> Double -> (Double, Double) -%call (int n) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_In_e(n, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_I0_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_I0_scaled(x); -%result (double y) - -%fun bessel_I0_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_I0_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_I1_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_I1_scaled(x); -%result (double y) - -%fun bessel_I1_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_I1_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_In_scaled :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_In_scaled(n, x); -%result (double y) - -%fun bessel_In_scaled_e :: Int -> Double -> (Double, Double) -%call (int n) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_In_scaled_e(n, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_K0 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_K0(x); -%result (double y) - -%fun bessel_K0_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_K0_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_K1 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_K1(x); -%result (double y) - -%fun bessel_K1_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_K1_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Kn :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_Kn(n, x); -%result (double y) - -%fun bessel_Kn_e :: Int -> Double -> (Double, Double) -%call (int n) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Kn_e(n, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) -------------------------------------------------------------------------------- - -%fun bessel_K0_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_K0_scaled(x); -%result (double y) - -%fun bessel_K0_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_K0_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_K1_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_K1_scaled(x); -%result (double y) - -%fun bessel_K1_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_K1_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Kn_scaled :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_Kn_scaled(n, x); -%result (double y) - -%fun bessel_Kn_scaled_e :: Int -> Double -> (Double, Double) -%call (int n) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Kn_scaled_e(n, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_j0 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_j0(x); -%result (double y) - -%fun bessel_j0_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_j0_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_j1 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_j1(x); -%result (double y) - -%fun bessel_j1_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_j1_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_jl :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_jl(n, x); -%result (double y) - -%fun bessel_jl_e :: Int -> Double -> (Double, Double) -%call (int l) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_jl_e(l, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_y0 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_y0(x); -%result (double y) - -%fun bessel_y0_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_y0_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_y1 :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_y1(x); -%result (double y) - -%fun bessel_y1_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_y1_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_yl :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_yl(n, x); -%result (double y) - -%fun bessel_yl_e :: Int -> Double -> (Double, Double) -%call (int l) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_yl_e(l, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_i0_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_i0_scaled(x); -%result (double y) - -%fun bessel_i0_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_i0_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_i1_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_i1_scaled(x); -%result (double y) - -%fun bessel_i1_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_i1_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_il_scaled :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_il_scaled(n, x); -%result (double y) - -%fun bessel_il_scaled_e :: Int -> Double -> (Double, Double) -%call (int l) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_il_scaled_e(l, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_k0_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_k0_scaled(x); -%result (double y) - -%fun bessel_k0_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_k0_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_k1_scaled :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_bessel_k1_scaled(x); -%result (double y) - -%fun bessel_k1_scaled_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_k1_scaled_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_kl_scaled :: Int -> Double -> Double -%call (int n) (double x) -%code double y; -% y = gsl_sf_bessel_kl_scaled(n, x); -%result (double y) - -%fun bessel_kl_scaled_e :: Int -> Double -> (Double, Double) -%call (int l) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_kl_scaled_e(l, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Jnu :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_Jnu(nu, x); -%result (double y) - -%fun bessel_Jnu_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Jnu_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Ynu :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_Ynu(nu, x); -%result (double y) - -%fun bessel_Ynu_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Ynu_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Inu :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_Inu(nu, x); -%result (double y) - -%fun bessel_Inu_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Inu_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Inu_scaled :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_Inu_scaled(nu, x); -%result (double y) - -%fun bessel_Inu_scaled_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Inu_scaled_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Knu :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_Knu(nu, x); -%result (double y) - -%fun bessel_Knu_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Knu_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_Knu_scaled :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_Knu_scaled(nu, x); -%result (double y) - -%fun bessel_Knu_scaled_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_Knu_scaled_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_lnKnu :: Double -> Double -> Double -%call (double nu) (double x) -%code double y; -% y = gsl_sf_bessel_lnKnu(nu, x); -%result (double y) - -%fun bessel_lnKnu_e :: Double -> Double -> (Double, Double) -%call (double nu) (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_lnKnu_e(nu, x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_zero_J0 :: Int -> Double -%call (int s) -%code double y; -% y = gsl_sf_bessel_zero_J0(s); -%result (double y) - -%fun bessel_zero_J0_e :: Int -> (Double, Double) -%call (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_zero_J0_e(s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_zero_J1 :: Int -> Double -%call (int s) -%code double y; -% y = gsl_sf_bessel_zero_J1(s); -%result (double y) - -%fun bessel_zero_J1_e :: Int -> (Double, Double) -%call (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_zero_J1_e(s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun bessel_zero_Jnu :: Double -> Int -> Double -%call (double nu) (int s) -%code double y; -% y = gsl_sf_bessel_zero_Jnu(nu, s); -%result (double y) - -%fun bessel_zero_Jnu_e :: Double -> Int -> (Double, Double) -%call (double nu) (int s) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_bessel_zero_Jnu_e(nu, s, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) diff --git a/Numeric/Special/Clausen.gc b/Numeric/Special/Clausen.gc deleted file mode 100644 --- a/Numeric/Special/Clausen.gc +++ /dev/null @@ -1,45 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Clausen --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Clausen functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Clausen (clausen, clausen_e, - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun clausen :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_clausen(x); -%result (double y) - -%fun clausen_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_clausen_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) diff --git a/Numeric/Special/Ellint.gc b/Numeric/Special/Ellint.gc deleted file mode 100644 --- a/Numeric/Special/Ellint.gc +++ /dev/null @@ -1,234 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Ellint --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Ellint functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Ellint (ellint_Kcomp, ellint_Kcomp_e, - ellint_Ecomp, ellint_Ecomp_e, - ellint_F, ellint_F_e, - ellint_E, ellint_E_e, - ellint_P, ellint_P_e, - ellint_D, ellint_D_e, - ellint_RC, ellint_RC_e, - ellint_RD, ellint_RD_e, - ellint_RF, ellint_RF_e, - ellint_RJ, ellint_RJ_e - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun ellint_Kcomp :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_ellint_Kcomp(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun ellint_Kcomp_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_Kcomp_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_Ecomp :: Double -> Double -%call (double k) -%code double y; -% y = gsl_sf_ellint_Ecomp(k, GSL_PREC_DOUBLE); -%result (double y) - -%fun ellint_Ecomp_e :: Double -> (Double, Double) -%call (double k) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_Ecomp_e(k, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_F :: Double -> Double -> Double -%call (double phi) (double k) -%code double y; -% y = gsl_sf_ellint_F(phi, k, GSL_PREC_DOUBLE); -%result (double y) - -%fun ellint_F_e :: Double -> Double -> (Double, Double) -%call (double phi) (double k) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_F_e(phi, k, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_E :: Double -> Double -> Double -%call (double phi) (double k) -%code double y; -% y = gsl_sf_ellint_E(phi, k, GSL_PREC_DOUBLE); -%result (double y) - -%fun ellint_E_e :: Double -> Double -> (Double, Double) -%call (double phi) (double k) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_E_e(phi, k, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_P :: Double -> Double -> Double -> Double -%call (double phi) (double k) (double n) -%code double y; -% y = gsl_sf_ellint_P(phi, k, n, GSL_PREC_DOUBLE); -%result (double y) - -%fun ellint_P_e :: Double -> Double -> Double -> (Double, Double) -%call (double phi) (double k) (double n) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_P_e(phi, k, n, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_D :: Double -> Double -> Double -> Double -%call (double phi) (double k) (double n) -%code double y; -% y = gsl_sf_ellint_D(phi, k, n, GSL_PREC_DOUBLE); -%result (double y) - -%fun ellint_D_e :: Double -> Double -> Double -> (Double, Double) -%call (double phi) (double k) (double n) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_D_e(phi, k, n, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_RC :: Double -> Double -> Double -%call (double x) (double y) -%code double it; -% it = gsl_sf_ellint_RC(x, y, GSL_PREC_DOUBLE); -%result (double it) - -%fun ellint_RC_e :: Double -> Double -> (Double, Double) -%call (double x) (double y) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_RC_e(x, y, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_RD :: Double -> Double -> Double -> Double -%call (double x) (double y) (double z) -%code double it; -% it = gsl_sf_ellint_RD(x, y, z, GSL_PREC_DOUBLE); -%result (double it) - -%fun ellint_RD_e :: Double -> Double -> Double -> (Double, Double) -%call (double x) (double y) (double z) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_RD_e(x, y, z, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_RF :: Double -> Double -> Double -> Double -%call (double x) (double y) (double z) -%code double it; -% it = gsl_sf_ellint_RF(x, y, z, GSL_PREC_DOUBLE); -%result (double it) - -%fun ellint_RF_e :: Double -> Double -> Double -> (Double, Double) -%call (double x) (double y) (double z) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_RF_e(x, y, z, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun ellint_RJ :: Double -> Double -> Double -> Double -> Double -%call (double x) (double y) (double z) (double p) -%code double it; -% it = gsl_sf_ellint_RJ(x, y, z, p, GSL_PREC_DOUBLE); -%result (double it) - -%fun ellint_RJ_e :: Double -> Double -> Double -> Double -> (Double, Double) -%call (double x) (double y) (double z) (double p) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_ellint_RJ_e(x, y, z, p, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) diff --git a/Numeric/Special/Elljac.gc b/Numeric/Special/Elljac.gc deleted file mode 100644 --- a/Numeric/Special/Elljac.gc +++ /dev/null @@ -1,93 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Elljac --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Elljac functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Elljac (elljac_e, - elljac_sn_e, elljac_cn_e, elljac_dn_e, - elljac_cd_e, elljac_dc_e, elljac_ns_e, - elljac_sd_e, elljac_nc_e, elljac_ds_e, - elljac_nd_e, elljac_sc_e, elljac_cs_e - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun elljac_e :: Double -> Double -> (Double, Double, Double) -%call (double u) (double m) -%code int status; -% double sn; -% double cn; -% double dn; -% status = gsl_sf_elljac_e(u, m, &sn, &cn, &dn); -%fail {status != 0} {gsl_strerror(status)} -%result (double sn, double cn, double dn) - -------------------------------------------------------------------------------- - --- Abramowitz & Stegun, Sec 16.3 - -elljac_sn_e :: Double -> Double -> Double -elljac_sn_e u m = sn - where (sn,_,_) = elljac_e u m - -elljac_cn_e :: Double -> Double -> Double -elljac_cn_e u m = cn - where (_,cn,_) = elljac_e u m - -elljac_dn_e :: Double -> Double -> Double -elljac_dn_e u m = dn - where (_,_,dn) = elljac_e u m - -elljac_cd_e :: Double -> Double -> Double -elljac_cd_e u m = cn / dn - where (_,cn,dn) = elljac_e u m - -elljac_sd_e :: Double -> Double -> Double -elljac_sd_e u m = sn / dn - where (sn,_,dn) = elljac_e u m - -elljac_nd_e :: Double -> Double -> Double -elljac_nd_e u m = 1 / dn - where (_,_,dn) = elljac_e u m - -elljac_dc_e :: Double -> Double -> Double -elljac_dc_e u m = dn / cn - where (_,cn,dn) = elljac_e u m - -elljac_nc_e :: Double -> Double -> Double -elljac_nc_e u m = 1 / cn - where (_,cn,_) = elljac_e u m - -elljac_sc_e :: Double -> Double -> Double -elljac_sc_e u m = sn / cn - where (sn,cn,_) = elljac_e u m - -elljac_ns_e :: Double -> Double -> Double -elljac_ns_e u m = 1 / sn - where (sn,_,_) = elljac_e u m - -elljac_ds_e :: Double -> Double -> Double -elljac_ds_e u m = dn / sn - where (sn,_,dn) = elljac_e u m - -elljac_cs_e :: Double -> Double -> Double -elljac_cs_e u m = cn / sn - where (sn,cn,_) = elljac_e u m diff --git a/Numeric/Special/Erf.gc b/Numeric/Special/Erf.gc deleted file mode 100644 --- a/Numeric/Special/Erf.gc +++ /dev/null @@ -1,129 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Erf --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Erf functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Erf (erfc, erfc_e, - log_erfc, log_erfc_e, - erf, erf_e, - erf_Z, erf_Z_e, - erf_Q, erf_Q_e, - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun erfc :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_erfc(x); -%result (double y) - -%fun erfc_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_erfc_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun log_erfc :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_log_erfc(x); -%result (double y) - -%fun log_erfc_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_log_erfc_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun erf :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_erf(x); -%result (double y) - -%fun erf_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_erf_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun erf_Z :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_erf_Z(x); -%result (double y) - -%fun erf_Z_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_erf_Z_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) - -------------------------------------------------------------------------------- - -%fun erf_Q :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_erf_Q(x); -%result (double y) - -%fun erf_Q_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_erf_Q_e(x, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) diff --git a/Numeric/Special/Foo.gc b/Numeric/Special/Foo.gc deleted file mode 100644 --- a/Numeric/Special/Foo.gc +++ /dev/null @@ -1,45 +0,0 @@ -{-# OPTIONS -fffi -fvia-C #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Special.Foo --- Copyright : (c) Matthew Donadio 2003 --- License : GPL --- --- Maintainer : m.p.donadio@ieee.org --- Stability : experimental --- Portability : portable --- --- FFI to GSL for the Foo functions --- ------------------------------------------------------------------------------ - -module Numeric.Special.Foo (foo_Bar, foo_Bar_e, - ) where - -import StdDIS - -import Foreign - -%#include -%#include - -------------------------------------------------------------------------------- - -%fun foo_Bar :: Double -> Double -%call (double x) -%code double y; -% y = gsl_sf_foo_Bar(x, GSL_PREC_DOUBLE); -%result (double y) - -%fun foo_Bar_e :: Double -> (Double, Double) -%call (double x) -%code int status; -% double val; -% double err; -% gsl_sf_result result; -% status = gsl_sf_foo_Bar_e(x, GSL_PREC_DOUBLE, &result); -% val = result.val; -% err = result.err; -%fail {status != 0} {gsl_strerror(status)} -%result (double val, double err) diff --git a/Numeric/Statistics/Covariance.hs b/Numeric/Statistics/Covariance.hs --- a/Numeric/Statistics/Covariance.hs +++ b/Numeric/Statistics/Covariance.hs @@ -20,8 +20,6 @@ module Numeric.Statistics.Covariance (cov) where -import Data.List - import Numeric.Statistics.Moment cov :: (Fractional a) => [a] -> [a] -> a diff --git a/Numeric/Statistics/Median.hs b/Numeric/Statistics/Median.hs --- a/Numeric/Statistics/Median.hs +++ b/Numeric/Statistics/Median.hs @@ -16,7 +16,7 @@ module Numeric.Statistics.Median (median, medianFast) where -import Data.List +import Data.List (sort) -- | Compute the median of a list diff --git a/Numeric/Statistics/Moment.hs b/Numeric/Statistics/Moment.hs --- a/Numeric/Statistics/Moment.hs +++ b/Numeric/Statistics/Moment.hs @@ -24,8 +24,6 @@ -- TODO: does var passes through the list twice, once to compute the mean of -- the squares, and the other to compute the mean? -import Data.List - -- * Functions -- | Compute the mean of a list diff --git a/Numeric/Statistics/TTest.hs b/Numeric/Statistics/TTest.hs --- a/Numeric/Statistics/TTest.hs +++ b/Numeric/Statistics/TTest.hs @@ -18,8 +18,6 @@ module Numeric.Statistics.TTest (ttest, tutest, tptest) where -import Data.List - import Numeric.Statistics.Covariance import Numeric.Statistics.Moment diff --git a/Numeric/Transform/Fourier/CT.hs b/Numeric/Transform/Fourier/CT.hs --- a/Numeric/Transform/Fourier/CT.hs +++ b/Numeric/Transform/Fourier/CT.hs @@ -14,7 +14,7 @@ module Numeric.Transform.Fourier.CT (fft_ct1, fft_ct2) where -import Data.List +import Data.List (transpose) import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/Eigensystem.hs b/Numeric/Transform/Fourier/Eigensystem.hs new file mode 100644 --- /dev/null +++ b/Numeric/Transform/Fourier/Eigensystem.hs @@ -0,0 +1,142 @@ +module Numeric.Transform.Fourier.Eigensystem where + +import qualified Numeric.Transform.Fourier.FFT as FT +import qualified Data.Complex as Complex +import qualified Data.Array as Array +import Data.Array (Array, Ix, ) + + +(^!) :: Num a => a -> Int -> a +(^!) = (^) + +{-# specialize gaussianPlain :: Int -> Float -> Array Int Float #-} +{-# specialize gaussianPlain :: Int -> Double -> Array Int Double #-} + +{- | +For small values of @narrow@ (large width) +it is faster to compute the spectrum of the Gaussian +and apply the Fourier transform on it. +-} +gaussianFourier :: + (Ix a, Integral a, RealFloat b) => + a -- ^ array size + -> b -- ^ reciprocal of width + -> Array a b -- ^ X[k] +gaussianFourier n narrow = + fmap ((/ (narrow*sqrt (fromIntegral n))) . Complex.realPart) $ + FT.rfft $ gaussianPlain n (recip narrow) + +gaussianPlain :: + (Ix a, Integral a, RealFloat b) => + a -- ^ array size + -> b -- ^ reciprocal of width, 1 is the width of the eigenvector + -> Array a b -- ^ X[k] +gaussianPlain n narrow = + let nr = fromIntegral n + in Array.listArray (0,n-1) $ + map (wrappedGaussian 1e-20 (narrow^!2 * nr * pi)) $ + map (/nr) $ + iterate (1+) 0 + +{- | +For small @c@ the convergence is slow, +but the result will be close to @recip (sqrt c)@. +-} +wrappedGaussian :: + (RealFloat b) => + b -> b -> b -> b +wrappedGaussian tol c x = + let xpos = iterate (1+) x + xneg = tail $ iterate (subtract 1) x + applyTrunc = + takeWhile ((>=tol) . abs) . + map (\xi -> exp (- c*xi*xi)) + in sum (applyTrunc xpos) + + sum (applyTrunc xneg) + + +example0 :: ([Double], [Complex.Complex Double]) +example0 = + let c=0.06; n=13 in (map (\k -> (sqrt (pi/c) *) $ sum $ map (\l -> exp (-pi^!2/(n^!2*c)*(k+n*l)^!2)) [-100..100]) [0..n-1], map (\j -> sum $ map (\k -> exp (((-c*k^!2) Complex.:+ 2*pi*j*k/n))) [-100..100]) [0..n-1]) + +example1 :: (Double, Double) +example1 = + let c=0.3 in (sum $ map (\k -> exp(-pi*k^!2/c)) [-100..100], sqrt c * sum (map (\k -> exp(-pi*k^!2*c)) [-100..100])) + +exampleAlgebraicDouble :: (Double, Double) +exampleAlgebraicDouble = + (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sqrt (4 - sqrt 8) * sum (map (\k -> exp(-pi*2*k^!2)) [-100..100])) + +{- +I checked polynomials up to degree 9 +and there does not seem to be one, +that has the ratio of the numbers as root. + +exampleAlgebraicTriple :: (Double, Double) +exampleAlgebraicTriple = + (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sum (map (\k -> exp(-pi*3*k^!2)) [-100..100])) +-} + +exampleAlgebraicQuadro :: (Double, Double) +exampleAlgebraicQuadro = + let a = (2-sqrt(2*sqrt 2)) * (2 + sqrt 2) + _b = 4 / (2 + sqrt(sqrt 8)) :: Double + in (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], a * sum (map (\k -> exp(-pi*4*k^!2)) [-100..100])) + +exampleAlgebraicHalf :: (Double, Double) +exampleAlgebraicHalf = + (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], sqrt (2 - sqrt 2) * sum (map (\k -> exp(-pi/2*k^!2)) [-100..100])) + +exampleAlgebraicFourth :: (Double, Double) +exampleAlgebraicFourth = + let a = (2-sqrt(2*sqrt 2)) * (2+sqrt 2)/2 + _b = 2/(2+sqrt(sqrt 8)) :: Double + in (sum $ map (\k -> exp(-pi*k^!2)) [-100..100], a * sum (map (\k -> exp(-pi/4*k^!2)) [-100..100])) + +exampleAlgebraic2 :: (Double, Double) +exampleAlgebraic2 = + let n=2 in (sum $ map (\k -> exp(-pi*n*k^!2)) [-100..100], (1+sqrt 2) * sum (map (\k -> exp(-pi*n*(k+1/n)^!2)) [-100..100])) +{- +a = 1+sqrt 2 +a - 1 = sqrt 2 +(a - 1)^2 = 2 +0 = -1 - 2*a + a^2 +-} + +exampleAlgebraic3 :: (Double, Double) +exampleAlgebraic3 = + let n=3 in (sum $ map (\k -> exp(-pi*n*k^!2)) [-100..100], (1+sqrt 3) * sum (map (\k -> exp(-pi*n*(k+1/n)^!2)) [-100..100])) +{- +0 = -2 - 2*a + a^2 +-} + +exampleAlgebraic4 :: [Double] +exampleAlgebraic4 = + let n=4 in zipWith (*) [1, 1+sqrt(sqrt 2), (1+sqrt(sqrt 2))/(sqrt(sqrt 2)-1), 1+sqrt(sqrt 2)] $ map (\j -> sum (map (\k -> exp(-pi*n*(k+j/n)^!2)) [-100..100])) [0..n-1] +{- +0 = -1 - 4*a + 6*a^2 - 4*a^3 + a^4 +0 = 1 - 12*b + 6*b^2 - 12*b^3 + b^4 +-} + +exampleAlgebraic5 :: [Double] +exampleAlgebraic5 = + let n=5; a=1.8743053964841243; b=11.833898536015248 in zipWith (*) [1, a, b, a, 1] $ map (\j -> sum (map (\k -> exp(-pi*n*(k+j/n)^!2)) [-100..100])) [0..n-1] +{- +This must be a linear combination of +(goldenRatio, 1, 0, 0, 1) +(goldenRatio, 0, 1, 1, 0) + +with goldenRatio = (1+sqrt 5)/2 + +0=-4*a^0-4*a^1+26*a^2-14*a^3+a^4 +0=-4*b^0-4*b^1+26*b^2-14*b^3+b^4 +-} + + + +{- +ToDo: + - efficient computation of 2-norm of the discrete Gaussian + - investigate eigenfunctions including Gaussian and Hermite polynomials + - is there an algebraic ratio between exampleAlgebraicTriple and exampleAlgebraicSixtimes +-} diff --git a/Numeric/Transform/Fourier/FFT.hs b/Numeric/Transform/Fourier/FFT.hs --- a/Numeric/Transform/Fourier/FFT.hs +++ b/Numeric/Transform/Fourier/FFT.hs @@ -16,7 +16,6 @@ module Numeric.Transform.Fourier.FFT (fft, ifft, rfft, irfft, r2fft) where -import Data.List import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/FFTUtils.hs b/Numeric/Transform/Fourier/FFTUtils.hs --- a/Numeric/Transform/Fourier/FFTUtils.hs +++ b/Numeric/Transform/Fourier/FFTUtils.hs @@ -28,7 +28,7 @@ magsq :: RealFloat a => Complex a -> a magsq (x:+y) = x*x + y*y -log10 :: Floating a => a -> a +log10 :: (Floating a, Eq a) => a -> a log10 0 = -1.0e9 log10 x = logBase 10 x @@ -111,12 +111,11 @@ FilePath -> Array i e -> IO () write_cvector f x = do let n = (snd $ bounds x) + 1 - h <- openFile f WriteMode - sequence $ map (hPrintIndex h n) $ assocs $ x - hClose h + withFile f WriteMode $ \h -> + mapM_ (hPrintIndex h n) $ assocs x write_fft_info :: (Ix i, Integral i) => - [Char] -> Array i (Complex Double) -> IO () + String -> Array i (Complex Double) -> IO () write_fft_info b x = do let (mag,db,arg,grd) = fft_info x write_cvector (b ++ "_mag.out") mag @@ -127,11 +126,10 @@ write_rvector :: Show e => FilePath -> Array Int e -> IO () write_rvector f x = do let n = (snd $ bounds x) + 1 - h <- openFile f WriteMode - sequence $ map (hPrintIndex h n) $ take (n `div` 2) $ assocs $ x - hClose h + withFile f WriteMode $ \h -> + mapM_ (hPrintIndex h n) $ take (n `div` 2) $ assocs x -write_rfft_info :: [Char] -> Array Int Double -> IO () +write_rfft_info :: String -> Array Int Double -> IO () write_rfft_info b x = do let (mag,db,arg,grd) = rfft_info x write_rvector (b ++ "_mag.out") mag diff --git a/Numeric/Transform/Fourier/PFA.hs b/Numeric/Transform/Fourier/PFA.hs --- a/Numeric/Transform/Fourier/PFA.hs +++ b/Numeric/Transform/Fourier/PFA.hs @@ -14,7 +14,7 @@ module Numeric.Transform.Fourier.PFA (fft_pfa) where -import Data.List +import Data.List (transpose) import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/R2DIF.hs b/Numeric/Transform/Fourier/R2DIF.hs --- a/Numeric/Transform/Fourier/R2DIF.hs +++ b/Numeric/Transform/Fourier/R2DIF.hs @@ -15,7 +15,6 @@ module Numeric.Transform.Fourier.R2DIF (fft_r2dif) where import DSP.Basic (interleave) -import Data.List import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/R2DIT.hs b/Numeric/Transform/Fourier/R2DIT.hs --- a/Numeric/Transform/Fourier/R2DIT.hs +++ b/Numeric/Transform/Fourier/R2DIT.hs @@ -14,7 +14,6 @@ module Numeric.Transform.Fourier.R2DIT (fft_r2dit) where -import Data.List import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/R4DIF.hs b/Numeric/Transform/Fourier/R4DIF.hs --- a/Numeric/Transform/Fourier/R4DIF.hs +++ b/Numeric/Transform/Fourier/R4DIF.hs @@ -15,7 +15,6 @@ module Numeric.Transform.Fourier.R4DIF (fft_r4dif) where import DSP.Basic (interleave) -import Data.List import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/Rader.hs b/Numeric/Transform/Fourier/Rader.hs --- a/Numeric/Transform/Fourier/Rader.hs +++ b/Numeric/Transform/Fourier/Rader.hs @@ -14,7 +14,6 @@ module Numeric.Transform.Fourier.Rader (fft_rader1, fft_rader2) where -import Data.List import Data.Array import Data.Complex diff --git a/Numeric/Transform/Fourier/SRDIF.hs b/Numeric/Transform/Fourier/SRDIF.hs --- a/Numeric/Transform/Fourier/SRDIF.hs +++ b/Numeric/Transform/Fourier/SRDIF.hs @@ -15,7 +15,6 @@ module Numeric.Transform.Fourier.SRDIF (fft_srdif) where import DSP.Basic (interleave) --- import Data.List import Data.Array import Data.Complex diff --git a/dsp.cabal b/dsp.cabal --- a/dsp.cabal +++ b/dsp.cabal @@ -1,46 +1,51 @@ Name: dsp -Version: 0.2.1 +Version: 0.2.2 License: GPL License-File: COPYING Copyright: Matt Donadio, 2003 Author: Matt Donadio Maintainer: Henning Thielemann Stability: Experimental -Homepage: http://haskelldsp.sourceforge.net/ -Package-URL: http://darcs.haskell.org/dsp/ +Homepage: http://www.haskell.org/haskellwiki/DSP +-- Homepage: http://haskelldsp.sourceforge.net/ Synopsis: Haskell Digital Signal Processing Description: Digital Signal Processing, Fourier Transform, Linear Algebra, Interpolation -Category: Sound +Category: Sound, Math Tested-With: GHC==6.4.1, GHC==6.8.2 -Cabal-Version: >=1.2 +Tested-With: GHC==7.4.2, GHC==7.6.3 +Cabal-Version: >=1.6 Build-Type: Simple Extra-Source-Files: - Numeric/Special/Airy.gc - Numeric/Special/Erf.gc - Numeric/Special/Foo.gc - Numeric/Special/Clausen.gc - Numeric/Special/Bessel.gc - Numeric/Special/Elljac.gc - Numeric/Special/Ellint.gc demo/Article.hs demo/FFTBench.hs demo/FFTTest.hs demo/FreqDemo.hs demo/IIRDemo.hs demo/NoiseDemo.hs - Makefile +Source-Repository this + Tag: 0.2.2 + Type: darcs + Location: http://code.haskell.org/~thielema/dsp/ + +Source-Repository head + Type: darcs + Location: http://code.haskell.org/~thielema/dsp/ + Flag splitBase description: Choose the new smaller, split-up base package. Library - Build-Depends: mtl If flag(splitBase) - Build-Depends: base >= 2, array, random + Build-Depends: + array >=0.1 && <0.5, + random >=1.0 && <1.1, + base >= 2 && <5 Else - Build-Depends: base >= 1.0 && < 2 + Build-Depends: + base >=1.0 && <2 GHC-Options: -Wall Exposed-modules: @@ -81,6 +86,7 @@ DSP.Source.Basic DSP.Source.Oscillator DSP.Unwrap + DSP.Window Matrix.Cholesky Matrix.LU Matrix.Levinson @@ -122,6 +128,8 @@ Polynomial.Maclaurin Polynomial.Roots DSP.Filter.IIR.Cookbook + Other-Modules: + Numeric.Transform.Fourier.Eigensystem -- Executable: -- Article.hs