diff --git a/DSP/Basic.hs b/DSP/Basic.hs --- a/DSP/Basic.hs +++ b/DSP/Basic.hs @@ -14,58 +14,100 @@ module DSP.Basic where -import Data.Array +import DSP.Source.Basic (zeros) -import DSP.Source.Basic +import Data.Array (Array, Ix, listArray, elems) -- * Functions --- | @z@ is the unit delay function, eg, +-- | 'delay' is the unit delay function, eg, -- --- @z [ 1, 2, 3 ] == [ 0, 1, 2, 3 ]@ +-- @delay1 [ 1, 2, 3 ] == [ 0, 1, 2, 3 ]@ -z :: (Num a) => [a] -> [a] -z a = 0 : a +delay1 :: (Num a) => [a] -> [a] +delay1 a = 0 : a --- | zn is the n sample delay function, eg, --- --- @zn 3 [ 1, 2, 3 ] == [ 0, 0, 0, 1, 2, 3 ]@ +-- | 'delay' is the n sample delay function, eg, +-- +-- @delay 3 [ 1, 2, 3 ] == [ 0, 0, 0, 1, 2, 3 ]@ -zn :: (Num a) => Int -> [a] -> [a] -zn 0 a = a -zn n a = 0 : zn (n - 1) a +delay :: (Num a) => Int -> [a] -> [a] +delay n a = replicate n 0 ++ a -- | @downsample@ throws away every n'th sample, eg, -- -- @downsample 2 [ 1, 2, 3, 4, 5, 6 ] == [ 1, 3, 5 ]@ -downsample :: (Num a) => Int -> [a] -> [a] -downsample n [] = [] -downsample n (x:xs) = x : downsample n (drop (n - 1) xs) +downsample :: Int -> [a] -> [a] +downsample n = + map head . takeWhile (not . null) . iterate (drop n) +downsampleRec :: Int -> [a] -> [a] +downsampleRec _ [] = [] +downsampleRec n (x:xs) = x : downsample n (drop (n - 1) xs) + -- | @upsample@ inserts n-1 zeros between each sample, eg, --- +-- -- @upsample 2 [ 1, 2, 3 ] == [ 1, 0, 2, 0, 3, 0 ]@ upsample :: (Num a) => Int -> [a] -> [a] -upsample _ [] = [] -upsample n (x:xs) = x : zero n n xs - where zero n 1 xs = upsample n xs - zero n i xs = 0 : zero n (i-1) xs +upsample n = concatMap (: replicate (n-1) 0) +upsampleRec :: (Num a) => Int -> [a] -> [a] +upsampleRec _ [] = [] +upsampleRec n (x:xs) = x : zero n xs + where zero 1 ys = upsample n ys + zero i ys = 0 : zero (i-1) ys + -- | @upsampleAndHold@ replicates each sample n times, eg, -- -- @upsampleAndHold 3 [ 1, 2, 3 ] == [ 1, 1, 1, 2, 2, 2, 3, 3, 3 ]@ -upsampleAndHold :: (Num a) => Int -> [a] -> [a] -upsampleAndHold n xs = hold' n n xs - where hold' _ _ [] = [] - hold' n 1 (x:xs) = x : hold' n n xs - hold' n i (x:xs) = x : hold' n (i-1) (x:xs) +upsampleAndHold :: Int -> [a] -> [a] +upsampleAndHold n = concatMap (replicate n) + +-- | merges elements from two lists into one list in an alternating way +-- +-- @interleave [0,1,2,3] [10,11,12,13] == [0,10,1,11,2,12,3,13]@ + +interleave :: [a] -> [a] -> [a] +interleave (e:es) (o:os) = e : o : interleave es os +interleave _ _ = [] + +-- | split a list into two lists in an alternating way +-- +-- @uninterleave [1,2,3,4,5,6] == ([1,3,5],[2,4,6])@ +-- +-- It's a special case of 'Numeric.Random.Spectrum.Pink.split'. + +uninterleave :: [a] -> ([a],[a]) +uninterleave = foldr (\x ~(xs,ys) -> (x:ys,xs)) ([],[]) + + -- | pad a sequence with zeros to length n -- -- @pad [ 1, 2, 3 ] 6 == [ 1, 2, 3, 0, 0, 0 ]@ pad :: (Ix a, Integral a, Num b) => Array a b -> a -> Array a b pad x n = listArray (0,n-1) $ elems x ++ zeros + + +-- | generates a 'Just' if the given condition holds + +toMaybe :: Bool -> a -> Maybe a +toMaybe False _ = Nothing +toMaybe True x = Just x + +-- | Computes the square of the Euclidean norm of a 2D point + +norm2sqr :: Num a => (a,a) -> a +norm2sqr (x,y) = x^!2 + y^!2 + +-- | Power with fixed exponent type. +-- This eliminates warnings about using default types. + +infixr 8 ^! + +(^!) :: Num a => a -> Int -> a +(^!) x n = x^n diff --git a/DSP/Convolution.hs b/DSP/Convolution.hs --- a/DSP/Convolution.hs +++ b/DSP/Convolution.hs @@ -12,7 +12,7 @@ -- ----------------------------------------------------------------------------- -module DSP.Convolution (conv) where +module DSP.Convolution (conv, test) where import Data.Array @@ -21,15 +21,21 @@ -- | @conv@ convolves two finite sequences conv :: (Ix a, Integral a, Num b) => Array a b -> Array a b -> Array a b -conv h1 h2 = h3 - where m1 = snd $ bounds h1 - m2 = snd $ bounds h2 +conv x1 x2 = x3 + where m1 = snd $ bounds x1 + m2 = snd $ bounds x2 m3 = m1 + m2 - h3 = listArray (0,m3) [ sum [ h1!k * h2!(n-k) | k <- [max 0 (n-m2)..min n m1] ] | n <- [0..m3] ] + x3 = listArray (0,m3) [ + sum [ x1!k * x2!(n-k) | k <- [max 0 (n-m2)..min n m1] ] + | n <- [0..m3] ] -- Test vectors. Linear convolution is also equivalent to polynomial -- multiplication. +h1, h2, h3 :: Array Int Integer h1 = listArray (0,3) [ 1, 2, 3, 4 ] h2 = listArray (0,4) [ 1, 2, 3, 4, 5 ] h3 = listArray (0,7) [ 1, 4, 10, 20, 30, 34, 31, 20 ] + +test :: Bool +test = conv h1 h2 == h3 diff --git a/DSP/Correlation.hs b/DSP/Correlation.hs --- a/DSP/Correlation.hs +++ b/DSP/Correlation.hs @@ -10,13 +10,13 @@ -- -- This module contains routines to perform cross- and auto-correlation. -- These formulas can be found in most DSP textbooks. --- +-- -- In the following routines, x and y are assumed to be of the same -- length. -- ----------------------------------------------------------------------------- -module DSP.Correlation (rxy, rxy_b, rxy_u, rxx, rxx_b, rxx_u) where +module DSP.Correlation (rxy, rxy_b, rxy_u, rxx, rxx_b, rxx_u, test) where import Data.Array import Data.Complex @@ -33,8 +33,10 @@ -> a -- ^ k -> Complex b -- ^ R_xy[k] -rxy x y k | k >= 0 = sum [ x!(i+k) * (conjugate (y!i)) | i <- [0..(n-1-k)] ] - | k < 0 = conjugate (rxy y x (-k)) +rxy x y k = + if k >= 0 + then sum [ x!(i+k) * conjugate (y!i) | i <- [0..(n-1-k)] ] + else conjugate (rxy y x (-k)) where n = snd (bounds x) + 1 -- | biased cross-correllation @@ -44,7 +46,7 @@ -> a -- ^ k -> Complex b -- ^ R_xy[k] \/ N -rxy_b x y k = (rxy x y k) / (fromIntegral n) +rxy_b x y k = rxy x y k / fromIntegral n where n = snd (bounds x) + 1 -- | unbiased cross-correllation @@ -54,7 +56,7 @@ -> a -- ^ k -> Complex b -- ^ R_xy[k] \/ (N-k) -rxy_u x y k = (rxy x y k) / (fromIntegral (n-(abs k))) +rxy_u x y k = rxy x y k / fromIntegral (n - abs k) where n = snd (bounds x) + 1 -- autocorrellation @@ -89,18 +91,23 @@ -- test routines ---------------------------------------------------------------------------- -x = array (0,4) [ (0, 1 :+ 0), - (1, 0 :+ 1), - (2, (-1) :+ 0), - (3, 0 :+ (-1)), - (4, 1 :+ 0) ] +xt, yt :: Array Int (Complex Double) +xt = array (0,4) + [ (0, 1 :+ 0), + (1, 0 :+ 1), + (2, (-1) :+ 0), + (3, 0 :+ (-1)), + (4, 1 :+ 0) ] -y = array (0,4) [ (0, 1 :+ 0), - (1, (-1) :+ 0), - (2, 1 :+ 0), - (3, (-1) :+ 0), - (4, 1 :+ 0) ] +yt = array (0,4) + [ (0, 1 :+ 0), + (1, (-1) :+ 0), + (2, 1 :+ 0), + (3, (-1) :+ 0), + (4, 1 :+ 0) ] -r = map (rxy_b x y) [ 0, 1, 2 ] +rt :: [Complex Double] +rt = map (rxy_b xt yt) [ 0, 1, 2 ] -verify = r == [ (0.2 :+ 0.0), (0.0 :+ 0.0), (0.0 :+ 0.2) ] +test :: Bool +test = rt == [ (0.2 :+ 0.0), (0.0 :+ 0.0), (0.0 :+ 0.2) ] diff --git a/DSP/Covariance.hs b/DSP/Covariance.hs --- a/DSP/Covariance.hs +++ b/DSP/Covariance.hs @@ -10,7 +10,7 @@ -- -- This module contains routines to perform cross- and auto-covariance -- These formulas can be found in most DSP textbooks. --- +-- -- In the following routines, x and y are assumed to be of the same -- length. -- @@ -34,7 +34,7 @@ -- -- We define covariance in terms of correlation. -- --- Cxy(X,Y) = E[(X - E[X])(Y - E[Y])] +-- Cxy(X,Y) = E[(X - E[X])(Y - E[Y])] -- = E[XY] - E[X]E[Y] -- = Rxy(X,Y) - E[X]E[Y] @@ -45,15 +45,16 @@ -> a -- ^ k -> Complex b -- ^ C_xy[k] -cxy x y k | k >= 0 = rxy x y k - xm * ym - | k < 0 = conjugate (cxy y x (-k)) +cxy x y k = + if k >= 0 + then rxy x y k - xm * ym + else conjugate (cxy y x (-k)) where xm = mean (elems x) ym = mean (map conjugate (elems y)) - n = snd (bounds x) + 1 -- | raw auto-covariance -- --- Cxx(X,X) = E[(X - E[X])(X - E[X])] +-- Cxx(X,X) = E[(X - E[X])(X - E[X])] -- = E[XX] - E[X]E[X] -- = Rxy(X,X) - E[X]^2 @@ -66,10 +67,10 @@ -> a -- ^ k -> Complex b -- ^ C_xx[k] -cxx x k | k >= 0 = rxx x k - xm^2 - | k < 0 = conjugate (cxx x (-k)) - where xm = mean (elems x) - n = snd (bounds x) + 1 +cxx x k = + if k >= 0 + then rxx x k - mean (elems x) ^ (2::Int) + else conjugate (cxx x (-k)) -- Define the biased and unbiased versions in terms of the above. @@ -80,7 +81,7 @@ -> a -- ^ k -> Complex b -- ^ C_xy[k] \/ N -cxy_b x y k = (cxy x y k) / (fromIntegral n) +cxy_b x y k = cxy x y k / fromIntegral n where n = snd (bounds x) + 1 -- | unbiased cross-covariance @@ -90,7 +91,7 @@ -> a -- ^ k -> Complex b -- ^ C_xy[k] \/ (N-k) -cxy_u x y k = (cxy x y k) / (fromIntegral (n-(abs k))) +cxy_u x y k = cxy x y k / fromIntegral (n - abs k) where n = snd (bounds x) + 1 -- | biased auto-covariance @@ -99,7 +100,7 @@ -> a -- ^ k -> Complex b -- ^ C_xx[k] \/ N -cxx_b x k = (cxx x k) / (fromIntegral n) +cxx_b x k = cxx x k / fromIntegral n where n = snd (bounds x) + 1 -- | unbiased auto-covariance @@ -108,5 +109,5 @@ -> a -- ^ k -> Complex b -- ^ C_xx[k] \/ (N-k) -cxx_u x k = (cxx x k) / (fromIntegral (n-(abs k))) +cxx_u x k = cxx x k / fromIntegral (n - abs k) where n = snd (bounds x) + 1 diff --git a/DSP/Estimation/Frequency/FCI.hs b/DSP/Estimation/Frequency/FCI.hs --- a/DSP/Estimation/Frequency/FCI.hs +++ b/DSP/Estimation/Frequency/FCI.hs @@ -17,10 +17,12 @@ module DSP.Estimation.Frequency.FCI (quinn1, quinn2, quinn3, jacobsen, macleod3, macleod5, rv) where +import DSP.Basic((^!)) import Data.Array import Data.Complex -log10 x = log x / log 10 +log10 :: Floating a => a -> a +log10 = logBase 10 -- | Quinn's First Estimator (FCI1) @@ -44,12 +46,12 @@ -> b -- ^ w quinn2 x k = 2 * pi * ((fromIntegral k) + d) / (fromIntegral n) - where d = (dp + dm) / 2 + tau(dp^2) - tau(dm^2) + where d = (dp + dm) / 2 + tau(dp^!2) - tau(dm^!2) dp = -ap / (1 - ap) dm = am / (1 - am) ap = magnitude (x!(k+1)) / magnitude (x!k) am = magnitude (x!(k-1)) / magnitude (x!k) - tau x = 0.25 * log10(3*x^2 + 6 * x + 1) - (sqrt 6) / 24 * log10 ((x + 1 - sqrt (2/3)) / (x + 1 + sqrt (2/3))) + tau y = 0.25 * log10(3*y^!2 + 6 * y + 1) - (sqrt 6) / 24 * log10 ((y + 1 - sqrt (2/3)) / (y + 1 + sqrt (2/3))) n = snd (bounds x) + 1 -- | Quinn's Third Estimator (FCI3) @@ -59,7 +61,7 @@ -> b -- ^ w quinn3 x k = 2 * pi * ((fromIntegral k) + d) / (fromIntegral n) - where d = (dm + dp) / 2 + (dp - dm) * (3*dt^3 + 2*dt) / (3*dt^4+6*dt^2+1) + where d = (dm + dp) / 2 + (dp - dm) * (3*dt^!3 + 2*dt) / (3*dt^!4+6*dt^!2+1) dt | dm > 0 && dp > 0 = dp | otherwise = dm dp = -ap / (1 - ap) @@ -88,7 +90,7 @@ where rm1 = realPart (x!(k-1) * conjugate (x!k)) r = realPart (x!k * conjugate (x!k)) rp1 = realPart (x!(k+1) * conjugate (x!k)) - d = (sqrt (1 + 8 * g^2) - 1) / 4 / g + d = (sqrt (1 + 8 * g^!2) - 1) / 4 / g g = (rm1 - rp1) / (2 * r + rm1 + rp1) n = snd (bounds x) + 1 @@ -116,6 +118,6 @@ rv x k = 2 * pi * ((fromIntegral k) + d) / (fromIntegral n) where d = fromIntegral at * magnitude (x!(k+at) / x!k) / (1 + magnitude (x!(k+at) / x!k)) - at | (magnitude (x!(k+1)))^2 > (magnitude (x!(k-1)))^2 = 1 + at | (magnitude (x!(k+1)))^!2 > (magnitude (x!(k-1)))^!2 = 1 | otherwise = -1 n = snd (bounds x) + 1 diff --git a/DSP/Estimation/Frequency/PerMax.hs b/DSP/Estimation/Frequency/PerMax.hs --- a/DSP/Estimation/Frequency/PerMax.hs +++ b/DSP/Estimation/Frequency/PerMax.hs @@ -11,23 +11,23 @@ -- This module implements an algorithm to maximize the peak value of a -- DFT\/FFT. It is based off an aticle by Mark Sullivan from Personal -- Engineering Magazine. --- +-- -- Maximizes --- +-- -- @S(w) = 1\/N * sum(k=0,N-1) |x[k] * e^(-jwk)|^2@ --- +-- -- which is equivalent to solving --- +-- -- @S'(w) = Im{X(w) * ~Y(w)} = 0@ --- +-- -- where --- +-- -- @X(w) = sum(k=0,N-1) (x[k] * e^(-jwk))@ -- @Y(w) = X'(w) = sum(k=0,N-1) (k * x[k] * e^(-jwk))@ --- +-- -- This algorithm used the bisection method for finding the zero of a -- function. The search area is +- half a bin width. --- +-- -- Regula falsi requires an additional (x,f(x)) pair which is expensive -- in this case. Newton's method could be used but requires S''(w), -- which takes twice as long to caculate as S'(w). Brent's method may be @@ -48,15 +48,16 @@ -- TODO: the twiddle factor in calc_x,calc_y can be shared -sign x | x < 0 = -1 - | x == 0 = 0 - | x > 0 = 1 -- calc_x x w = sum [ x!k * cis (-w * fromIntegral k) | k <- [0..(n-1)] ] -- where n = snd (bounds x) + 1 +calc_x :: (RealFloat a, Ix i) => + Array i (Complex a) -> a -> Complex a calc_x x w = sum $ zipWith (*) (elems x) (iterate (cis (-w) *) 1) +calc_y :: (RealFloat b, Ix i, Integral i) => + Array i (Complex b) -> b -> Complex b calc_y x w = sum [ fromIntegral k * x!k * cis (-w * fromIntegral k) | k <- [0..(n-1)] ] where n = snd (bounds x) + 1 @@ -69,16 +70,18 @@ permax x k = permax' x (w-d) (w+d) where w = 2 * pi * fromIntegral k / fromIntegral n d = 1 / fromIntegral (2*n) -- half a bin width - n = snd (bounds x) + 1 + n = snd (bounds x) + 1 +permax' :: (RealFloat b, Ix i, Integral i) => + Array i (Complex b) -> b -> b -> b permax' x w0 w1 | w1-w0 < eps = wmid - | otherwise = if sign t0 == sign tm - then permax' x wmid w1 - else permax' x w0 wmid + | otherwise = if signum t0 == signum tm + then permax' x wmid w1 + else permax' x w0 wmid where t0 = imagPart ((calc_x x w0) * (conjugate (calc_y x w0))) - tm = imagPart ((calc_x x wmid) * (conjugate (calc_y x wmid))) - t1 = imagPart ((calc_x x w1) * (conjugate (calc_y x w1))) + tm = imagPart ((calc_x x wmid) * (conjugate (calc_y x wmid))) +-- t1 = imagPart ((calc_x x w1) * (conjugate (calc_y x w1))) wmid = (w0 + w1) / 2 -- bisection method -- wmid = w1 - t1 * (w1 - w0) / (t1 - t0) -- regula falsi eps = 1.0e-6 - n = snd (bounds x) + 1 +-- n = snd (bounds x) + 1 diff --git a/DSP/Estimation/Frequency/Pisarenko.hs b/DSP/Estimation/Frequency/Pisarenko.hs --- a/DSP/Estimation/Frequency/Pisarenko.hs +++ b/DSP/Estimation/Frequency/Pisarenko.hs @@ -19,8 +19,14 @@ module DSP.Estimation.Frequency.Pisarenko (pisarenko) where +import DSP.Basic((^!)) import Data.Array +rss :: (Ix a, Integral a, Num b) => + Array a b + -> a + -> b + rss x k = sum [ x!(i+k) * x!i | i <- [0..(n-1-k)] ] where n = snd (bounds x) + 1 @@ -30,7 +36,7 @@ -> b -- ^ w pisarenko x = acos (alpha / 2) - where alpha = (rss2 + sqrt (rss2^2 + 8*rss1^2)) / (rss1 + eps) / 2 + where alpha = (rss2 + sqrt (rss2^!2 + 8*rss1^!2)) / (2*(rss1 + eps)) rss1 = rss x 1 rss2 = rss x 2 eps = 1.0e-15 diff --git a/DSP/Estimation/Frequency/QuinnFernandes.hs b/DSP/Estimation/Frequency/QuinnFernandes.hs --- a/DSP/Estimation/Frequency/QuinnFernandes.hs +++ b/DSP/Estimation/Frequency/QuinnFernandes.hs @@ -25,9 +25,13 @@ qf y w = qf' y (2 * cos w) +qf' :: (Ix a, Integral a, RealFloat b) => + Array a b + -> b + -> b qf' y a | abs (a-b) < eps = acos(0.5 * b) | otherwise = qf' y b where z = array (-2,n-1) ([ (-2, 0), (-1, 0) ] ++ [ (i, y!i + a * z!(i-1) - z!(i-2)) | i <- [0..(n-1)] ]) - b = sum [ (z!i + z!(i-2)) * z!(i-1) | i <- [0..(n-1)] ] / sum [ (z!(i-1))^2 | i <- [0..(n-1)] ] + b = sum [ (z!i + z!(i-2)) * z!(i-1) | i <- [0..(n-1)] ] / sum [ (z!(i-1))^(2::Int) | i <- [0..(n-1)] ] eps = 1.0e-6 n = snd (bounds y) + 1 diff --git a/DSP/Estimation/Frequency/WLP.hs b/DSP/Estimation/Frequency/WLP.hs --- a/DSP/Estimation/Frequency/WLP.hs +++ b/DSP/Estimation/Frequency/WLP.hs @@ -13,10 +13,9 @@ -- ----------------------------------------------------------------------------- --- Boy, fromIntegral makes these look really messy. - -module DSP.Estimation.Frequency.WLP where +module DSP.Estimation.Frequency.WLP (wlp, lrp, kay, lw, ckq,) where +import DSP.Basic((^!)) import Data.Array import Data.Complex @@ -34,24 +33,27 @@ lrp :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z -> b -- ^ w -lrp z = wlp (array (1,n-1) [ (t, 1 / fromIntegral (n-1)) | t <- [1..(n-1)] ]) z - where n = snd (bounds z) + 1 +lrp = processArray (\n _ _ -> recip (n-1)) -- | WLP using kay's window kay :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z -> b -- ^ w -kay z = wlp (array (1,n-1) [ (t, fromIntegral (6*t*(n-t)) / fromIntegral (n*(n^2-1))) | t <- [1..(n-1)] ]) z - where n = snd (bounds z) + 1 +kay = processArray (\n _ t -> kayWin n t) +kayWin :: Fractional b => b -> b -> b +kayWin n t = 6*t*(n-t) / (n*(n^!2-1)) + -- | WLP using Lovell and Williamson's window lw :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ z -> b -- ^ w -lw z = wlp (array (1,n-1) [ (t, fromIntegral (6*t*(n-t)) / (fromIntegral (n*(n^2-1)) * magnitude (z!t) * magnitude (conjugate (z!(t-1))))) | t <- [1..(n-1)] ]) z - where n = snd (bounds z) + 1 +lw z = + processArray + (\ n ti t -> kayWin n t / + (magnitude (z!ti) * magnitude (z!(ti-1)))) z -- | WLP using Clarkson, Kootsookos, and Quinn's window @@ -60,8 +62,23 @@ -> b -- ^ sigma -> b -- ^ w -ckq z rho sig = wlp (array (1,n-1) [ (t, num t / den) | t <- [1..(n-1)] ]) z - where num t = sinh (fromIntegral n * th) - sinh (fromIntegral t * th) - sinh (fromIntegral (n-t) * th) - den = fromIntegral (n-1) * sinh (fromIntegral n * th) - 2 * sinh (0.5 * fromIntegral n * th) * sinh (0.5 * fromIntegral (n-1) * th) / sinh (0.5 * th) - th = log (1 + sig^2 / rho^2 + sqrt (sig^4 / rho^4 + sig^2 / rho^2)) - n = snd (bounds z) + 1 +ckq z rho sig = + processArray + (\n -> + let num t = sinh (n * th) - sinh (t * th) - sinh ((n-t) * th) + den = (n-1) * sinh (n * th) - + 2 * sinh (0.5 * n * th) + * sinh (0.5 * (n-1) * th) / sinh (0.5 * th) + sigRho2 = (sig / rho) ^! 2 + th = log (1 + sig^!2 / rho^!2 + sqrt ((sigRho2+1) * sigRho2)) + in \_ t -> num t / den) z + +processArray :: + (Integral a, Ix a, RealFloat b) => + (b -> a -> b -> b) -> Array a (Complex b) -> b +processArray f z = + let n = snd (bounds z) + 1 + g = f (fromIntegral n) + bnd = (1,n-1) + in wlp (listArray bnd + (map (\t -> g t (fromIntegral t)) (range bnd))) z diff --git a/DSP/Estimation/Spectral/AR.hs b/DSP/Estimation/Spectral/AR.hs --- a/DSP/Estimation/Spectral/AR.hs +++ b/DSP/Estimation/Spectral/AR.hs @@ -18,13 +18,14 @@ module DSP.Estimation.Spectral.AR where -import Data.Array -import Data.Complex - import DSP.Correlation import Matrix.Levinson import Matrix.Cholesky +import DSP.Basic((^!)) +import Data.Array +import Data.Complex + -- * Functions ------------------------------------------------------------------------------- @@ -106,8 +107,8 @@ where a = array ((1,1),(p,p)) [ ((k,i), ak k i) | k <- [1..p], i <- [1..k] ] ak k i | i==k = kk!k | otherwise = a!(k-1,i) + kk!k * (conjugate (a!(k-1,k-i))) - kk = array (1,p) [ (k, -2 * sum [ ef!((k-1),i) * (conjugate (eb!(k-1,i-1))) | i <- [k..(n-1)] ] / sum [ (abs (ef!(k-1,i)))^2 + (abs (eb!(k-1,i-1)))^2 | i <- [k..(n-1)] ]) | k <- [1..p] ] - rho = array (0,p) ((0, rxx_b x 0) : [ (k, (1 - ((abs (kk!k))^2)) * rho!(k-1)) | k <- [1..p] ]) + kk = array (1,p) [ (k, -2 * sum [ ef!((k-1),i) * (conjugate (eb!(k-1,i-1))) | i <- [k..(n-1)] ] / sum [ abs (ef!(k-1,i)) ^! 2 + abs (eb!(k-1,i-1)) ^! 2 | i <- [k..(n-1)] ]) | k <- [1..p] ] + rho = array (0,p) ((0, rxx_b x 0) : [ (k, (1 - abs (kk!k) ^! 2) * rho!(k-1)) | k <- [1..p] ]) ef = array ((0,1),(p,n-1)) [ ((k,i), efki k i) | k <- [0..p], i <- [(k+1)..(n-1)] ] eb = array ((0,0),(p,n-2)) [ ((k,i), ebki k i) | k <- [0..p], i <- [k..(n-2)] ] efki 0 i = x!i diff --git a/DSP/Estimation/Spectral/ARMA.hs b/DSP/Estimation/Spectral/ARMA.hs --- a/DSP/Estimation/Spectral/ARMA.hs +++ b/DSP/Estimation/Spectral/ARMA.hs @@ -24,14 +24,18 @@ import Data.Complex import DSP.Correlation -import DSP.Estimation.Spectral.MA +-- import DSP.Estimation.Spectral.MA -import Matrix.LU +-- import Matrix.LU + -- * Functions --- THIS DOES NOT WORK +-- | THIS DOES NOT WORK +arma_mywe :: + (RealFloat b, Integral i, Ix i) => + Array i (Complex b) -> i -> i -> Array i (Complex b) arma_mywe x p q = a' where r = array (q-2*p+1,q+p) [ (k, rxx_u x k) | k <- [(q-2*p+1)..(q+p)] ] a' = array (1,p) [ (k, a!(p,k)) | k <- [1..p] ] diff --git a/DSP/Estimation/Spectral/KayData.hs b/DSP/Estimation/Spectral/KayData.hs --- a/DSP/Estimation/Spectral/KayData.hs +++ b/DSP/Estimation/Spectral/KayData.hs @@ -14,76 +14,78 @@ module DSP.Estimation.Spectral.KayData (xc,xr) where -import Data.Array -import Data.Complex +import Data.Array (Array, array) +import Data.Complex (Complex((:+))) -- | Complex test data xc :: Array Int (Complex Double) -xc = array (0,31) [ (0, (( 6.3307) :+ (-0.174915))), - (1, ((-1.33539) :+ (-0.03044))), - (2, (( 3.61896) :+ (-0.260459))), - (3, (( 1.87513) :+ (-0.323974))), - (4, ((-1.08561) :+ (-0.136055))), - (5, (( 3.99114) :+ (-0.101864))), - (6, ((-4.10184) :+ ( 0.130571))), - (7, (( 1.55399) :+ ( 0.0977916))), - (8, ((-2.1258) :+ (-0.306485))), - (9, ((-3.27873) :+ (-0.0544436))), - (10, (( 0.241218) :+ ( 0.0962379))), - (11, ((-5.74708 ) :+ ( 0.0186908))), - (12, ((-0.0165977) :+ ( 0.237493))), - (13, ((-3.28921) :+ (-0.188478))), - (14, ((-1.31227) :+ (-0.120636))), - (15, (( 0.745251) :+ (-0.0679575))), - (16, ((-1.77199) :+ (-0.416229))), - (17, (( 2.56419) :+ (-0.270373))), - (18, (( 0.21325) :+ (-0.232544))), - (19, (( 2.23409) :+ ( 0.236383))), - (20, (( 2.2949) :+ ( 0.173061))), - (21, (( 1.09186) :+ ( 0.140938))), - (22, (( 2.29353) :+ ( 0.442044))), - (23, (( 0.695823) :+ ( 0.509325))), - (24, (( 0.759858) :+ ( 0.417967))), - (25, ((-0.354267) :+ ( 0.506891))), - (26, ((-0.594517) :+ ( 0.39708))), - (27, ((-1.88618) :+ ( 0.649179))), - (28, ((-1.39041) :+ ( 0.867086))), - (29, ((-3.06381) :+ ( 0.422965))), - (30, ((-2.0433) :+ ( 0.0825514))), - (31, ((-2.1628) :+ (-0.0933218))) ] +xc = array (0,31) [ + (0, (( 6.3307) :+ (-0.174915))), + (1, ((-1.33539) :+ (-0.03044))), + (2, (( 3.61896) :+ (-0.260459))), + (3, (( 1.87513) :+ (-0.323974))), + (4, ((-1.08561) :+ (-0.136055))), + (5, (( 3.99114) :+ (-0.101864))), + (6, ((-4.10184) :+ ( 0.130571))), + (7, (( 1.55399) :+ ( 0.0977916))), + (8, ((-2.1258) :+ (-0.306485))), + (9, ((-3.27873) :+ (-0.0544436))), + (10, (( 0.241218) :+ ( 0.0962379))), + (11, ((-5.74708 ) :+ ( 0.0186908))), + (12, ((-0.0165977) :+ ( 0.237493))), + (13, ((-3.28921) :+ (-0.188478))), + (14, ((-1.31227) :+ (-0.120636))), + (15, (( 0.745251) :+ (-0.0679575))), + (16, ((-1.77199) :+ (-0.416229))), + (17, (( 2.56419) :+ (-0.270373))), + (18, (( 0.21325) :+ (-0.232544))), + (19, (( 2.23409) :+ ( 0.236383))), + (20, (( 2.2949) :+ ( 0.173061))), + (21, (( 1.09186) :+ ( 0.140938))), + (22, (( 2.29353) :+ ( 0.442044))), + (23, (( 0.695823) :+ ( 0.509325))), + (24, (( 0.759858) :+ ( 0.417967))), + (25, ((-0.354267) :+ ( 0.506891))), + (26, ((-0.594517) :+ ( 0.39708))), + (27, ((-1.88618) :+ ( 0.649179))), + (28, ((-1.39041) :+ ( 0.867086))), + (29, ((-3.06381) :+ ( 0.422965))), + (30, ((-2.0433) :+ ( 0.0825514))), + (31, ((-2.1628) :+ (-0.0933218))) ] -- | Real test data xr :: Array Int Double -xr = array (0,31) [ (0, 6.46768), - (1, -1.28024), - (2, 3.74788), - (3, 1.96092), - (4, -0.768349), - (5, 4.14569), - (6, -4.05277), - (7, 1.65836), - (8, -2.06405), - (9, -3.33397), - (10, 0.085145), - (11, -6.06562), - (12, -0.411658), - (13, -3.61831), - (14, -1.53352), - (15, 0.481522), - (16, -1.93653), - (17, 2.35532), - (18, 0.145624), - (19, 2.21991), - (20, 2.25884), - (21, 1.07373), - (22, 2.26531), - (23, 0.685007), - (24, 0.762859), - (25, -0.501008), - (26, -0.640518), - (27, -1.99263), - (28, -1.60416), - (29, -3.22751), - (30, -2.21946), - (31, -2.42246) ] +xr = array (0,31) [ + (0, 6.46768), + (1, -1.28024), + (2, 3.74788), + (3, 1.96092), + (4, -0.768349), + (5, 4.14569), + (6, -4.05277), + (7, 1.65836), + (8, -2.06405), + (9, -3.33397), + (10, 0.085145), + (11, -6.06562), + (12, -0.411658), + (13, -3.61831), + (14, -1.53352), + (15, 0.481522), + (16, -1.93653), + (17, 2.35532), + (18, 0.145624), + (19, 2.21991), + (20, 2.25884), + (21, 1.07373), + (22, 2.26531), + (23, 0.685007), + (24, 0.762859), + (25, -0.501008), + (26, -0.640518), + (27, -1.99263), + (28, -1.60416), + (29, -3.22751), + (30, -2.21946), + (31, -2.42246) ] diff --git a/DSP/Estimation/Spectral/MA.hs b/DSP/Estimation/Spectral/MA.hs --- a/DSP/Estimation/Spectral/MA.hs +++ b/DSP/Estimation/Spectral/MA.hs @@ -44,4 +44,3 @@ where (b,_) = ar_yw a' q a' = array (0,l) ((0,1) : [ (i, a''!i) | i <- [1..l] ]) (a'', sig2) = ar_yw x l - n = snd (bounds x) + 1 diff --git a/DSP/Filter/Analog/Prototype.hs b/DSP/Filter/Analog/Prototype.hs --- a/DSP/Filter/Analog/Prototype.hs +++ b/DSP/Filter/Analog/Prototype.hs @@ -33,9 +33,9 @@ module DSP.Filter.Analog.Prototype where -import Data.Complex +import Data.Complex (Complex((:+)), realPart) -import Polynomial.Basic +import Polynomial.Basic (roots2poly) -- | Generates Butterworth filter prototype @@ -46,7 +46,7 @@ where poles = [ (-u k) :+ (w k) | k <- [0..(n-1)] ] u k = sin (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) w k = cos (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) - num = [ 1 ] + num = [ 1 ] den = map realPart $ roots2poly $ poles -- | Generates Chebyshev filter prototype @@ -62,8 +62,10 @@ num = [ gain ] den = map realPart $ roots2poly $ poles v0 = asinh (1/eps) / fromIntegral n - gain | even n = abs $ head den / sqrt (1 + eps^2) - | odd n = abs $ head den + gain = + if even n + then abs $ head den / sqrt (1 + eps^(2::Int)) + else abs $ head den -- | Generates Inverse Chebyshev filter prototype diff --git a/DSP/Filter/Analog/Response.hs b/DSP/Filter/Analog/Response.hs --- a/DSP/Filter/Analog/Response.hs +++ b/DSP/Filter/Analog/Response.hs @@ -16,6 +16,7 @@ module DSP.Filter.Analog.Response where +import DSP.Basic ((^!)) import Polynomial.Basic import Polynomial.Chebyshev @@ -36,8 +37,8 @@ -> Double -- ^ w -> Double -- ^ |H_c(w)|^2 -chebyshev1_H n eps wc w = 1 / (1 + eps^2 * vn(w/wc)^2) - where vn w = polyeval (cheby n) w +chebyshev1_H n eps wc w = 1 / (1 + eps^!2 * vn(w/wc)^!2) + where vn = polyeval (cheby n) -- | Inverse Chebyshev filter response function -- @@ -49,5 +50,5 @@ -> Double -- ^ w -> Double -- ^ |H_c(w)|^2 -chebyshev2_H n eps wc w = 1 / (1 + (eps^2 * vn(wc/w)^2)**(-1)) - where vn w = polyeval (cheby n) w +chebyshev2_H n eps wc w = 1 / (1 + (eps^!2 * vn(wc/w)^!2)**(-1)) + where vn = polyeval (cheby n) diff --git a/DSP/Filter/Analog/Transform.hs b/DSP/Filter/Analog/Transform.hs --- a/DSP/Filter/Analog/Transform.hs +++ b/DSP/Filter/Analog/Transform.hs @@ -9,29 +9,31 @@ -- Portability : portable -- -- Analog prototype filter transforms ---- +--- -- Reference: R&G, pg 258; P&M, pg 698 -- ----------------------------------------------------------------------------- -module DSP.Filter.Analog.Transform (a_lp2lp, a_lp2hp, a_lp2bp, a_lp2bs) where - -import Data.Complex +module DSP.Filter.Analog.Transform ( + a_lp2lp, a_lp2hp, a_lp2bp, a_lp2bs, + substitute, propSubstituteRecip, propSubstituteAlt, + ) where import Polynomial.Basic -- Normalizes a filter +normalize :: ([Double],[Double]) -> ([Double],[Double]) normalize (num,den) = (num',den') where a0 = last den - num' = map (/ a0) num - den' = map (/ a0) den + num' = map (/ a0) num + den' = map (/ a0) den -- | Lowpass to lowpass: @s --> s\/wc@ a_lp2lp :: Double -- ^ wc - -> ([Double],[Double]) -- ^ (b,a) - -> ([Double],[Double]) -- ^ (b',a') + -> ([Double],[Double]) -- ^ (b,a) + -> ([Double],[Double]) -- ^ (b',a') a_lp2lp wu (num,den) = normalize (num',den') where num' = polysubst [ 0, 1/wu ] num @@ -40,46 +42,65 @@ -- | Lowpass to highpass: @s --> wc\/s@ a_lp2hp :: Double -- ^ wc - -> ([Double],[Double]) -- ^ (b,a) - -> ([Double],[Double]) -- ^ (b',a') + -> ([Double],[Double]) -- ^ (b,a) + -> ([Double],[Double]) -- ^ (b',a') a_lp2hp wu (num,den) = normalize (num',den') where nn = length num - nd = length den - n = max nn nd - num' = polysubst [ 0, 1/wu ] $ reverse $ num ++ replicate (n-nn) 0 - den' = polysubst [ 0, 1/wu ] $ reverse $ den ++ replicate (n-nd) 0 + nd = length den + n = max nn nd + num' = polysubst [ 0, 1/wu ] $ reverse $ num ++ replicate (n-nn) 0 + den' = polysubst [ 0, 1/wu ] $ reverse $ den ++ replicate (n-nd) 0 + -- | Lowpass to bandpass: @s --> (s^2 + wl*wu) \/ (s(wu-wl))@ a_lp2bp :: Double -- ^ wl - -> Double -- ^ wu - -> ([Double],[Double]) -- ^ (b,a) - -> ([Double],[Double]) -- ^ (b',a') + -> Double -- ^ wu + -> ([Double],[Double]) -- ^ (b,a) + -> ([Double],[Double]) -- ^ (b',a') -a_lp2bp wl wu (num,den) = normalize (num',den') - where n = max (length num - 1) (length den - 1) - num' = step3 $ step2 n [ 0, wu-wl ] $ step1 0 [ wl*wu, 0, 1 ] $ num - den' = step3 $ step2 n [ 0, wu-wl ] $ step1 0 [ wl*wu, 0, 1 ] $ den - step1 _ _ [] = [] - step1 n w (x:xs) = map (x*) (polypow w n) : step1 (n+1) w xs - step2 _ _ [] = [] - step2 n w (x:xs) = polymult (polypow w n) x : step2 (n-1) w xs - step3 x = foldr polyadd [0] x +a_lp2bp wl wu = substitute ([ wl*wu, 0, 1 ], [ 0, wu-wl ]) + -- | Lowpass to bandstop: @s --> (s(wu-wl)) \/ (s^2 + wl*wu)@ a_lp2bs :: Double -- ^ wl - -> Double -- ^ wu - -> ([Double],[Double]) -- ^ (b,a) - -> ([Double],[Double]) -- ^ (b',a') + -> Double -- ^ wu + -> ([Double],[Double]) -- ^ (b,a) + -> ([Double],[Double]) -- ^ (b',a') -a_lp2bs wl wu (num,den) = normalize (num',den') - where n = max (length num - 1) (length den - 1) - num' = step3 $ step2 n [ wu*wl, 0, 1 ] $ step1 0 [ 0, wu-wl ] $ num - den' = step3 $ step2 n [ wu*wl, 0, 1 ] $ step1 0 [ 0, wu-wl ] $ den - step1 _ _ [] = [] - step1 n w (x:xs) = map (x*) (polypow w n) : step1 (n+1) w xs - step2 _ _ [] = [] - step2 n w (x:xs) = polymult (polypow w n) x : step2 (n-1) w xs - step3 x = foldr polyadd [0] x +a_lp2bs wl wu = substitute ([ 0, wu-wl ], [ wu*wl, 0, 1 ]) + + + +substitute :: + ([Double],[Double]) -> ([Double],[Double]) -> ([Double],[Double]) +substitute (nsub,dsub) (num,den) = normalize (num',den') + where num' = polyPolySubst nsub $ weightedPowers $ num + den' = polyPolySubst nsub $ weightedPowers $ den + weightedPowers = flip (zipWith polyscale) dsubPowers + dsubPowers = reverse $ take m $ iterate (polymult dsub) [1] + m = max (length num) (length den) + +substituteAlt :: + ([Double],[Double]) -> ([Double],[Double]) -> ([Double],[Double]) +substituteAlt (nsub,dsub) (num,den) = normalize (num',den') + where m = max (length num - 1) (length den - 1) + num' = step3 $ step2 (0::Int) $ step1 m $ num + den' = step3 $ step2 (0::Int) $ step1 m $ den + step1 _ [] = [] + step1 n (x:xs) = map (x*) (polypow dsub n) : step1 (n-1) xs + step2 _ [] = [] + step2 n (x:xs) = polymult (polypow nsub n) x : step2 (n+1) xs + step3 x = foldr polyadd [0] x + + +propSubstituteRecip :: ([Double],[Double]) -> ([Double],[Double]) -> Bool +propSubstituteRecip (nsub,dsub) (num,den) = + let (x,y) = substitute (nsub,dsub) (num,den) + in (y,x) == substitute (dsub,nsub) (den,num) + + +propSubstituteAlt :: ([Double],[Double]) -> ([Double],[Double]) -> Bool +propSubstituteAlt q p = substitute q p == substituteAlt q p 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 @@ -12,12 +12,12 @@ -- ----------------------------------------------------------------------------- -module DSP.Filter.FIR.FIR (fir) where +module DSP.Filter.FIR.FIR (fir, test) where import Data.Array -- | Implements the following function, which is a FIR filter --- +-- -- @y[n] = sum(k=0,M) h[k]*x[n-k]@ -- -- We implement the fir function with five helper functions, depending on @@ -32,6 +32,7 @@ -> [a] -- ^ x[n] -> [a] -- ^ y[n] +fir _ [] = [] fir h (x:xs) | isFIRType1 h = fir'1 h w xs | isFIRType2 h = fir'2 h w xs | isFIRType3 h = fir'3 h w xs @@ -42,6 +43,8 @@ -- This is for testing the symetric helpers. +fir0 :: Num a => Array Int a -> [a] -> [a] +fir0 _ [] = [] fir0 h (x:xs) = fir'0 h w xs where w = listArray (0,m) $ x : replicate m 0 m = snd $ bounds h @@ -147,36 +150,37 @@ {-# specialize isFIRType3 :: Array Int Double -> Bool #-} isFIRType3 :: Num a => Array Int a -> Bool -isFIRType3 h = even m && h1 == reverse h2 +isFIRType3 h = even m && ha == reverse hb where m = snd $ bounds h h' = elems h - h1 = take n h' - h2 = map negate (drop (n+1) h') + ha = take n h' + hb = map negate (drop (n+1) h') n = m `div` 2 {-# specialize isFIRType4 :: Array Int Float -> Bool #-} {-# specialize isFIRType4 :: Array Int Double -> Bool #-} isFIRType4 :: Num a => Array Int a -> Bool -isFIRType4 h = odd m && h1 == reverse h2 +isFIRType4 h = odd m && ha == reverse hb where m = snd $ bounds h - h1 = elems h - h2 = fmap negate $ h1 + ha = elems h + hb = fmap negate $ ha -- Test routines -- This tests out fir'0 -h :: Array Int Double -h = listArray (0,4) [ 1, 2, 0, -1, 1 ] +ht :: Array Int Double +ht = listArray (0,4) [ 1, 2, 0, -1, 1 ] -x :: [Double] -x = [1, 3, -1, -2, 0, 0, 0, 0 ] +xt :: [Double] +xt = [1, 3, -1, -2, 0, 0, 0, 0 ] -y :: [Double] -y = [1, 5, 5, -5, -6, 4, 1, -2] +yt :: [Double] +yt = [1, 5, 5, -5, -6, 4, 1, -2] -y' = fir h x +yt' :: [Double] +yt' = fir ht xt -- This checks the symetric routines against fir'0 @@ -189,16 +193,19 @@ h4 :: Array Int Double h4 = listArray (0,5) [ 1, 2, 3, -3, -2, -1 ] -y1 = fir0 h1 x -y2 = fir0 h2 x -y3 = fir0 h3 x -y4 = fir0 h4 x +y1, y2, y3, y4 :: [Double] +y1 = fir0 h1 xt +y2 = fir0 h2 xt +y3 = fir0 h3 xt +y4 = fir0 h4 xt -y1' = fir h1 x -y2' = fir h2 x -y3' = fir h3 x -y4' = fir h4 x +y1', y2', y3', y4' :: [Double] +y1' = fir h1 xt +y2' = fir h2 xt +y3' = fir h3 xt +y4' = fir h4 xt -- If everything works, then test == True -test = foldr (&&) True [ y == y', y1 == y1', y2 == y2', y3 == y3', y4 == y4' ] +test :: Bool +test = and [ yt == yt', y1 == y1', y2 == y2', y3 == y3', y4 == y4' ] 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 @@ -14,7 +14,7 @@ ----------------------------------------------------------------------------- -- Reference: --- +-- -- @Book{dsp, -- author = "Alan V. Oppenheim and Ronald W. Schafer", -- title = "Discrete-Time Signal Processing", @@ -34,24 +34,29 @@ -- Set the cutoff frequency to the middle of the transition band. This -- equation isn't numbered. +calc_wc :: Fractional a => a -> a -> a calc_wc wp ws = (wp + ws) / 2 -- Equation 7.90 +calc_dw :: Num a => a -> a -> a calc_dw wp ws = abs (ws - wp) -- Equation 7.91 +calc_A :: (Floating a, Ord a) => a -> a -> a calc_A d1 d2 = -20 * logBase 10 (min d1 d2) -- xEquation 7.92 +calc_beta :: (Ord a, Floating a) => a -> a calc_beta a | a > 50 = 0.1102 * (a - 8.7) | a >= 21 = 0.5842 * ((a-21) ** 0.4) + 0.07886 * (a-21) | otherwise = 0.0 -- Equation 7.93 +calc_M :: (Integral b, RealFrac a) => a -> a -> b calc_M a dw = ceiling ((a - 8) / (2.285 * dw)) -- Procedure on pg 455. We should really check the peak approximation @@ -91,5 +96,7 @@ dw = calc_dw wp ws a = calc_A d1 d2 beta = calc_beta a - m | odd (calc_M a dw) = (calc_M a dw) + 1 - | otherwise = (calc_M a dw) + m = ceilingEven (calc_M a dw) + +ceilingEven :: Integral b => b -> b +ceilingEven x = x + mod (-x) 2 diff --git a/DSP/Filter/FIR/PolyInterp.hs b/DSP/Filter/FIR/PolyInterp.hs --- a/DSP/Filter/FIR/PolyInterp.hs +++ b/DSP/Filter/FIR/PolyInterp.hs @@ -9,7 +9,7 @@ -- Portability : portable -- -- Polynomial interpolators. Taken from: --- +-- -- Olli Niemitalo (ollinie\@freenet.hut.fi), "Polynomial Interpolators for -- High-Quality Resampling of Oversampled Audio" Search for "deip.pdf" with -- Google and you will find it. @@ -57,22 +57,22 @@ -- B-Splines -bspline_1p0o :: (Ord a, Fractional a) => a -> a +bspline_1p0o :: (Ord a, Fractional a) => a -> a bspline_1p0o x | 0 <= x && x < 1 = polyeval [ 1 ] x | otherwise = 0 -bspline_2p1o :: (Ord a, Fractional a) => a -> a +bspline_2p1o :: (Ord a, Fractional a) => a -> a bspline_2p1o x | 0 <= x && x < 1 = polyeval [ 1, -1 ] x | 1 <= x = 0 | otherwise = bspline_2p1o (-x) -bspline_4p3o :: (Ord a, Fractional a) => a -> a +bspline_4p3o :: (Ord a, Fractional a) => a -> a bspline_4p3o x | 0 <= x && x < 1 = polyeval [ 2/3, 0, -1, 1/2 ] x | 1 <= x && x < 2 = polyeval [ 4/3, -2, 1, -1/6 ] x | 2 <= x = 0 | otherwise = bspline_4p3o (-x) -bspline_6p5o :: (Ord a, Fractional a) => a -> a +bspline_6p5o :: (Ord a, Fractional a) => a -> a bspline_6p5o x | 0 <= x && x < 1 = polyeval [ 11/20, 0, -1/2, 0, 1/4, -1/12 ] x | 1 <= x && x < 2 = polyeval [ 17/40, 5/8, -7/4, 5/4, -3/8, 1/24 ] x | 2 <= x && x < 3 = polyeval [ 81/40, -27/8, 9/4, -3/4, 1/8, -1/120 ] x @@ -83,13 +83,13 @@ -- Lagrange polynomials -lagrange_4p3o :: (Ord a, Fractional a) => a -> a +lagrange_4p3o :: (Ord a, Fractional a) => a -> a lagrange_4p3o x | 0 <= x && x < 1 = polyeval [ 1, -1/2, -1, 1/2 ] x | 1 <= x && x < 2 = polyeval [ 1, -11/6, 1, -1/6 ] x | 2 <= x = 0 | otherwise = lagrange_4p3o (-x) -lagrange_6p5o :: (Ord a, Fractional a) => a -> a +lagrange_6p5o :: (Ord a, Fractional a) => a -> a lagrange_6p5o x | 0 <= x && x < 1 = polyeval [ 1, -1/3, -5/4, 5/12, 1/4, -1/12 ] x | 1 <= x && x < 2 = polyeval [ 1, -13/12, -5/8, 25/24, -3/8, 1/24 ] x | 2 <= x && x < 3 = polyeval [ 1, -137/60, 15/8, -17/24, 1/8, -1/120 ] x @@ -100,20 +100,20 @@ -- Hermite (1st-order-osculating) polynomials -hermite_4p3o :: (Ord a, Fractional a) => a -> a +hermite_4p3o :: (Ord a, Fractional a) => a -> a hermite_4p3o x | 0 <= x && x < 1 = polyeval [ 1, 0, -5/2, 3/2 ] x | 1 <= x && x < 2 = polyeval [ 2, -4, 5/2, -1/2 ] x | 2 <= x = 0 | otherwise = hermite_4p3o (-x) -hermite_6p3o :: (Ord a, Fractional a) => a -> a +hermite_6p3o :: (Ord a, Fractional a) => a -> a hermite_6p3o x | 0 <= x && x < 1 = polyeval [ 1, 0, -7/3, 4/3 ] x | 1 <= x && x < 2 = polyeval [ 5/2, -59/12, 3, -7/12 ] x | 2 <= x && x < 3 = polyeval [ -3/2, 7/4, -2/3, 1/12 ] x | 3 <= x = 0 | otherwise = hermite_6p3o (-x) -hermite_6p5o :: (Ord a, Fractional a) => a -> a +hermite_6p5o :: (Ord a, Fractional a) => a -> a hermite_6p5o x | 0 <= x && x < 1 = polyeval [ 1, 0, -25/12, 5/12, 13/12, -5/12 ] x | 1 <= x && x < 2 = polyeval [ 1, 5/12, -35/8, 35/8, -13/8, 5/24 ] x | 2 <= x && x < 3 = polyeval [ 3, -29/4, 155/24, -65/24, 13/24, -1/24 ] x @@ -124,13 +124,13 @@ -- 2nd-order-osculating polynomials -sndosc_4p5o :: (Ord a, Fractional a) => a -> a +sndosc_4p5o :: (Ord a, Fractional a) => a -> a sndosc_4p5o x | 0 <= x && x < 1 = polyeval [ 1, 0, -1, -9/2, 15/2, -3 ] x | 1 <= x && x < 2 = polyeval [ -4, 18, -29, 43/2, -15/2, 1 ] x | 2 <= x = 0 | otherwise = sndosc_4p5o (-x) -sndosc_6p5o :: (Ord a, Fractional a) => a -> a +sndosc_6p5o :: (Ord a, Fractional a) => a -> a sndosc_6p5o x | 0 <= x && x < 1 = polyeval [ 1, 0, -5/4, -35/12, 21/4, -25/12 ] x | 1 <= x && x < 2 = polyeval [ -4, 75/4, -245/8, 545/24, -63/8, 25/24 ] x | 2 <= x && x < 3 = polyeval [ 18, -153/4, 255/8, -313/24, 21/8, -5/24 ] x @@ -141,13 +141,13 @@ -- Misc -watte_4p2o :: (Ord a, Fractional a) => a -> a +watte_4p2o :: (Ord a, Fractional a) => a -> a watte_4p2o x | 0 <= x && x < 1 = polyeval [ 1, -1/2, -1/2 ] x | 1 <= x && x < 2 = polyeval [ 1, -3/2, 1/2 ] x | 2 <= x = 0 | otherwise = watte_4p2o (-x) -parabolic2x_4p2o :: (Ord a, Fractional a) => a -> a +parabolic2x_4p2o :: (Ord a, Fractional a) => a -> a parabolic2x_4p2o x | 0 <= x && x < 1 = polyeval [ 1/2, 0, -1/4 ] x | 1 <= x && x < 2 = polyeval [ 1, -1, 1/4 ] x | 2 <= x = 0 @@ -157,292 +157,292 @@ -- Optimal designs -optimal_2p3o2x :: (Ord a, Fractional a) => a -> a +optimal_2p3o2x :: (Ord a, Fractional a) => a -> a optimal_2p3o2x x | 0 <= x && x < 1 = polyeval [ 0.80607906469176971, 0.17594740788514596, -2.35977550974341630, 1.57015627178718420 ] x | 1 <= x = 0 | otherwise = optimal_2p3o2x (-x) -optimal_2p3o4x :: (Ord a, Fractional a) => a -> a -optimal_2p3o4x x | 0 <= x && x < 1 = polyeval [ 0.88207975731800936, -0.10012219395448523, +optimal_2p3o4x :: (Ord a, Fractional a) => a -> a +optimal_2p3o4x x | 0 <= x && x < 1 = polyeval [ 0.88207975731800936, -0.10012219395448523, -1.99054787320203810, 1.32598918957298410 ] x | 1 <= x = 0 | otherwise = optimal_2p3o4x (-x) -optimal_2p3o8x :: (Ord a, Fractional a) => a -> a -optimal_2p3o8x x | 0 <= x && x < 1 = polyeval [ 0.94001491168487883, -0.51213628865925998, +optimal_2p3o8x :: (Ord a, Fractional a) => a -> a +optimal_2p3o8x x | 0 <= x && x < 1 = polyeval [ 0.94001491168487883, -0.51213628865925998, -1.10319974084152170, 0.73514591836770027 ] x | 1 <= x = 0 | otherwise = optimal_2p3o8x (-x) -optimal_2p3o16x :: (Ord a, Fractional a) => a -> a -optimal_2p3o16x x | 0 <= x && x < 1 = polyeval [ 0.96964782067188493, -0.74617479745643256, +optimal_2p3o16x :: (Ord a, Fractional a) => a -> a +optimal_2p3o16x x | 0 <= x && x < 1 = polyeval [ 0.96964782067188493, -0.74617479745643256, -0.57923093055631791, 0.38606621963374965 ] x | 1 <= x = 0 | otherwise = optimal_2p3o16x (-x) -optimal_2p3o32x :: (Ord a, Fractional a) => a -> a -optimal_2p3o32x x | 0 <= x && x < 1 = polyeval [ 0.98472017575676363, -0.87053863725307623, +optimal_2p3o32x :: (Ord a, Fractional a) => a -> a +optimal_2p3o32x x | 0 <= x && x < 1 = polyeval [ 0.98472017575676363, -0.87053863725307623, -0.29667081825572522, 0.19775766248673177 ] x | 1 <= x = 0 | otherwise = optimal_2p3o32x (-x) -optimal_4p2o2x :: (Ord a, Fractional a) => a -> a -optimal_4p2o2x x | 0 <= x && x < 1 = polyeval [ 0.50061662213752656, -0.04782068534965925, +optimal_4p2o2x :: (Ord a, Fractional a) => a -> a +optimal_4p2o2x x | 0 <= x && x < 1 = polyeval [ 0.50061662213752656, -0.04782068534965925, -0.21343978756177684 ] x - | 1 <= x && x < 2 = polyeval [ 0.92770135528027386, -0.88689658749623701, + | 1 <= x && x < 2 = polyeval [ 0.92770135528027386, -0.88689658749623701, 0.21303593243799016 ] x | 2 <= x = 0 | otherwise = optimal_4p2o2x (-x) -optimal_4p2o4x :: (Ord a, Fractional a) => a -> a -optimal_4p2o4x x | 0 <= x && x < 1 = polyeval [ 0.33820365736567115, 0.2114449807519728, +optimal_4p2o4x :: (Ord a, Fractional a) => a -> a +optimal_4p2o4x x | 0 <= x && x < 1 = polyeval [ 0.33820365736567115, 0.2114449807519728, -0.22865399531858188 ] x - | 1 <= x && x < 2 = polyeval [ 1.12014639874555470, -1.01414466618792900, + | 1 <= x && x < 2 = polyeval [ 1.12014639874555470, -1.01414466618792900, 0.22858390767180370 ] x | 2 <= x = 0 | otherwise = optimal_4p2o4x (-x) -optimal_4p2o8x :: (Ord a, Fractional a) => a -> a -optimal_4p2o8x x | 0 <= x && x < 1 = polyeval [ 0.09224718574204172, 0.59257579283164508, +optimal_4p2o8x :: (Ord a, Fractional a) => a -> a +optimal_4p2o8x x | 0 <= x && x < 1 = polyeval [ 0.09224718574204172, 0.59257579283164508, -0.24005206207889518 ] x - | 1 <= x && x < 2 = polyeval [ 1.38828036063664320, -1.17126532964206100, + | 1 <= x && x < 2 = polyeval [ 1.38828036063664320, -1.17126532964206100, 0.24004281672637814 ] x | 2 <= x = 0 | otherwise = optimal_4p2o8x (-x) -optimal_4p2o16x :: (Ord a, Fractional a) => a -> a -optimal_4p2o16x x | 0 <= x && x < 1 = polyeval [ -0.41849525763976203, 1.36361593203840510, +optimal_4p2o16x :: (Ord a, Fractional a) => a -> a +optimal_4p2o16x x | 0 <= x && x < 1 = polyeval [ -0.41849525763976203, 1.36361593203840510, -0.24506117865474364 ] x | 1 <= x && x < 2 = polyeval [ 1.90873339502208310, -1.44144384373471430, 0.24506002360805534 ] x | 2 <= x = 0 | otherwise = optimal_4p2o16x (-x) -optimal_4p2o32x :: (Ord a, Fractional a) => a -> a -optimal_4p2o32x x | 0 <= x && x < 1 = polyeval [ -1.42170796824052890, 2.87083485132510450, +optimal_4p2o32x :: (Ord a, Fractional a) => a -> a +optimal_4p2o32x x | 0 <= x && x < 1 = polyeval [ -1.42170796824052890, 2.87083485132510450, -0.24755243839713828 ] x | 1 <= x && x < 2 = polyeval [ 2.91684291662070860, -1.95043794419108290, 0.24755229501840223 ] x | 2 <= x = 0 | otherwise = optimal_4p2o32x (-x) -optimal_4p3o2x :: (Ord a, Fractional a) => a -> a -optimal_4p3o2x x | 0 <= x && x < 1 = polyeval [ 0.59244492420272321, 0.03573669883299365, +optimal_4p3o2x :: (Ord a, Fractional a) => a -> a +optimal_4p3o2x x | 0 <= x && x < 1 = polyeval [ 0.59244492420272321, 0.03573669883299365, -0.78664888597764893, 0.36030925263849456 ] x - | 1 <= x && x < 2 = polyeval [ 1.20220428331406090, -1.60101160971478710, + | 1 <= x && x < 2 = polyeval [ 1.20220428331406090, -1.60101160971478710, 0.70401463131621556, -0.10174985775982505 ] x | 2 <= x = 0 | otherwise = optimal_4p3o2x (-x) -optimal_4p3o4x :: (Ord a, Fractional a) => a -> a -optimal_4p3o4x x | 0 <= x && x < 1 = polyeval [ 0.60304009430474115, 0.05694012453786401, +optimal_4p3o4x :: (Ord a, Fractional a) => a -> a +optimal_4p3o4x x | 0 <= x && x < 1 = polyeval [ 0.60304009430474115, 0.05694012453786401, -0.89223007211175309, 0.42912649274763925 ] x | 1 <= x && x < 2 = polyeval [ 1.31228823423882930, -1.85072890189700660, 0.87687351895686727, -0.13963062613760227 ] x | 2 <= x = 0 | otherwise = optimal_4p3o4x (-x) -optimal_4p3o8x :: (Ord a, Fractional a) => a -> a -optimal_4p3o8x x | 0 <= x && x < 1 = polyeval [ 0.60658368706046584, 0.07280793921972525, +optimal_4p3o8x :: (Ord a, Fractional a) => a -> a +optimal_4p3o8x x | 0 <= x && x < 1 = polyeval [ 0.60658368706046584, 0.07280793921972525, -0.95149675410360302, 0.46789242171187317 ] x - | 1 <= x && x < 2 = polyeval [ 1.35919815911169020, -1.95618744839533010, + | 1 <= x && x < 2 = polyeval [ 1.35919815911169020, -1.95618744839533010, 0.94949311590826524, -0.15551896027602030 ] x | 2 <= x = 0 | otherwise = optimal_4p3o8x (-x) -optimal_4p3o16x :: (Ord a, Fractional a) => a -> a -optimal_4p3o16x x | 0 <= x && x < 1 = polyeval [ 0.60844825096346644, 0.07980169577604959, +optimal_4p3o16x :: (Ord a, Fractional a) => a -> a +optimal_4p3o16x x | 0 <= x && x < 1 = polyeval [ 0.60844825096346644, 0.07980169577604959, -0.97894238166068270, 0.48601256046234864 ] x - | 1 <= x && x < 2 = polyeval [ 1.37724137476464990, -1.99807048591354810, + | 1 <= x && x < 2 = polyeval [ 1.37724137476464990, -1.99807048591354810, 0.97870442828560433, -0.16195131297091253 ] x | 2 <= x = 0 | otherwise = optimal_4p3o16x (-x) -optimal_4p3o32x :: (Ord a, Fractional a) => a -> a -optimal_4p3o32x x | 0 <= x && x < 1 = polyeval [ 0.60908264223655417, 0.08298544053689563, +optimal_4p3o32x :: (Ord a, Fractional a) => a -> a +optimal_4p3o32x x | 0 <= x && x < 1 = polyeval [ 0.60908264223655417, 0.08298544053689563, -0.99052586766084594, 0.49369595780454456 ] x | 1 <= x && x < 2 = polyeval [ 1.38455689452848450, -2.01496368680360890, 0.99049753216621961, -0.16455902278580614 ] x | 2 <= x = 0 | otherwise = optimal_4p3o32x (-x) -optimal_4p4o2x :: (Ord a, Fractional a) => a -> a -optimal_4p4o2x x | 0 <= x && x < 1 = polyeval [ 0.58448510036125145, 0.04442540676862300, - -0.7586487041827807, 0.29412762852131868, +optimal_4p4o2x :: (Ord a, Fractional a) => a -> a +optimal_4p4o2x x | 0 <= x && x < 1 = polyeval [ 0.58448510036125145, 0.04442540676862300, + -0.7586487041827807, 0.29412762852131868, 0.04252164479749607 ] x - | 1 <= x && x < 2 = polyeval [ 1.06598379704160570, -1.16581445347275190, - 0.21256821036268256, 0.13781898240764315, + | 1 <= x && x < 2 = polyeval [ 1.06598379704160570, -1.16581445347275190, + 0.21256821036268256, 0.13781898240764315, -0.04289144034653719 ] x | 2 <= x = 0 | otherwise = optimal_4p4o2x (-x) -optimal_4p4o4x :: (Ord a, Fractional a) => a -> a -optimal_4p4o4x x | 0 <= x && x < 1 = polyeval [ 0.61340295990566229, 0.06128937679587994, - -0.94057832565094635, 0.44922093286355397, +optimal_4p4o4x :: (Ord a, Fractional a) => a -> a +optimal_4p4o4x x | 0 <= x && x < 1 = polyeval [ 0.61340295990566229, 0.06128937679587994, + -0.94057832565094635, 0.44922093286355397, 0.00986988334359864 ] x - | 1 <= x && x < 2 = polyeval [ 1.30835018075821670, -1.82814511658458520, - 0.81943257721092366, -0.09642760567543440, + | 1 <= x && x < 2 = polyeval [ 1.30835018075821670, -1.82814511658458520, + 0.81943257721092366, -0.09642760567543440, -0.00989340017126506 ] x | 2 <= x = 0 | otherwise = optimal_4p4o4x (-x) -optimal_4p4o8x :: (Ord a, Fractional a) => a -> a -optimal_4p4o8x x | 0 <= x && x < 1 = polyeval [ 0.62095991632974834, 0.06389302461261143, +optimal_4p4o8x :: (Ord a, Fractional a) => a -> a +optimal_4p4o8x x | 0 <= x && x < 1 = polyeval [ 0.62095991632974834, 0.06389302461261143, -0.98489647972932193, 0.48698871865064902, 0.00255074537015887 ] x | 1 <= x && x < 2 = polyeval [ 1.35943398999940390, -1.97277963497287720, - 0.95410568622888214, -0.14868053358928229, + 0.95410568622888214, -0.14868053358928229, -0.00255226912537286 ] x | 2 <= x = 0 | otherwise = optimal_4p4o8x (-x) -optimal_4p4o16x :: (Ord a, Fractional a) => a -> a -optimal_4p4o16x x | 0 <= x && x < 1 = polyeval [ 0.62293049365660191, 0.06443376638262904, - -0.99620011474430481, 0.49672182806667398, +optimal_4p4o16x :: (Ord a, Fractional a) => a -> a +optimal_4p4o16x x | 0 <= x && x < 1 = polyeval [ 0.62293049365660191, 0.06443376638262904, + -0.99620011474430481, 0.49672182806667398, 0.00064264050033187 ] x - | 1 <= x && x < 2 = polyeval [ 1.37216269878963180, -2.00931632449031920, - 0.98847675044522398, -0.16214364417487748, + | 1 <= x && x < 2 = polyeval [ 1.37216269878963180, -2.00931632449031920, + 0.98847675044522398, -0.16214364417487748, -0.00064273459469381 ] x | 2 <= x = 0 | otherwise = optimal_4p4o16x (-x) -optimal_4p4o32x :: (Ord a, Fractional a) => a -> a -optimal_4p4o32x x | 0 <= x && x < 1 = polyeval [ 0.62342449465938121, 0.06456923251842608, - -0.99904509583176049, 0.49917660509564427, +optimal_4p4o32x :: (Ord a, Fractional a) => a -> a +optimal_4p4o32x x | 0 <= x && x < 1 = polyeval [ 0.62342449465938121, 0.06456923251842608, + -0.99904509583176049, 0.49917660509564427, 0.00016095224137360 ] x - | 1 <= x && x < 2 = polyeval [ 1.37534629142898650, -2.01847637982642340, - 0.99711292321092770, -0.16553360612350931, + | 1 <= x && x < 2 = polyeval [ 1.37534629142898650, -2.01847637982642340, + 0.99711292321092770, -0.16553360612350931, -0.00016095810460478 ] x | 2 <= x = 0 | otherwise = optimal_4p4o32x (-x) -optimal_6p4o2x :: (Ord a, Fractional a) => a -> a -optimal_6p4o2x x | 0 <= x && x < 1 = polyeval [ 0.42640922432669054, -0.0052558029434142, - -0.20486985491012843, 0.00255494211547300, +optimal_6p4o2x :: (Ord a, Fractional a) => a -> a +optimal_6p4o2x x | 0 <= x && x < 1 = polyeval [ 0.42640922432669054, -0.0052558029434142, + -0.20486985491012843, 0.00255494211547300, 0.03134095684084392 ] x - | 1 <= x && x < 2 = polyeval [ 0.30902529029941583, 0.37868437559565432, - -0.70564644117967990, 0.31182026815653541, + | 1 <= x && x < 2 = polyeval [ 0.30902529029941583, 0.37868437559565432, + -0.70564644117967990, 0.31182026815653541, -0.04385804833432710 ] x - | 2 <= x && x < 3 = polyeval [ 1.51897639740576910, -1.83761742915820410, - 0.83217835730406542, -0.16695522597587154, + | 2 <= x && x < 3 = polyeval [ 1.51897639740576910, -1.83761742915820410, + 0.83217835730406542, -0.16695522597587154, 0.01249475765486819 ] x | 3 <= x = 0 | otherwise = optimal_6p4o2x (-x) -optimal_6p4o4x :: (Ord a, Fractional a) => a -> a -optimal_6p4o4x x | 0 <= x && x < 1 = polyeval [ 0.20167941634921072, -0.06119274485321008, - 0.56468711069379207, -0.42059475673758634, +optimal_6p4o4x :: (Ord a, Fractional a) => a -> a +optimal_6p4o4x x | 0 <= x && x < 1 = polyeval [ 0.20167941634921072, -0.06119274485321008, + 0.56468711069379207, -0.42059475673758634, 0.02881527997393852 ] x - | 1 <= x && x < 2 = polyeval [ -0.64579641436229407, 2.33580825807694700, - -1.85350543411307390, 0.51926458031522660, + | 1 <= x && x < 2 = polyeval [ -0.64579641436229407, 2.33580825807694700, + -1.85350543411307390, 0.51926458031522660, -0.04250898918476453 ] x - | 2 <= x && x < 3 = polyeval [ 2.76228852293285200, -3.09936092833253300, - 1.27147464005834010, -0.22283280665600644, + | 2 <= x && x < 3 = polyeval [ 2.76228852293285200, -3.09936092833253300, + 1.27147464005834010, -0.22283280665600644, 0.01369173779618459 ] x | 3 <= x = 0 | otherwise = optimal_6p4o4x (-x) -optimal_6p4o8x :: (Ord a, Fractional a) => a -> a -optimal_6p4o8x x | 0 <= x && x < 1 = polyeval [ -0.17436452172055789, -0.15190225510786248, - 1.87551558979819120, -1.15976496200057480, +optimal_6p4o8x :: (Ord a, Fractional a) => a -> a +optimal_6p4o8x x | 0 <= x && x < 1 = polyeval [ -0.17436452172055789, -0.15190225510786248, + 1.87551558979819120, -1.15976496200057480, 0.03401038103941584 ] x - | 1 <= x && x < 2 = polyeval [ -2.26955357035241170, 5.73320660746477540, - -3.92391712129699590, 0.93463067895166918, + | 1 <= x && x < 2 = polyeval [ -2.26955357035241170, 5.73320660746477540, + -3.92391712129699590, 0.93463067895166918, -0.05090907029392906 ] x - | 2 <= x && x < 3 = polyeval [ 4.84834508915762540, -5.25661448354449060, - 2.04584149450148180, -0.32814290420019698, + | 2 <= x && x < 3 = polyeval [ 4.84834508915762540, -5.25661448354449060, + 2.04584149450148180, -0.32814290420019698, 0.01689861603514873 ] x | 3 <= x = 0 | otherwise = optimal_6p4o8x (-x) -optimal_6p4o16x :: (Ord a, Fractional a) => a -> a -optimal_6p4o16x x | 0 <= x && x < 1 = polyeval [ -0.94730014688427577, -0.33649680079382827, - 4.53807483241466340, -2.64598691215356660, +optimal_6p4o16x :: (Ord a, Fractional a) => a -> a +optimal_6p4o16x x | 0 <= x && x < 1 = polyeval [ -0.94730014688427577, -0.33649680079382827, + 4.53807483241466340, -2.64598691215356660, 0.03755086455339280 ] x - | 1 <= x && x < 2 = polyeval [ -5.55035312316726960, 12.52871168241192600, - -7.98288364772738750, 1.70665858343069510, + | 1 <= x && x < 2 = polyeval [ -5.55035312316726960, 12.52871168241192600, + -7.98288364772738750, 1.70665858343069510, -0.05631219122315393 ] x - | 2 <= x && x < 3 = polyeval [ 8.94785524286246310, -9.37021675593126700, - 3.44447036756440590, -0.49470749109917245, + | 2 <= x && x < 3 = polyeval [ 8.94785524286246310, -9.37021675593126700, + 3.44447036756440590, -0.49470749109917245, 0.01876132424143207 ] x | 3 <= x = 0 | otherwise = optimal_6p4o16x (-x) -optimal_6p4o32x :: (Ord a, Fractional a) => a -> a -optimal_6p4o32x x | 0 <= x && x < 1 = polyeval [ -2.44391738331193720, -0.69468212315980082, - 9.67889243081689440, -5.50592307590218160, +optimal_6p4o32x :: (Ord a, Fractional a) => a -> a +optimal_6p4o32x x | 0 <= x && x < 1 = polyeval [ -2.44391738331193720, -0.69468212315980082, + 9.67889243081689440, -5.50592307590218160, 0.03957507923965987 ] x - | 1 <= x && x < 2 = polyeval [ -11.87524595267807600, 25.58633277328986500, - -15.73068663442630400, 3.15288929279855570, + | 1 <= x && x < 2 = polyeval [ -11.87524595267807600, 25.58633277328986500, + -15.73068663442630400, 3.15288929279855570, -0.05936083498715066 ] x - | 2 <= x && x < 3 = polyeval [ 16.79403235763479100, -17.17264148794549100, - 6.05175140696421730, -0.79053754554850286, + | 2 <= x && x < 3 = polyeval [ 16.79403235763479100, -17.17264148794549100, + 6.05175140696421730, -0.79053754554850286, 0.01978575568000696 ] x | 3 <= x = 0 | otherwise = optimal_6p4o32x (-x) -optimal_6p5o2x :: (Ord a, Fractional a) => a -> a -optimal_6p5o2x x | 0 <= x && x < 1 = polyeval [ 0.48217702203158502, -0.00127577239632662, - -0.3267507171395277, -0.02014846731685776, +optimal_6p5o2x :: (Ord a, Fractional a) => a -> a +optimal_6p5o2x x | 0 <= x && x < 1 = polyeval [ 0.48217702203158502, -0.00127577239632662, + -0.3267507171395277, -0.02014846731685776, 0.14640674192652170, -0.04317950185225609 ] x - | 1 <= x && x < 2 = polyeval [ 0.35095903476754237, 0.53534756396439365, - -1.22477236472789920, 0.74995484587342742, + | 1 <= x && x < 2 = polyeval [ 0.35095903476754237, 0.53534756396439365, + -1.22477236472789920, 0.74995484587342742, -0.19234043023690772, 0.01802814255926417 ] x - | 2 <= x && x < 3 = polyeval [ 1.62814578813495040, -2.26168360510917840, - 1.22220278720010690, -0.31577407091450355, + | 2 <= x && x < 3 = polyeval [ 1.62814578813495040, -2.26168360510917840, + 1.22220278720010690, -0.31577407091450355, 0.03768876199398620, -0.00152170021558204 ] x | 3 <= x = 0 | otherwise = optimal_6p5o2x (-x) -optimal_6p5o4x :: (Ord a, Fractional a) => a -> a -optimal_6p5o4x x | 0 <= x && x < 1 = polyeval [ 0.50164509338655083, -0.00256790184606694, - -0.36229943140977111, -0.04512026308730401, +optimal_6p5o4x :: (Ord a, Fractional a) => a -> a +optimal_6p5o4x x | 0 <= x && x < 1 = polyeval [ 0.50164509338655083, -0.00256790184606694, + -0.36229943140977111, -0.04512026308730401, 0.20620318519804220, -0.06607747864416924 ] x - | 1 <= x && x < 2 = polyeval [ 0.30718330223223800, 0.78336433172501685, - -1.66940481896969310, 1.08365113099941970, + | 1 <= x && x < 2 = polyeval [ 0.30718330223223800, 0.78336433172501685, + -1.66940481896969310, 1.08365113099941970, -0.30560854964737405, 0.03255079211953620 ] x - | 2 <= x && x < 3 = polyeval [ 2.05191571792256240, -3.19403437421534920, - 1.99766476840488070, -0.62765808573554227, + | 2 <= x && x < 3 = polyeval [ 2.05191571792256240, -3.19403437421534920, + 1.99766476840488070, -0.62765808573554227, 0.09909173357642603, -0.00628989632244913 ] x | 3 <= x = 0 | otherwise = optimal_6p5o4x (-x) -optimal_6p5o8x :: (Ord a, Fractional a) => a -> a -optimal_6p5o8x x | 0 <= x && x < 1 = polyeval [ 0.50513183702821474, -0.00368143670114908, - -0.36434084624989699, -0.06070462616102962, +optimal_6p5o8x :: (Ord a, Fractional a) => a -> a +optimal_6p5o8x x | 0 <= x && x < 1 = polyeval [ 0.50513183702821474, -0.00368143670114908, + -0.36434084624989699, -0.06070462616102962, 0.22942797169644802, -0.07517133281176167 ] x - | 1 <= x && x < 2 = polyeval [ 0.28281884957695946, 0.88385964850687193, - -1.82581238657617080, 1.19588167464050650, + | 1 <= x && x < 2 = polyeval [ 0.28281884957695946, 0.88385964850687193, + -1.82581238657617080, 1.19588167464050650, -0.34363487882262922, 0.03751837438141215 ] x - | 2 <= x && x < 3 = polyeval [ 2.15756386503245070, -3.42137079071284810, - 2.18592382088982260, -0.70370361187427199, + | 2 <= x && x < 3 = polyeval [ 2.15756386503245070, -3.42137079071284810, + 2.18592382088982260, -0.70370361187427199, 0.11419603882898799, -0.00747588873055296 ] x | 3 <= x = 0 | otherwise = optimal_6p5o8x (-x) -optimal_6p5o16x :: (Ord a, Fractional a) => a -> a -optimal_6p5o16x x | 0 <= x && x < 1 = polyeval [ 0.50819303579369868, -0.00387117789818541, - -0.36990908725555449, -0.06616250180411522, +optimal_6p5o16x :: (Ord a, Fractional a) => a -> a +optimal_6p5o16x x | 0 <= x && x < 1 = polyeval [ 0.50819303579369868, -0.00387117789818541, + -0.36990908725555449, -0.06616250180411522, 0.24139298776307896, -0.07990500783668089 ] x - | 1 <= x && x < 2 = polyeval [ 0.27758734130911511, 0.91870010875159547, - -1.89281840112089440, 1.24834464824612510, + | 1 <= x && x < 2 = polyeval [ 0.27758734130911511, 0.91870010875159547, + -1.89281840112089440, 1.24834464824612510, -0.36203450650610985, 0.03994519162531633 ] x - | 2 <= x && x < 3 = polyeval [ 2.19284545406407450, -3.50786533926449100, - 2.26228244623301580, -0.73559668875725392, + | 2 <= x && x < 3 = polyeval [ 2.19284545406407450, -3.50786533926449100, + 2.26228244623301580, -0.73559668875725392, 0.12064126711558003, -0.00798609327859495 ] x | 3 <= x = 0 | otherwise = optimal_6p5o16x (-x) -optimal_6p5o32x :: (Ord a, Fractional a) => a -> a -optimal_6p5o32x x | 0 <= x && x < 1 = polyeval [ 0.52558916128536759, 0.00010896283126635, - -0.42682321682847008, -0.04095676092513167, +optimal_6p5o32x :: (Ord a, Fractional a) => a -> a +optimal_6p5o32x x | 0 <= x && x < 1 = polyeval [ 0.52558916128536759, 0.00010896283126635, + -0.42682321682847008, -0.04095676092513167, 0.25041444762720882, -0.08349799235675044 ] x - | 1 <= x && x < 2 = polyeval [ 0.33937904183610190, 0.80946953063234006, - -1.86228986389877100, 1.27215033630638800, + | 1 <= x && x < 2 = polyeval [ 0.33937904183610190, 0.80946953063234006, + -1.86228986389877100, 1.27215033630638800, -0.37562266426589430, 0.04174912841630993 ] x - | 2 <= x && x < 3 = polyeval [ 2.13606003964474490, -3.48774662195185850, - 2.28912105276248390, -0.75510203509083995, + | 2 <= x && x < 3 = polyeval [ 2.13606003964474490, -3.48774662195185850, + 2.28912105276248390, -0.75510203509083995, 0.12520821766375972, -0.00834987866042734 ] x | 3 <= x = 0 | otherwise = optimal_6p5o32x (-x) @@ -507,10 +507,10 @@ interpolate y h = sum $ zipWith (*) y (elems h) -x1 = sin $ 0.345 + 0.1234 * 1.2 +x1 = sin $ 0.345 + 0.1234 * 1.2 x1' = map (interpolate y) [ h1, h2, h3, h4, h5, h6, h7 ] -x2 = sin $ 0.345 + 0.1234 * 2.2 +x2 = sin $ 0.345 + 0.1234 * 2.2 x2' = map (interpolate y) [ h8, h9, h10, h11, h12 ] The values of all these lists should be one, or nearly one. They 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 @@ -9,9 +9,9 @@ -- Portability : portable -- -- Module to sharpen FIR filters --- +-- -- Reference: Hamming, Sect 6.6 --- +-- -- @H'(z) = 3 * H(z)^2 - s * H(z)^3@ -- @ = H(z)^2 * (3 - 2 * H(z))@ -- @@ -31,8 +31,7 @@ import Data.Array -import DSP.Basic -import DSP.Convolution +import qualified DSP.Basic as Basic import DSP.Filter.FIR.FIR -- | Filter shaprening routine @@ -43,7 +42,7 @@ sharpen h x = step4 where step1 = fir h x step2 = map (2*) step1 - step3 = zipWith (-) (map (3*) (zn delay x)) step2 + step3 = zipWith (-) (map (3*) (Basic.delay delay x)) step2 step4 = fir h $ fir h $ step3 -- step4 = fir $ conv h h $ step3 m = snd $ bounds h 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 @@ -10,7 +10,7 @@ -- -- Herrmann type smooth FIR filters, from Hamming, Chapter 7, also -- known as maximally flat FIR filters --- +-- -- If x is the -3 dB point, then p\/q = -(x+1)\/(x-1) -- ----------------------------------------------------------------------------- @@ -27,20 +27,24 @@ -- Normalize is the step to set g(1) = 1 (pg 123) +normalize :: Fractional a => [a] -> [a] normalize x = map (/ a) x where a = sum x -- Expand performs the algorithm in Sect 7.3 +expand :: Fractional a => [a] -> [a] expand (x1:x2:[]) = [ x1, x2 ] expand (x:xs) = expand' x $ expand xs -expand' x ys = zipWith (+) (m1 x ys) (p1 ys) - where m1 x (y:ys) = x : y : map (0.5*) ys - p1 (y:ys) = map (0.5*) ys ++ [ 0, 0 ] +expand' :: Fractional a => a -> [a] -> [a] +expand' x ys0 = zipWith (+) (x : m1 ys0) (p1 ys0) + 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 :: Fractional a => [a] -> [a] reflect (x:xs) = (map (0.5*) $ reverse xs) ++ x : (map (0.5*) xs) -- The actual function. Note that we use (1+t)^p * (1-t)^q directly 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 @@ -32,6 +32,7 @@ -- indexes generates the list of indexes that we will map the prototype -- functions onto +indexes :: (Integral a, Num b, Enum b) => a -> [b] indexes m = [ 0 .. fromIntegral m ] -- the _tap functions generate one tap for the given function @@ -42,32 +43,37 @@ -- Lowpass tap function +lpf_tap :: (Integral a, Floating 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 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 (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" -- Raised-cosine tap function. This does _not_ have 0 dB DC gain. -- ws = symbol rate in normalized radians -- b = filter beta +rc_tap :: (Integral a, Floating a1) => a1 -> a1 -> a -> a1 -> a1 rc_tap ws b m n | n-a == 0 = 1 | den == 0 = 0 | otherwise = sin sarg / sarg * cos carg / den where sarg = ws * (n-a) / 2 carg = b * ws * (n-a) / 2 - den = 1 - 4 * ((b*ws*(n-a)) / (2*pi)) ^ 2 + den = 1 - 4 * ((b*ws*(n-a)) / (2*pi)) ^ (2::Int) a = (fromIntegral m) / 2 -- The following functions generate a list of the taps for a given set of 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 @@ -25,11 +25,11 @@ Reference: @Book{dsp, - author = "Alan V. Oppenheim and Ronald W. Schafer", - title = "Discrete-Time Signal Processing", - publisher = "Pretice-Hall", - year = 1989, - address = "Englewood Cliffs", + 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} } @@ -44,11 +44,14 @@ -} -module DSP.Filter.FIR.Window (window, rectangular, bartlett, hanning, hamming, blackman, - kaiser, gen_hamming, parzen) where +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] @@ -61,78 +64,90 @@ -- | rectangular window rectangular :: Int -- ^ M - -> Array Int Double -- ^ w[n] + -> Array Int Double -- ^ w[n] -rectangular m = listArray (0,m) $ replicate (m+1) 1.0 +rectangular m = listArray (0,m) $ repeat 1 -- | Bartlett window bartlett :: Int -- ^ M - -> Array Int Double -- ^ w[n] + -> Array Int Double -- ^ w[n] -bartlett m = listArray (0,m) $ map (bartlett' md) [ 0.0 .. md ] - where bartlett' m n | n <= m / 2 = 2 * n / m - | otherwise = 2 - 2 * n / m - md = fromIntegral m +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] + -> Array Int Double -- ^ w[n] -hanning m = listArray (0,m) $ map (hanning' md) [ 0.0 .. md ] - where hanning' m n = 0.5 - 0.5 * cos(2 * pi * n / m) - md = fromIntegral m +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] + -> Array Int Double -- ^ w[n] -hamming m = listArray (0,m) $ map (hamming' md) [ 0.0 .. md ] - where hamming' m n = 0.54 - 0.46 * cos(2 * pi * n / m) - md = fromIntegral m +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] + -> Array Int Double -- ^ w[n] -blackman m = listArray (0,m) $ map (blackman' md) [ 0.0 .. md ] - where blackman' m n = 0.42 - 0.5 * cos(2 * pi * n / m) + - 0.08 * cos (4 * pi * n / m) - md = fromIntegral m +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] + -> Int -- ^ M + -> Array Int Double -- ^ w[n] -gen_hamming a m = listArray (0,m) $ map (hamming' a md) [ 0.0 .. md ] - where hamming' a m n = a - (1 - a) * cos(2 * pi * n / m) - md = fromIntegral m +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 b m = listArray (0,m) $ map (kaiser' b md) [ 0.0 .. md ] - where kaiser' b m n = i0 (b * sqrt (1 -((n-a)/a)^2)) / i0 b - md = fromIntegral m - a = md / 2 +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' :: 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)) + | 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... @@ -142,7 +157,16 @@ parzen :: Int -- ^ M -> Array Int Double -- ^ w[n] -parzen m = listArray (0,m) $ map (parzen' md) [ 0.0 .. md ] - where parzen' m n | n <= m / 2 = 2 * (1-n/m) ^ 3 - (1-2*n/m) ^ 3 - | otherwise = 2 * (1-n/m) ^ 3 - md = fromIntegral m +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/DSP/Filter/IIR/Bilinear.hs b/DSP/Filter/IIR/Bilinear.hs --- a/DSP/Filter/IIR/Bilinear.hs +++ b/DSP/Filter/IIR/Bilinear.hs @@ -9,21 +9,21 @@ -- Portability : portable -- -- The module contains a function for performing the bilinear transform. --- +-- -- The input is a rational polynomial representation of the s-domain -- function to be transformed. --- +-- -- In the bilinear transform, we substitute --- --- @ 2 1 - z^-1@ -- +-- @ 2 1 - z^-1@ +-- -- @s \<-- -- * --------@ -- -- @ ts 1 + z^-1@ --- +-- -- into the rational polynomial, where ts is the sampling period. To get -- a rational polynomial back, we use the following method: --- +-- -- (1) Substitute s^n with (2\/ts * (1-z^-1))^n == [ -2\/ts, 2\/ts ]^n -- -- 2. Multiply the results by (1+z^-1)^n == [ 1, 1 ]^n @@ -31,7 +31,7 @@ -- 3. Add up all of the common terms -- -- 4. Normalize all of the coeficients by a0 --- +-- -- where n is the maximum order of the numerator and denominator -- ----------------------------------------------------------------------------- @@ -41,36 +41,42 @@ -- TODO: Do we want to include prewarping? -module DSP.Filter.IIR.Bilinear (bilinear, prewarp) where +module DSP.Filter.IIR.Bilinear {- (bilinear, prewarp) -} where import Polynomial.Basic -- Computes (2\/ts * (1-z^-1))^n == [ -2\/ts, 2\/ts ]^n +zm :: (Integral b, Fractional a) => a -> b -> [a] zm ts n = polypow [ -2/ts, 2/ts ] n -- Computes (1+z^-1)^n == [ 1, 1 ]^n +zp :: (Integral b, Num a) => b -> [a] zp n = polypow [ 1, 1 ] n -- Step 1: Substitute s^n with (2\/ts * (1-z^-1))^n == [ -2\/ts, 2\/ts ]^n -- in num and den -step1 ts x = step1' ts 0 x - where step1' _ _ [] = [] - step1' ts n (x:xs) = map (x*) (zm ts n) : step1' ts (n+1) xs +step1 :: Fractional a => a -> [a] -> [[a]] +step1 ts = step1' (0::Int) + where step1' _ [] = [] + step1' n (x:xs) = polyscale x (zm ts n) : step1' (n+1) xs -- Step 2: Multiply the num and den by (1+z^-1)^n == [ 1, 1 ]^n +step2 :: (Num a, Integral b) => b -> [[a]] -> [[a]] step2 _ [] = [] step2 n (x:xs) = polymult (zp n) x : step2 (n-1) xs -- Step 3: Add up all of the common terms +step3 :: Num a => [[a]] -> [a] step3 x = foldr polyadd [0] x -- Step 4: Normalize all of the coeficients by a0 +step4 :: Fractional a => a -> [a] -> [a] step4 a0 x = map (/a0) x -- Glue it all together diff --git a/DSP/Filter/IIR/Cookbook.lhs b/DSP/Filter/IIR/Cookbook.lhs --- a/DSP/Filter/IIR/Cookbook.lhs +++ b/DSP/Filter/IIR/Cookbook.lhs @@ -27,6 +27,8 @@ > import DSP.Filter.IIR.IIR +> import DSP.Basic((^!)) + Cookbook formulae for audio EQ biquad filter coefficients ----------------------------------------------------------------------------- by Robert Bristow-Johnson @@ -300,18 +302,18 @@ a1 = -2*[ (A-1) + (A+1)*cos ] a2 = (A+1) + (A-1)*cos - beta*sin -> {-# specialize lowShelf :: Float -> Float -> Float -> Float -> [Float] -> [Float] #-} -> {-# specialize lowShelf :: Double -> Double -> Double -> Double -> [Double] -> [Double] #-} +> {-# specialize lowShelf :: Float -> Float -> Float -> [Float] -> [Float] #-} +> {-# specialize lowShelf :: Double -> Double -> Double -> [Double] -> [Double] #-} -> lowShelf :: Floating a => a -> a -> a -> a -> [a] -> [a] -> lowShelf bw s dBgain w = biquad_df1 (a1/a0) (a2/a0) (b0/a0) (b1/a0) (b2/a0) +> lowShelf :: Floating a => a -> a -> a -> [a] -> [a] +> lowShelf s dBgain w = biquad_df1 (a1/a0) (a2/a0) (b0/a0) (b1/a0) (b2/a0) > where b0 = a*( (a+1) - (a-1) * cos w + beta * sin w) > b1 = 2*a*( (a-1) - (a+1) * cos w ) > b2 = a*( (a+1) - (a-1) * cos w - beta * sin w) > a0 = (a+1) + (a-1) * cos w + beta * sin w > a1 = -2*( (a-1) + (a+1) * cos w ) > a2 = (a+1) + (a-1) * cos w - beta * sin w -> beta = sqrt ((a^2 + 1) / s - (a-1)^2) +> beta = sqrt ((a^!2 + 1) / s - (a-1)^!2) > a = 10 ** (dBgain / 40) highShelf: H(s) = A * (A*s^2 + (sqrt(A)/Q)*s + 1) / (s^2 + (sqrt(A)/Q)*s + A) @@ -323,18 +325,18 @@ a1 = 2*[ (A-1) - (A+1)*cos ] a2 = (A+1) - (A-1)*cos - beta*sin -> {-# specialize highShelf :: Float -> Float -> Float -> Float -> [Float] -> [Float] #-} -> {-# specialize highShelf :: Double -> Double -> Double -> Double -> [Double] -> [Double] #-} +> {-# specialize highShelf :: Float -> Float -> Float -> [Float] -> [Float] #-} +> {-# specialize highShelf :: Double -> Double -> Double -> [Double] -> [Double] #-} -> highShelf :: Floating a => a -> a -> a -> a -> [a] -> [a] -> highShelf bw s dBgain w = biquad_df1 (a1/a0) (a2/a0) (b0/a0) (b1/a0) (b2/a0) +> highShelf :: Floating a => a -> a -> a -> [a] -> [a] +> highShelf s dBgain w = biquad_df1 (a1/a0) (a2/a0) (b0/a0) (b1/a0) (b2/a0) > where b0 = a*( (a+1) - (a-1) * cos w + beta * sin w) > b1 = -2*a*( (a-1) - (a+1) * cos w ) > b2 = a*( (a+1) - (a-1) * cos w - beta * sin w) > a0 = (a+1) + (a-1) * cos w + beta * sin w > a1 = -2*( (a-1) + (a+1) * cos w ) > a2 = (a+1) + (a-1) * cos w - beta * sin w -> beta = sqrt ((a^2 + 1) / s - (a-1)^2) +> beta = sqrt ((a^!2 + 1) / s - (a-1)^!2) > a = 10 ** (dBgain / 40) (This text-only file is best viewed or printed with a mono-spaced font.) diff --git a/DSP/Filter/IIR/Design.hs b/DSP/Filter/IIR/Design.hs --- a/DSP/Filter/IIR/Design.hs +++ b/DSP/Filter/IIR/Design.hs @@ -22,12 +22,15 @@ module DSP.Filter.IIR.Design where -import Data.Array - import DSP.Filter.Analog.Prototype import DSP.Filter.Analog.Transform import DSP.Filter.IIR.Bilinear +import DSP.Basic ((^!)) + +import Data.Array + +poly2iir :: ([a], [b]) -> (Array Int a, Array Int b) poly2iir (b,a) = (b',a') where b' = listArray (0,m) $ reverse $ b a' = listArray (0,n) $ reverse $ a @@ -40,12 +43,12 @@ -> (Double, Double) -- ^ (ws,ds) -> (Array Int Double, Array Int Double) -- ^ (b,a) -mkButterworth (wp,dp) (ws,ds) = poly2iir $ - bilinear 1 $ - a_lp2lp wc $ +mkButterworth (wp,dp) (ws,ds) = poly2iir $ + bilinear 1 $ + a_lp2lp wc $ butterworth n - where n = ceiling $ log (((1/ds)^2-1) / ((1/(1-dp))^2-1)) / 2 / log (ws' / wp') - wc = ws' / ((1/ds)^2-1)**(1/2/fromIntegral n) + where n = ceiling $ log (((1/ds)^!2-1) / ((1/(1-dp))^!2-1)) / 2 / log (ws' / wp') + wc = ws' / ((1/ds)^!2-1)**(1/2/fromIntegral n) wp' = prewarp wp 1 ws' = prewarp ws 1 @@ -55,17 +58,17 @@ -> (Double, Double) -- ^ (ws,ds) -> (Array Int Double, Array Int Double) -- ^ (b,a) -mkChebyshev1 (wp,dp) (ws,ds) = poly2iir $ - bilinear 1 $ - a_lp2lp wp' $ +mkChebyshev1 (wp,dp) (ws,ds) = poly2iir $ + bilinear 1 $ + a_lp2lp wp' $ chebyshev1 eps n where wp' = prewarp wp 1 ws' = prewarp ws 1 eps = sqrt ((2 - dp)*dp) / (1 - dp) a = 1 / ds - k1 = eps / sqrt (a^2 - 1) + k1 = eps / sqrt (a^!2 - 1) k = wp' / ws' - n = ceiling $ acosh (1/k1) / log ((1 + sqrt (1 - k^2)) / k) + n = ceiling $ acosh (1/k1) / log ((1 + sqrt (1 - k^!2)) / k) -- | Generates lowpass Inverse Chebyshev IIR filters @@ -73,12 +76,12 @@ -> (Double, Double) -- ^ (ws,ds) -> (Array Int Double, Array Int Double) -- ^ (b,a) -mkChebyshev2 (wp,dp) (ws,ds) = poly2iir $ - bilinear 1 $ - a_lp2lp ws' $ +mkChebyshev2 (wp,dp) (ws,ds) = poly2iir $ + bilinear 1 $ + a_lp2lp ws' $ chebyshev2 eps n where wp' = prewarp wp 1 ws' = prewarp ws 1 - eps = ds / sqrt (1 - ds^2) + eps = ds / sqrt (1 - ds^!2) g = 1 - dp - n = ceiling $ acosh (g / eps / sqrt (1 - g^2)) / acosh (ws' / wp') + n = ceiling $ acosh (g / eps / sqrt (1 - g^!2)) / acosh (ws' / wp') 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 @@ -15,9 +15,9 @@ -- Except in integrator, we use the convention that -- -- @y[n] = sum(k=0..M) b_k*x[n-k] - sum(k=1..N) a_k*y[n-k]@ --- -- -- +-- -- @ sum(k=0..M) b_k*z^-1@ -- -- @H(z) = ------------------------@ @@ -47,16 +47,19 @@ -} module DSP.Filter.IIR.IIR (integrator, - fos_df1, fos_df2, fos_df2t, - biquad_df1, biquad_df2, biquad_df2t, - iir_df1, iir_df2) where + fos_df1, fos_df2, fos_df2t, + biquad_df1, biquad_df2, biquad_df2t, + iir_df1, iir_df2, + -- for testing + xt, yt, f1, f2, f3, f4, f5, + ) where import Data.Array import DSP.Filter.FIR.FIR -- | This is an integrator when a==1, and a leaky integrator when @0 \< a \< 1@. --- +-- -- @y[n] = a * y[n-1] + x[n]@ {-# specialize integrator :: Float -> [Float] -> [Float] #-} @@ -96,8 +99,8 @@ fos_df1' a1 b0 b1 x1 y1 (x:xs) = y : fos_df1' a1 b0 b1 x y xs where v = b0 * x + b1 * x1 y = v - a1 * y1 - + -- | First order section, DF2 -- -- @w[n] = -a1 * w[n-1] + x[n]@ @@ -137,7 +140,7 @@ -> a -- ^ b_1 -> [a] -- ^ x[n] -> [a] -- ^ y[n] - + fos_df2t a1 b0 b1 x = fos_df2t' a1 b0 b1 0 x fos_df2t' :: Num a => a -> a -> a -> a -> [a] -> [a] @@ -171,8 +174,8 @@ df1 a1 a2 b0 b1 b2 x1 x2 y1 y2 (x:xs) = y : df1 a1 a2 b0 b1 b2 x x1 y y1 xs where v = b0 * x + b1 * x1 + b2 * x2 y = v - a1 * y1 - a2 * y2 - + -- | Direct Form II for a second order section (biquad) -- -- @w[n] = -a1 * w[n-1] - a2 * w[n-2] + x[n]@ @@ -218,7 +221,7 @@ -> a -- ^ b_2 -> [a] -- ^ x[n] -> [a] -- ^ y[n] - + biquad_df2t a1 a2 b0 b1 b2 x = df2t a1 a2 b0 b1 b2 0 0 x df2t :: Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a] @@ -254,7 +257,7 @@ {- specialize iir'df1 :: Array Int Double -> Array Int Double -> [Double] -> [Double] -} iir'df1 :: (Num a) => Array Int a -> Array Int a -> [a] -> [a] -iir'df1 a w [] = [] +iir'df1 _ _ [] = [] iir'df1 a w (v:vs) = y : iir'df1 a w' vs where y = v - sum [ a!i * w!i | i <- [1..n] ] w' = listArray (1,n) $ y : elems w @@ -284,32 +287,41 @@ {- specialize iir'df2 :: Array Int Double -> Array Int Double -> [Double] -> [Double] -} iir'df2 :: (Num a) => (Array Int a,Array Int a) -> Array Int a -> [a] -> [a] -iir'df2 (b,a) w [] = [] +iir'df2 _ _ [] = [] iir'df2 (b,a) w (x:xs) = y : iir'df2 (b,a) w' xs where y = sum [ b!i * w'!i | i <- [0..m] ] w0 = x - sum [ a!i * w'!i | i <- [1..m] ] w' = listArray (0,mn) $ w0 : elems w m = snd $ bounds b - n = snd $ bounds a mn = snd $ bounds w --------- -- test -x = [ 1, 0, 0, 0, 0, 0, 0, 0 ] :: [Double] +xt :: [Double] +xt = [ 1, 0, 0, 0, 0, 0, 0, 0 ] :: [Double] -y = integrator 0.5 x +yt :: [Double] +yt = integrator 0.5 xt +f1 :: Fractional a => [a] -> [a] f1 x = biquad_df1 (-0.4) 0.3 0.5 0.4 (-0.3) x +f2 :: Fractional a => [a] -> [a] f2 x = biquad_df2 (-0.4) 0.3 0.5 0.4 (-0.3) x +f3 :: Fractional a => [a] -> [a] f3 x = biquad_df2t (-0.4) 0.3 0.5 0.4 (-0.3) x -a = listArray (1,2) [ -0.4, 0.3 ] -b = listArray (0,2) [ 0.5, 0.4, -0.3 ] +at :: Array Int Double +at = listArray (1,2) [ -0.4, 0.3 ] -f4 x = iir_df1 (b,a) x +bt :: Array Int Double +bt = listArray (0,2) [ 0.5, 0.4, -0.3 ] -f5 x = iir_df2 (b,a) x +f4 :: [Double] -> [Double] +f4 x = iir_df1 (bt,at) x + +f5 :: [Double] -> [Double] +f5 x = iir_df2 (bt,at) x diff --git a/DSP/Filter/IIR/Prony.hs b/DSP/Filter/IIR/Prony.hs --- a/DSP/Filter/IIR/Prony.hs +++ b/DSP/Filter/IIR/Prony.hs @@ -9,7 +9,7 @@ -- Portability : portable -- -- General case of Prony's Method where K > p+q --- +-- -- References: L&I, Sect 8.1; P&B, Sect 7.5; P&M, Sect 8.5.2 -- -- Notation follows L&I @@ -29,7 +29,7 @@ {------------------------------------------------------------------------------ -Case 1: K=p+q +Case 1: K=p+q a = array (0,p) b = array (0,q) @@ -55,7 +55,7 @@ ------------------------------------------------------------------------------} --- Case 2: K>p+q +-- Case 2: K>p+q -- a = array (0,p) -- b = array (0,q) diff --git a/DSP/Filter/IIR/Transform.hs b/DSP/Filter/IIR/Transform.hs --- a/DSP/Filter/IIR/Transform.hs +++ b/DSP/Filter/IIR/Transform.hs @@ -23,83 +23,61 @@ module DSP.Filter.IIR.Transform (d_lp2lp, d_lp2hp, d_lp2bp, d_lp2bs) where -import Data.Complex +import DSP.Filter.Analog.Transform (substitute) +import Numeric.Special.Trigonometric (cot) -import Polynomial.Basic -normalize :: ([Double],[Double]) -> ([Double],[Double]) -normalize (num,den) = (num',den') - where a0 = last den - num' = map (/ a0) num - den' = map (/ a0) den - -substitute :: ([Double],[Double]) -> ([Double],[Double]) -> ([Double],[Double]) -substitute (nsub,dsub) (num,den) = normalize (num',den') - where n = max (length num - 1) (length den - 1) - num' = step3 $ step2 0 dsub $ step1 n nsub $ num - den' = step3 $ step2 0 dsub $ step1 n nsub $ den - step1 _ _ [] = [] - step1 n w (x:xs) = map (x*) (polypow w n) : step1 (n-1) w xs - step2 _ _ [] = [] - step2 n w (x:xs) = polymult (polypow w n) x : step2 (n+1) w xs - step3 x = foldr polyadd [0] x - --- Cotangent - -cot :: Double -> Double -cot x = 1 / tan x - -- | Lowpass to lowpass: @z^-1 --> (z^-1 - a)\/(1 - a*z^-1)@ d_lp2lp :: Double -- ^ theta_p - -> Double -- ^ omega_p - -> ([Double], [Double]) -- ^ (b,a) - -> ([Double], [Double]) -- ^ (b',a') + -> Double -- ^ omega_p + -> ([Double], [Double]) -- ^ (b,a) + -> ([Double], [Double]) -- ^ (b',a') d_lp2lp tp wp (num,den) = substitute (nsub,dsub) (num,den) - where nsub = [1, -a] - dsub = [-a, 1] - a = sin ((tp-wp)/2) / sin ((tp+wp)/2) + where nsub = [-a, 1] + dsub = [1, -a] + a = sin ((tp-wp)/2) / sin ((tp+wp)/2) -- | Lowpass to Highpass: @z^-1 --> -(z^-1 + a)\/(1 + a*z^-1)@ d_lp2hp :: Double -- ^ theta_p - -> Double -- ^ omega_p - -> ([Double], [Double]) -- ^ (b,a) - -> ([Double], [Double]) -- ^ (b',a') + -> Double -- ^ omega_p + -> ([Double], [Double]) -- ^ (b,a) + -> ([Double], [Double]) -- ^ (b',a') d_lp2hp tp wp (num,den) = substitute (nsub,dsub) (num,den) - where nsub = [-1, -a] - dsub = [a, 1] - a = -cos ((tp+wp)/2) / cos ((tp-wp)/2) + where nsub = [a, 1] + dsub = [-1, -a] + a = -cos ((tp+wp)/2) / cos ((tp-wp)/2) --- | Lowpass to Bandpass: z^-1 --> +-- | Lowpass to Bandpass: z^-1 --> d_lp2bp :: Double -- ^ theta_p - -> Double -- ^ omega_p1 - -> Double -- ^ omega_p2 - -> ([Double], [Double]) -- ^ (b,a) - -> ([Double], [Double]) -- ^ (b',a') + -> Double -- ^ omega_p1 + -> Double -- ^ omega_p2 + -> ([Double], [Double]) -- ^ (b,a) + -> ([Double], [Double]) -- ^ (b',a') d_lp2bp tp wp1 wp2 (num,den) = substitute (nsub,dsub) (num,den) - where nsub = [ 1, -2*a*k/(k+1), (k-1)/(k+1) ] - dsub = [ (k-1)/(k+1), -2*a*k/(k+1), 1 ] - a = cos ((wp2+wp1)/2) / cos ((wp2-wp1)/2) - k = cot ((wp2-wp1)/2) * tan (tp/2) + where nsub = [ (k-1)/(k+1), -2*a*k/(k+1), 1 ] + dsub = [ 1, -2*a*k/(k+1), (k-1)/(k+1) ] + a = cos ((wp2+wp1)/2) / cos ((wp2-wp1)/2) + k = cot ((wp2-wp1)/2) * tan (tp/2) --- | Lowpass to Bandstop: z^-1 --> +-- | Lowpass to Bandstop: z^-1 --> d_lp2bs :: Double -- ^ theta_p - -> Double -- ^ omega_p1 - -> Double -- ^ omega_p2 - -> ([Double], [Double]) -- ^ (b,a) - -> ([Double], [Double]) -- ^ (b',a') + -> Double -- ^ omega_p1 + -> Double -- ^ omega_p2 + -> ([Double], [Double]) -- ^ (b,a) + -> ([Double], [Double]) -- ^ (b',a') d_lp2bs tp wp1 wp2 (num,den) = substitute (nsub,dsub) (num,den) - where nsub = [ 1, -2*a/(1+k), (1-k)/(1+k) ] - dsub = [ (1-k)/(1+k), -2*a/(1+k), 1 ] - a = cos ((wp2+wp1)/2) / cos ((wp2-wp1)/2) - k = cot ((wp2-wp1)/2) * tan (tp/2) + where nsub = [ (1-k)/(1+k), -2*a/(1+k), 1 ] + dsub = [ 1, -2*a/(1+k), (1-k)/(1+k) ] + a = cos ((wp2+wp1)/2) / cos ((wp2-wp1)/2) + k = cot ((wp2-wp1)/2) * tan (tp/2) {- diff --git a/DSP/Multirate/CIC.hs b/DSP/Multirate/CIC.hs --- a/DSP/Multirate/CIC.hs +++ b/DSP/Multirate/CIC.hs @@ -43,9 +43,9 @@ module DSP.Multirate.CIC (cic_interpolate, cic_decimate) where -import DSP.Basic +import DSP.Basic (delay1, delay, upsample, downsample) --- apply returns a function of n applications of a function, eg, +-- apply returns a function of n applications of a function, eg, -- apply f 3 = f . f . f @@ -61,13 +61,13 @@ -- integrate [ 1, 1, 1, 1 ] = [ 1, 2, 3, 4 ] integrate :: (Num a) => [a] -> [a] -integrate a = zipWith (+) a (z (integrate a)) +integrate a = zipWith (+) a (delay1 (integrate a)) -- comb implements the comb function described in the paper above. The m -- parameter is the length of the delay in the feed-forward element. comb :: (Num a) => Int -> [a] -> [a] -comb m a = zipWith (-) a (zn m a) +comb m a = zipWith (-) a (delay m a) {- diff --git a/DSP/Multirate/Halfband.hs b/DSP/Multirate/Halfband.hs --- a/DSP/Multirate/Halfband.hs +++ b/DSP/Multirate/Halfband.hs @@ -18,35 +18,23 @@ import Data.Array -import DSP.Basic +import DSP.Basic (delay, uninterleave, interleave) import DSP.Filter.FIR.FIR mkhalfband :: Num a => Array Int a -> Array Int a -mkhalfband h = listArray (0,m `div` 2) [ h!n | n <- [0..m], even n ] +mkhalfband h = listArray (0,m `div` 2) [ h!n | n <- [0,2..m] ] where m = snd $ bounds h -demux :: Num a => [a] -> ([a],[a]) -demux (x:xs) = (demux' (x:xs), demux' xs) - where demux' [] = [] - demux' (x:[]) = x : [] - demux' (x:_:xs) = x : demux' xs - -mux :: Num a => [a] -> [a] -> [a] -mux [] [] = [] -mux [] _ = [] -mux _ [] = [] -mux (x:xs) (y:ys) = x : y : mux xs ys - -- | Halfband interpolator hb_interp :: (Num a) => Array Int a -- ^ h[n] -> [a] -- ^ x[n] -> [a] -- ^ y[n] -hb_interp h x = mux y1 y2 - where (x1,x2) = demux x +hb_interp h x = interleave y1 y2 + where (x1,x2) = uninterleave x y1 = fir (mkhalfband h) x1 - y2 = map (h!m2 *) $ zn m2 $ x2 + y2 = map (h!m2 *) $ delay m2 $ x2 m2 = (snd $ bounds h) `div` 2 -- | Halfband decimator @@ -56,7 +44,7 @@ -> [a] -- ^ y[n] hb_decim h x = zipWith (+) y1 y2 - where (x1,x2) = demux x + where (x1,x2) = uninterleave x y1 = fir (mkhalfband h) x1 - y2 = map (h!m2 *) $ zn m2 $ x2 + y2 = map (h!m2 *) $ delay m2 $ x2 m2 = (snd $ bounds h) `div` 2 diff --git a/DSP/Multirate/Polyphase.hs b/DSP/Multirate/Polyphase.hs --- a/DSP/Multirate/Polyphase.hs +++ b/DSP/Multirate/Polyphase.hs @@ -16,6 +16,7 @@ module DSP.Multirate.Polyphase (poly_interp) where +import Data.List (transpose) import Data.Array import DSP.Filter.FIR.FIR @@ -33,11 +34,9 @@ -> [a] -- ^ x[n] -> [a] -- ^ y[n] -poly_interp l h x = commutate y +poly_interp l h x = concat $ transpose y where g = map (fir . mkpoly h l) [0..(l-1)] - y = map (\f -> f x) g - commutate [] = [] - commutate xs = [h | (h:t) <- xs] ++ commutate [t | (h:t) <- xs] + y = map ($ x) g {- diff --git a/DSP/Source/Basic.hs b/DSP/Source/Basic.hs --- a/DSP/Source/Basic.hs +++ b/DSP/Source/Basic.hs @@ -17,7 +17,7 @@ -- | all zeros zeros :: (Num a) => [a] -zeros = 0 : zeros +zeros = repeat 0 -- | single impulse @@ -27,9 +27,9 @@ -- | unit step step :: (Num a) => [a] -step = 1 : step +step = repeat 1 -- | ramp ramp :: (Num a) => [a] -ramp = 0 : zipWith (+) ramp (repeat 1) +ramp = iterate (1+) 0 diff --git a/DSP/Source/Oscillator.hs b/DSP/Source/Oscillator.hs --- a/DSP/Source/Oscillator.hs +++ b/DSP/Source/Oscillator.hs @@ -12,15 +12,16 @@ -- ----------------------------------------------------------------------------- -module DSP.Source.Oscillator (nco, ncom, - quadrature_nco, complex_ncom, - quadrature_ncom) where +module DSP.Source.Oscillator (nco, ncom, + quadrature_nco, complex_ncom, + quadrature_ncom, + agc) where import Data.Complex -- | 'nco' creates a sine wave with normalized frequency wn (numerically -- controlled oscillator, or NCO) using the recurrence relation y[n] = --- 2cos(wn)*y[n-1] - y[n-2]. Eventually, cumlative errors will creep +-- 2cos(wn)*y[n-1] - y[n-2]. Eventually, cumulative errors will creep -- into the data. This is unavoidable since performing AGC on this type -- of real data is hard. The good news is that the error is small with -- floating point data. @@ -55,7 +56,7 @@ -- (3-x)/2 for x ~= 1 to eliminate doing a sqrt for every point. agc :: RealFloat a => Complex a -> Complex a -agc z@(x:+y) = x * r :+ y * r +agc (x:+y) = x * r :+ y * r where r = (3 - x * x - y * y) / 2 -- | 'quadrature_nco' returns an infinite list representing a complex phasor diff --git a/DSP/Unwrap.hs b/DSP/Unwrap.hs --- a/DSP/Unwrap.hs +++ b/DSP/Unwrap.hs @@ -28,7 +28,7 @@ -> Array a b -- ^ arg unwrap eps phi = listArray b [ phi!i + 2 * pi * r!i | i <- range b ] - where r = listArray b [ ri i | i <- range b ] + where r = listArray b [ ri i | i <- range b ] ri 0 = 0 ri i | phi!i - phi!(i-1) > (2*pi-eps) = r!(i-1) - 1 | phi!i - phi!(i-1) < -(2*pi-eps) = r!(i-1) + 1 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/Matrix/Cholesky.hs b/Matrix/Cholesky.hs --- a/Matrix/Cholesky.hs +++ b/Matrix/Cholesky.hs @@ -32,5 +32,5 @@ l = array ((1,1),(n,n)) [ ((i,j), lij i j) | i <- [2..n], j <- [1..(i-1)] ] lij i j | j==1 = a!(i,1) / d!1 | otherwise = a!(i,j) / d!j - sum [ l!(i,k) * d!k * (conjugate (l!(j,k))) / d!j | k <- [1..(j-1)] ] - d = array (1,n) ((1, a!(1,1)) : [ (i, a!(i,i) - sum [ d!k * ((abs (l!(i,k)))^2) | k <- [1..(i-1)] ] ) | i <- [2..n]]) + d = array (1,n) ((1, a!(1,1)) : [ (i, a!(i,i) - sum [ d!k * ((abs (l!(i,k)))^(2::Int)) | k <- [1..(i-1)] ] ) | i <- [2..n]]) ((_,_),(n,_)) = bounds a diff --git a/Matrix/LU.hs b/Matrix/LU.hs --- a/Matrix/LU.hs +++ b/Matrix/LU.hs @@ -12,7 +12,7 @@ -- ----------------------------------------------------------------------------- -module Matrix.LU (lu, lu_solve, improve, inverse, lu_det, solve, det) where +module Matrix.LU (lu, lu_solve, improve, inverse, lu_det, solve, det) where import Data.Array @@ -29,9 +29,11 @@ lu a = a' where a' = array bnds [ ((i,j), luij i j) | (i,j) <- range bnds ] - luij i j | i>j = (a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(j-1)] ]) / a'!(j,j) - | i<=j = a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(i-1)] ] - bnds = bounds a + luij i j = + if i>j + then (a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(j-1)] ]) / a'!(j,j) + else a!(i,j) - sum [ a'!(i,k) * a'!(k,j) | k <- [1 ..(i-1)] ] + bnds = bounds a -- | Solution to Ax=b via LU decomposition @@ -39,32 +41,32 @@ -- backward is forumla (2.3.7) in NRIC lu_solve :: Array (Int,Int) Double -- ^ LU(A) - -> Array Int Double -- ^ b - -> Array Int Double -- ^ x + -> Array Int Double -- ^ b + -> Array Int Double -- ^ x lu_solve a b = x where x = array (1,n) ([(n,xn)] ++ [ (i, backward i) | i <- (reverse [1..(n-1)]) ]) - y = array (1,n) ([(1,y1)] ++ [ (i, forward i) | i <- [2..n] ]) - y1 = b!1 - forward i = (b!i - sum [ a!(i,j) * y!j | j <- [1..(i-1)] ]) - xn = y!n / a!(n,n) - backward i = (y!i - sum [ a!(i,j) * x!j | j <- [(i+1)..n] ]) / a!(i,i) - ((_,_),(n,_)) = bounds a + y = array (1,n) ([(1,y1)] ++ [ (i, forward i) | i <- [2..n] ]) + y1 = b!1 + forward i = (b!i - sum [ a!(i,j) * y!j | j <- [1..(i-1)] ]) + xn = y!n / a!(n,n) + backward i = (y!i - sum [ a!(i,j) * x!j | j <- [(i+1)..n] ]) / a!(i,i) + ((_,_),(n,_)) = bounds a -- | Improve a solution to Ax=b via LU decomposition -- formula (2.7.4) from NRIC improve :: Array (Int,Int) Double -- ^ A - -> Array (Int,Int) Double -- ^ LU(A) - -> Array Int Double -- ^ b - -> Array Int Double -- ^ x - -> Array Int Double -- ^ x' + -> Array (Int,Int) Double -- ^ LU(A) + -> Array Int Double -- ^ b + -> Array Int Double -- ^ x + -> Array Int Double -- ^ x' improve a a_lu b x = array (1,n) [ (i, x!i - err!i) | i <- [1..n] ] where err = lu_solve a_lu rhs - rhs = array (1,n) [ (i, sum [ a!(i,j) * x!j | j <- [1..n] ] - b!i) | i <- [1..n] ] - ((_,_),(n,_)) = bounds a + rhs = array (1,n) [ (i, sum [ a!(i,j) * x!j | j <- [1..n] ] - b!i) | i <- [1..n] ] + ((_,_),(n,_)) = bounds a -- | Matrix inversion via LU decomposition @@ -73,16 +75,16 @@ -- TODO: build in improve inverse :: Array (Int,Int) Double -- ^ A - -> Array (Int,Int) Double -- ^ A^-1 + -> Array (Int,Int) Double -- ^ A^-1 -inverse a = a' - where a' = array (bounds a) (arrange (makecols (lu a)) 1) - makecol i n = array (1,n) [ (j, (\i j->if i == j then 1.0 else 0.0) i j) | j <- [1..n] ] - makecols a = [ lu_solve a (makecol i n) | i <- [1..n] ] - ((_,_),(n,_)) = bounds a - arrange [] _ = [] - arrange (m:ms) j = (flatten m j) ++ (arrange ms (j+1)) - flatten m j = map (\(i,x) -> ((i,j),x)) (assocs m) +inverse a0 = a' + where a' = array (bounds a0) (arrange (makecols (lu a0)) 1) + makecol i n' = array (1,n') [ (j, if i == j then 1.0 else 0.0) | j <- [1..n'] ] + makecols a = [ lu_solve a (makecol i n) | i <- [1..n] ] + ((_,_),(n,_)) = bounds a0 + arrange [] _ = [] + arrange (m:ms) j = flatten m j ++ arrange ms (j+1) + flatten m j = map (\(i,x) -> ((i,j),x)) (assocs m) -- | Determinant of a matrix via LU decomposition @@ -98,8 +100,8 @@ -- | LU solver using original matrix solve :: Array (Int,Int) Double -- ^ A - -> Array Int Double -- ^ b - -> Array Int Double -- ^ x + -> Array Int Double -- ^ b + -> Array Int Double -- ^ x solve a b = (lu_solve . lu) a b @@ -117,11 +119,11 @@ {- a = array ((1,1),(3,3)) [ ((1,1), 1.0), ((1,2), 2.0), ((1,3), 3.0), - ((2,1), 2.0), ((2,2), 5.0), ((2,3), 3.0), - ((3,1), 1.0), ((3,2), 0.0), ((3,3), 8.0) ] + ((2,1), 2.0), ((2,2), 5.0), ((2,3), 3.0), + ((3,1), 1.0), ((3,2), 0.0), ((3,3), 8.0) ] a' = array ((1,1),(3,3)) [ ((1,1), -40.0), ((1,2), 16.0), ((1,3), 9.0), - ((2,1), 13.0), ((2,2), -5.0), ((2,3), -3.0), - ((3,1), 5.0), ((3,2), -2.0), ((3,3), -1.0) ] + ((2,1), 13.0), ((2,2), -5.0), ((2,3), -3.0), + ((3,1), 5.0), ((3,2), -2.0), ((3,3), -1.0) ] a_lu = lu a b = array (1,3) [ (1, 1.0), (2, 2.0), (3, 5.0) ] @@ -130,8 +132,8 @@ x'' = improve a a_lu b x' verify = a' == inverse a && -- tests lu, lu_solve, and inverse - det a == -1 && -- tests lu_det - x == x' && -- tests improve - x' == x'' + det a == -1 && -- tests lu_det + x == x' && -- tests improve + x' == x'' -} diff --git a/Matrix/Levinson.hs b/Matrix/Levinson.hs --- a/Matrix/Levinson.hs +++ b/Matrix/Levinson.hs @@ -39,8 +39,8 @@ ak 1 1 = -r!1 / r!0 ak k i | k==i = -(r!k + sum [ a!(k-1,l) * r!(k-l) | l <- [1..(k-1)] ]) / rho!(k-1) | otherwise = a!(k-1,i) + a!(k,k) * (conjugate (a!(k-1,k-i))) - rhok 1 = (1 - (abs (a!(1,1)))^2) * r!0 - rhok k = (1 - (abs (a!(k,k)))^2) * rho!(k-1) + rhok 1 = (1 - (abs (a!(1,1)))^(2::Int)) * r!0 + rhok k = (1 - (abs (a!(k,k)))^(2::Int)) * rho!(k-1) -- r = array (0,2) [ (0, (2.0 :+ 0.0)), (1, ((-1.0) :+ 1.0)), (2, (0.0 :+ 0.0)) ] -- a = fst (levinson r 2) diff --git a/Matrix/Matrix.hs b/Matrix/Matrix.hs --- a/Matrix/Matrix.hs +++ b/Matrix/Matrix.hs @@ -23,10 +23,10 @@ -> Array (a,a) b -- ^ B -> Array (a,a) b -- ^ C -mm_mult a b = if ac /= br +mm_mult a b = if ac /= br then error "mm_mult: inside dimensions inconsistent" else array bnds [ ((i,j), mult i j) | (i,j) <- range bnds ] - where mult i j = sum [ a!(i,k) * b!(k,j) | k <- [1..ac] ] + where mult i j = sum [ a!(i,k) * b!(k,j) | k <- [1..ac] ] ((_,_),(ar,ac)) = bounds a ((_,_),(br,bc)) = bounds b bnds = ((1,1),(ar,bc)) @@ -37,7 +37,7 @@ -> Array a b -- ^ b -> Array a b -- ^ c -mv_mult a b = if ac /= br +mv_mult a b = if ac /= br then error "mv_mult: dimensions inconsistent" else array bnds [ (i, mult i) | i <- range bnds ] where mult i = sum [ a!(i,k) * b!(k) | k <- [1..ac] ] diff --git a/Matrix/Simplex.hs b/Matrix/Simplex.hs --- a/Matrix/Simplex.hs +++ b/Matrix/Simplex.hs @@ -36,40 +36,47 @@ -- Pivot around a!(p,q) pivot :: Int -> Int -> Array (Int,Int) Double -> Array (Int,Int) Double -pivot p q a = step4 p q $ step3 p q $ step2 p q $ step1 p q $ a - where step1 p q a = a // [ ((j,k), a!(j,k) - a!(p,k) * a!(j,q) / a!(p,q)) | k <- [0..m], j <- [ph..n], j /= p && k /= q ] - step2 p q a = a // [ ((j,q),0) | j <- [ph..n], j /= p ] - step3 p q a = a // [ ((p,k), a!(p,k) / a!(p,q)) | k <- [0..m], k /= q ] - step4 p q a = a // [ ((p,q),1) ] - ((ph,_),(n,m)) = bounds a +pivot p q a0 = step4 $ step3 $ step2 $ step1 $ a0 + where step1 a = a // [ ((j,k), a!(j,k) - a!(p,k) * a!(j,q) / a!(p,q)) | k <- [0..m], j <- [ph..n], j /= p && k /= q ] + step2 a = a // [ ((j,q),0) | j <- [ph..n], j /= p ] + step3 a = a // [ ((p,k), a!(p,k) / a!(p,q)) | k <- [0..m], k /= q ] + step4 a = a // [ ((p,q),1) ] + ((ph,_),(n,m)) = bounds a0 -- chooseq picks the lowest numbered favorable column. If there are no -- favorable columns, then q==m is returned, and we have reached an -- optimum. -chooseq a = chooseq' 1 a - where chooseq' q a | q > m = q - | a!(0,q) < -eps = q - | otherwise = chooseq' (q+1) a - ((_,_),(n,m)) = bounds a +chooseq :: (Ord b, Num b, Ix a, Ix b, Num a) => + Array (a, b) Double -> b +chooseq a0 = chooseq' 1 a0 + where chooseq' q a | q > m = q + | a!(0,q) < -eps = q + | otherwise = chooseq' (q+1) a + ((_,_),(_,m)) = bounds a0 + -- choosep picks a row with a positive element in column q. If no such -- element exists, then the p==n is returned, and the problem is -- unfeasible. -choosep q a = choosep' 1 q a - where choosep' p q a | p > n = p - | a!(p,q) > eps = p - | otherwise = choosep' (p+1) q a - ((_,_),(n,m)) = bounds a +choosep :: (Ord a, Num a, Ix a, Ix b) => + b -> Array (a, b) Double -> a +choosep q a0 = choosep' 1 a0 + where choosep' p a | p > n = p + | a!(p,q) > eps = p + | otherwise = choosep' (p+1) a + ((_,_),(n,_)) = bounds a0 -- refinep picks the row using the ratio test. -refinep p q a = refinep' (p+1) p q a - where refinep' i p q a | i > n = p - | a!(i,q) > eps && a!(i,0) / a!(i,q) < a!(p,0) / a!(p,q) = refinep' (i+1) i q a - | otherwise = refinep' (i+1) p q a - ((_,_),(n,m)) = bounds a +refinep :: (Ord a, Num a, Num b, Ix a, Ix b) => + a -> b -> Array (a, b) Double -> a +refinep p0 q a0 = refinep' (p0+1) p0 a0 + where refinep' i p a | i > n = p + | a!(i,q) > eps && a!(i,0) / a!(i,q) < a!(p,0) / a!(p,q) = refinep' (i+1) i a + | otherwise = refinep' (i+1) p a + ((_,_),(n,_)) = bounds a0 -- * Types @@ -77,12 +84,10 @@ data Simplex a = Unbounded | Infeasible | Optimal a deriving (Read,Show) -gettab (Optimal a) = a - -- * Functions -- | The simplex algorithm for standard form: --- +-- -- min c'x -- -- where Ax = b, x >= 0 @@ -96,46 +101,58 @@ -- a!(i,j) = A_ij simplex :: Array (Int,Int) Double -- ^ stating tableau - -> Simplex (Array (Int,Int) Double) -- ^ solution + -> Simplex (Array (Int,Int) Double) -- ^ solution simplex a | q > m = Optimal a - | p > n = Unbounded + | p > n = Unbounded | otherwise = simplex $ pivot p' q $ a where q = chooseq a p = choosep q a p' = refinep p q a - ((_,_),(n,m)) = bounds a + ((_,_),(n,m)) = bounds a ------------------------------------------------------------------------------- +addart :: (Num e, Enum a, Ix a, Num a) => + Array (a, a) e -> Array (a, a) e addart a = array ((-1,0),(n,m+n)) $ z ++ xsi ++ b ++ art ++ x where z = ((-1,0), a!(0,0)) : [ ((-1,j),0) | j <- [1..n] ] ++ [ ((-1,j+n),a!(0,j)) | j <- [1..m] ] - xsi = ((0,0), -colsum a 0) : [ ((0,j),0) | j <- [1..n] ] ++ [ ((0,j+n), -colsum a j) | j <- [1..m] ] - b = [ ((i,0), a!(i,0)) | i <- [1..n] ] - art = [ ((i,j), if i == j then 1 else 0) | i <- [1..n], j <- [1..n] ] - x = [ ((i,j+n), a!(i,j)) | i <- [1..n], j <- [1..m] ] + xsi = ((0,0), -colsum a 0) : [ ((0,j),0) | j <- [1..n] ] ++ [ ((0,j+n), -colsum a j) | j <- [1..m] ] + b = [ ((i,0), a!(i,0)) | i <- [1..n] ] + art = [ ((i,j), if i == j then 1 else 0) | i <- [1..n], j <- [1..n] ] + x = [ ((i,j+n), a!(i,j)) | i <- [1..n], j <- [1..m] ] ((_,_),(n,m)) = bounds a +colsum :: (Num e, Num a, Enum a, Ix a, Ix b) => + Array (a, b) e -> b -> e colsum a j = sum [ a!(i,j) | i <- [1..n] ] - where ((_,_),(n,m)) = bounds a + where ((_,_),(n,_)) = bounds a +delart :: (Enum a, Ix a, Num a) => + Array (a, a) e -> Array (a, a) e -> Array (a, a) e delart a a' = array ((0,0),(n,m)) $ z ++ b ++ x where z = ((0,0), a'!(-1,0)) : [ ((0,j), a!(0,j)) | j <- [1..m] ] - b = [ ((i,0), a'!(i,0)) | i <- [1..n] ] - x = [ ((i,j), a'!(i,j+n)) | i <- [1..n], j <- [1..m] ] + b = [ ((i,0), a'!(i,0)) | i <- [1..n] ] + x = [ ((i,j), a'!(i,j+n)) | i <- [1..n], j <- [1..m] ] ((_,_),(n,m)) = bounds a -- | The two-phase simplex algorithm twophase :: Array (Int,Int) Double -- ^ stating tableau - -> Simplex (Array (Int,Int) Double) -- ^ solution + -> Simplex (Array (Int,Int) Double) -- ^ solution twophase a | cost a' > eps = Infeasible | otherwise = simplex $ delart a (gettab a') where a' = simplex $ addart $ a + +-- How to handle cases where 'simplex' does not return Optimal? +gettab :: Simplex a -> a +gettab (Optimal a) = a + +cost :: (Num e, Ix a, Ix b, Num a, Num b) => + Simplex (Array (a, b) e) -> e cost (Optimal a) = negate $ a!(0,0) - where ((_,_),(n,m)) = bounds a ------------------------------------------------------------------------------- @@ -146,40 +163,40 @@ This is from Sedgewick > x1 = listArray ((0,0),(5,8)) [ 0, -1, -1, -1, 0, 0, 0, 0, 0, -> 5, -1, 1, 0, 1, 0, 0, 0, 0, -> 45, 1, 4, 0, 0, 1, 0, 0, 0, -> 27, 2, 1, 0, 0, 0, 1, 0, 0, -> 24, 3, -4, 0, 0, 0, 0, 1, 0, -> 4, 0, 0, 1, 0, 0, 0, 0, 1 ] :: Array (Int,Int) Double +> 5, -1, 1, 0, 1, 0, 0, 0, 0, +> 45, 1, 4, 0, 0, 1, 0, 0, 0, +> 27, 2, 1, 0, 0, 0, 1, 0, 0, +> 24, 3, -4, 0, 0, 0, 0, 1, 0, +> 4, 0, 0, 1, 0, 0, 0, 0, 1 ] :: Array (Int,Int) Double P&S, Example 2.6 > x2 = listArray ((0,0),(3,5)) [ 0, 1, 1, 1, 1, 1, -> 1, 3, 2, 1, 0, 0, -> 3, 5, 1, 1, 1, 0, -> 4, 2, 5, 1, 0, 1 ] :: Array (Int,Int) Double +> 1, 3, 2, 1, 0, 0, +> 3, 5, 1, 1, 1, 0, +> 4, 2, 5, 1, 0, 1 ] :: Array (Int,Int) Double P&S, Example 2.6 (after BFS selection) > x2' = listArray ((0,0),(3,5)) [ -6, -3, -3, 0, 0, 0, -> 1, 3, 2, 1, 0, 0, -> 2, 2, -1, 0, 1, 0, -> 3, -1, 3, 0, 0, 1 ] :: Array (Int,Int) Double +> 1, 3, 2, 1, 0, 0, +> 2, 2, -1, 0, 1, 0, +> 3, -1, 3, 0, 0, 1 ] :: Array (Int,Int) Double P&S, Example 2.2 / Section 2.9 > x3 = listArray ((0,0),(4,7)) [ -34, -1, -14, -6, 0, 0, 0, 0, -> 4, 1, 1, 1, 1, 0, 0, 0, -> 2, 1, 0, 0, 0, 1, 0, 0, -> 3, 0, 0, 1, 0, 0, 1, 0, -> 6, 0, 3, 1, 0, 0, 0, 1 ] :: Array (Int,Int) Double +> 4, 1, 1, 1, 1, 0, 0, 0, +> 2, 1, 0, 0, 0, 1, 0, 0, +> 3, 0, 0, 1, 0, 0, 1, 0, +> 6, 0, 3, 1, 0, 0, 0, 1 ] :: Array (Int,Int) Double P&S, Example 2.7 > x4 = listArray ((0,0),(3,7)) [ 3, -3/4, 20, -1/2, 6, 0, 0, 0, -> 0, 1/4, -8, -1, 9, 1, 0, 0, -> 0, 1/2, -12, -1/2, 3, 0, 1, 0, -> 1, 0, 0, 1, 0, 0, 0, 1 ] :: Array (Int,Int) Double +> 0, 1/4, -8, -1, 9, 1, 0, 0, +> 0, 1/2, -12, -1/2, 3, 0, 1, 0, +> 1, 0, 0, 1, 0, 0, 0, 1 ] :: Array (Int,Int) Double These come in handy for testing @@ -194,8 +211,8 @@ > find a j = findone' a 1 j > where findone' a i j | i > n = 0 -> | a!(i,j) == 1.0 = b!i -> | otherwise = findone' a (i+1) j +> | a!(i,j) == 1.0 = b!i +> | otherwise = findone' a (i+1) j > b = listArray (1,n) [ a!(i,0) | i <- [1..n] ] > ((_,_),(n,m)) = bounds a diff --git a/Numeric/Random/Distribution/Binomial.hs b/Numeric/Random/Distribution/Binomial.hs --- a/Numeric/Random/Distribution/Binomial.hs +++ b/Numeric/Random/Distribution/Binomial.hs @@ -28,7 +28,7 @@ -> Double -- ^ p -> [Double] -- ^ U -> [Double] -- ^ X - + binomial n p us = sum xi : binomial n p (drop n us) where xi = map (\u -> if u < p then 1 else 0) (take n us) diff --git a/Numeric/Random/Distribution/Gamma.hs b/Numeric/Random/Distribution/Gamma.hs --- a/Numeric/Random/Distribution/Gamma.hs +++ b/Numeric/Random/Distribution/Gamma.hs @@ -30,7 +30,7 @@ -> Double -- ^ lambda -> [Double] -- ^ U -> [Double] -- ^ X - + gamma n lambda u = x : gamma n lambda u' where x = -log (product (take n u)) / lambda u' = drop n u diff --git a/Numeric/Random/Distribution/Geometric.hs b/Numeric/Random/Distribution/Geometric.hs --- a/Numeric/Random/Distribution/Geometric.hs +++ b/Numeric/Random/Distribution/Geometric.hs @@ -29,6 +29,6 @@ geometric :: Double -- ^ p -> [Double] -- ^ U -> [Double] -- ^ X - + geometric p us = map (\u -> 1 + log u / log (1 - p)) us diff --git a/Numeric/Random/Distribution/Normal.hs b/Numeric/Random/Distribution/Normal.hs --- a/Numeric/Random/Distribution/Normal.hs +++ b/Numeric/Random/Distribution/Normal.hs @@ -21,9 +21,12 @@ -- TODO: Ahrens-Dieter method -module Numeric.Random.Distribution.Normal (normal_clt, normal_bm, - normal_ar, normal_r) where +module Numeric.Random.Distribution.Normal (normal_clt, normal_bm, + normal_ar, normal_r) where +import DSP.Basic (interleave, uninterleave, norm2sqr, toMaybe) +import Data.Maybe (mapMaybe) + -- * Functions -- adjust takes a unit normal random variable and sets the mean and @@ -39,66 +42,70 @@ -- [-n/2,n/2] normal_clt :: Int -- ^ Number of uniforms to sum - -> (Double,Double) -- ^ (mu,sigma) - -> [Double] -- ^ U - -> [Double] -- ^ X + -> (Double,Double) -- ^ (mu,sigma) + -> [Double] -- ^ U + -> [Double] -- ^ X -normal_clt n (mu,sigma) u = map (adjust (mu,sigma)) $ normal' u +normal_clt n muSigma u = map (adjust muSigma) $ normal' u where normal' us = var_adj * ((sum $ take n us) - mean_adj) : (normal' $ drop n us) - var_adj = sqrt $ 12 / fromIntegral n - mean_adj = fromIntegral n / 2 + var_adj = sqrt $ 12 / fromIntegral n + mean_adj = fromIntegral n / 2 -- | Normal random variables via the Box-Mueller Polar Method (Ross, pp -- 450--452) --- +-- -- If mu=0 and sigma=1, then this will generate numbers in the range -- [-8.57,8.57] assuing that the uniform RNG is really giving full -- precision for doubles. normal_bm :: (Double,Double) -- ^ (mu,sigma) - -> [Double] -- ^ U - -> [Double] -- ^ X + -> [Double] -- ^ U + -> [Double] -- ^ X -normal_bm (mu,sigma) u = map (adjust (mu,sigma)) $ normal' u - where normal' (u1:u2:us) | w <= 1 = x : y : normal' us - | otherwise = normal' us - where v1 = 2 * u1 - 1 - v2 = 2 * u2 - 1 - w = v1 * v1 + v2 * v2 - x = v1 * sqrt (-2 * log w / w) - y = v2 * sqrt (-2 * log w / w) +normal_bm muSigma = + map (adjust muSigma) . + uncurry interleave . unzip . mapMaybe normalDist . + uncurry zip . uninterleave . map (\u -> 2*u-1) +normalDist :: (Floating a, Ord a) => (a,a) -> Maybe (a,a) +normalDist z@(x,y) = + let norm2 = norm2sqr z + p = sqrt (-2 * log norm2) / norm2 + in toMaybe (norm2<=1) (p*x, p*y) + + -- | Acceptance-Rejection Method (Ross, pp 448--450) --- +-- -- If mu=0 and sigma=1, then this will generate numbers in the range -- [-36.74,36.74] assuming that the uniform RNG is really giving full -- precision for doubles. normal_ar :: (Double,Double) -- ^ (mu,sigma) - -> [Double] -- ^ U - -> [Double] -- ^ X + -> [Double] -- ^ U + -> [Double] -- ^ X -normal_ar (mu,sigma) u = map (adjust (mu,sigma)) $ normal' u +normal_ar muSigma u = map (adjust muSigma) $ normal' u where normal' (u1:u2:u3:us) | y > 0 = z : normal' us - | otherwise = normal' (u3:us) - where y1 = -log u1 - y2 = -log u2 - y = y2 - (y1 - 1)^2 / 2 - z | u3 <= 0.5 = y1 - | u3 > 0.5 = -y1 + | otherwise = normal' (u3:us) + where y1 = -log u1 + y2 = -log u2 + y = y2 - (y1 - 1)^(2::Int) / 2 + z = if u3 <= 0.5 then y1 else -y1 + normal' _ = error "normal_ar: infinite list of random variables expected" -- | Ratio Method (Kinderman-Monahan) (Knuth, v2, 2ed, pp 125--127) --- +-- -- If mu=0 and sigma=1, then this will generate numbers in the range -- [-1e15,1e15] (?) assuming that the uniform RNG is really giving full -- precision for doubles. normal_r :: (Double,Double) -- ^ (mu,sigma) - -> [Double] -- ^ U - -> [Double] -- ^ X + -> [Double] -- ^ U + -> [Double] -- ^ X -normal_r (mu,sigma) u = map (adjust (mu,sigma)) $ normal' u - where normal' (u:v:us) | x^2 <= -4 * log u = x : normal' us - | otherwise = normal' us - where x = a * (v - 0.5) / u - a = 1.71552776992141359295 -- sqrt $ 8 / e +normal_r muSigma = map (adjust muSigma) . normal' + where normal' (u:v:us) | x^(2::Int) <= -4 * log u = x : normal' us + | otherwise = normal' us + where x = a * (v - 0.5) / u + a = sqrt $ 8 / exp 1 -- 1.71552776992141359295 + normal' _ = error "normal_r: infinite list of random variables expected" diff --git a/Numeric/Random/Distribution/Poisson.hs b/Numeric/Random/Distribution/Poisson.hs --- a/Numeric/Random/Distribution/Poisson.hs +++ b/Numeric/Random/Distribution/Poisson.hs @@ -10,25 +10,76 @@ -- -- UNTESTED -- --- Module for transforming a list of uniform random variables into a --- list of Poisson random variables. +-- Module for transforming a list of uniform random variables +-- into a list of Poisson random variables. -- -- Reference: Ross +-- Donald E. Knuth (1969). Seminumerical Algorithms, The Art of Computer Programming, Volume 2 -- ---------------------------------------------------------------------------- -module Numeric.Random.Distribution.Poisson (poisson) where +module Numeric.Random.Distribution.Poisson (poisson, test, testHead) where +import Numeric.Statistics.Moment (mean) + +import Data.List (mapAccumL) +import System.Random (randomRs, mkStdGen) + + -- * Functions --- | Generates a list of poisson random variables from a list --- of uniforms +{- | +Generates a list of poisson random variables from a list of uniforms. +-} -poisson :: Double -- ^ lambda - -> [Double] -- ^ U - -> [Double] -- ^ X - -poisson lambda (u:us) = poisson' u us - where poisson' n (u:us) | n < e = n-1 : poisson lambda (u:us) - | otherwise = poisson' (n*u) us - e = exp (-lambda) +poisson :: Double -- ^ lambda - expectation value, should be non-negative. + -> [Double] -- ^ uniformly distributed values from the interval [0,1] + -> [Int] -- ^ Poisson distributed outputs + +poisson lambda (u:us) = + let e = exp (-lambda) + {- 'group' cannot replace segmentAfter here, + because it merges adjacent False values. -} + in map (length . tail) . segmentAfter not . snd $ + mapAccumL + (\p ui -> + let b = p >= e + in (if b then p*ui else ui, b)) + u us +poisson _ [] = + error "poisson: list of uniformly distributed values must not be empty" + + +{- | +Split after every element that satisfies the predicate. + +A candidate for a Utility module. +-} +segmentAfter :: (a -> Bool) -> [a] -> [[a]] +segmentAfter p = + foldr (\ x ~yt@(y:ys) -> if p x then [x]:yt else (x:y):ys) [[]] + + + +{- +The expectation value, +and thus the mean of a sequence of Poisson distributed values, +approximates lambda. +-} + +test :: Int -> Double -> Double +test n lambda = + mean $ map fromIntegral $ + take n $ poisson lambda $ + randomRs (0,1) $ mkStdGen 1 + +{- +Only test the leading number of several Poisson lists. +-} +testHead :: Int -> Double -> Double +testHead n lambda = + mean $ map fromIntegral $ + map + (head . poisson lambda . + randomRs (0,1) . mkStdGen) + [1..n] diff --git a/Numeric/Random/Distribution/Uniform.hs b/Numeric/Random/Distribution/Uniform.hs --- a/Numeric/Random/Distribution/Uniform.hs +++ b/Numeric/Random/Distribution/Uniform.hs @@ -25,7 +25,7 @@ -- 4294967295 = 2^32 - 1 uniform32cc :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform32cc xs = map ((/ 4294967295.0) . fromIntegral) $ xs @@ -34,21 +34,21 @@ -- 4294967296 = 2^32 uniform32co :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform32co xs = map ((/ 4294967296.0) . fromIntegral) $ xs -- | 32 bits in (0,1] uniform32oc :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform32oc xs = filter (/= 0) $ uniform32cc $ xs -- | 32 bits in (0,1) uniform32oo :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform32oo xs = filter (/= 1) $ uniform32oc $ xs @@ -58,12 +58,13 @@ -- 9007199254740991 = 2^53 - 1 uniform53cc :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U -uniform53cc xs = uniform' $ xs +uniform53cc xs = uniform' xs where uniform' (u1:u2:us) = (a * 67108864.0 + b) / 9007199254740991.0 : uniform' us - where a = fromIntegral u1 / 32.0 -- 27 bits - b = fromIntegral u2 / 64.0 -- 26 bits + where a = fromIntegral u1 / 32.0 -- 27 bits + b = fromIntegral u2 / 64.0 -- 26 bits + uniform' _ = error "uniform53cc: input list must be infinite" -- | 53 bits in [0,1), ie 64-bit IEEE 754 in [0,1) @@ -71,32 +72,33 @@ -- 9007199254740992 = 2^53 uniform53co :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform53co xs = uniform' $ xs where uniform' (u1:u2:us) = (a * 67108864.0 + b) / 9007199254740992.0 : uniform' us - where a = fromIntegral u1 / 32.0 -- 27 bits - b = fromIntegral u2 / 64.0 -- 26 bits + where a = fromIntegral u1 / 32.0 -- 27 bits + b = fromIntegral u2 / 64.0 -- 26 bits + uniform' _ = error "uniform53co: input list must be infinite" -- | 53 bits in (0,1] uniform53oc :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform53oc xs = filter (/= 0) $ uniform53cc $ xs -- | 53 bits in (0,1) uniform53oo :: [Word32] -- ^ X - -> [Double] -- ^ U + -> [Double] -- ^ U uniform53oo xs = filter (/= 1) $ uniform53oc $ xs -- | transforms uniform [0,1] to [a,b] uniform :: Double -- ^ a - -> Double -- ^ b - -> [Double] -- ^ U - -> [Double] -- ^ U' + -> Double -- ^ b + -> [Double] -- ^ U + -> [Double] -- ^ U' uniform a b us = map (\u -> (b-a)*u + a) us diff --git a/Numeric/Random/Generator/MT19937.hs b/Numeric/Random/Generator/MT19937.hs --- a/Numeric/Random/Generator/MT19937.hs +++ b/Numeric/Random/Generator/MT19937.hs @@ -19,12 +19,13 @@ -- but I can't get to the site anymore. As much as the orginal -- formatting has been retained as possible. --mpd +-- TODO: Make an instance of RandomGen class {- - Function genrand generates an infinite list of pseudorandom + Function genrand generates an infinite list of pseudorandom unsigned integers (32bit) which are uniformly distributed among 0 to 2^32-1. sgenrand(seed) uses an algorithm of Knuth - to provide 624 initial values to genrand(). + to provide 624 initial values to genrand(). Rewritten in Haskell by Matt Harden from original code in C by Takuji Nishimura. @@ -63,7 +64,7 @@ Vol. 8, No. 1, January 1998, pp 3--30. -} -module Numeric.Random.Generator.MT19937 (W, genrand) where +module Numeric.Random.Generator.MT19937 (W, genrand, test) where import Data.Word import Data.Bits @@ -77,47 +78,59 @@ type W = Word32 -- Period parameters -parm_N = 624 :: Int -parm_M = 397 :: Int -parm_A = 0x9908b0df :: W -uPPER_MASK = (bit 31) :: W -lOWER_MASK = (complement uPPER_MASK) :: W +parmN :: Int +parmN = 624 +parmM :: Int +parmM = 397 +parmA :: W +parmA = 0x9908b0df +upperMask :: W +upperMask = (bit 31) +lowerMask :: W +lowerMask = (complement upperMask) -- Tempering parameters -tEMPERING_MASK_B = (.&. 0x9d2c5680) :: W -> W -tEMPERING_MASK_C = (.&. 0xefc60000) :: W -> W -tEMPERING_SHIFT_U = (.>>. 11) :: W -> W -tEMPERING_SHIFT_S = (.<<. 7) :: W -> W -tEMPERING_SHIFT_T = (.<<. 15) :: W -> W -tEMPERING_SHIFT_L = (.>>. 18) :: W -> W +temperingMaskB :: W -> W +temperingMaskB = (.&. 0x9d2c5680) +temperingMaskC :: W -> W +temperingMaskC = (.&. 0xefc60000) +temperingShiftU :: W -> W +temperingShiftU = (.>>. 11) +temperingShiftS :: W -> W +temperingShiftS = (.<<. 7) +temperingShiftT :: W -> W +temperingShiftT = (.<<. 15) +temperingShiftL :: W -> W +temperingShiftL = (.>>. 18) -- A Knuth algorithm just to seed the seed... -- Line 25 of table 1 -- in [KNUTH 1981, The Art of Computer Programming Vol. 2 (2nd Ed.), pp102] sgenrand :: W -> [W] sgenrand 0 = sgenrand 4357 -- 0 not acceptable. Why 4357? I dunno. -sgenrand seed = take parm_N (iterate (69069 *) seed) +sgenrand seed = take parmN (iterate (69069 *) seed) -mag01 :: W -> W -mag01 0 = 0 -mag01 1 = parm_A +mag01 :: Bool -> W +mag01 False = 0 +mag01 True = parmA tempering :: W -> W tempering = let (^=) x f = xor x (f x) in - (^= (tEMPERING_SHIFT_L)) . - (^= (tEMPERING_MASK_C . tEMPERING_SHIFT_T)) . - (^= (tEMPERING_MASK_B . tEMPERING_SHIFT_S)) . - (^= (tEMPERING_SHIFT_U)) + (^= (temperingShiftL)) . + (^= (temperingMaskC . temperingShiftT)) . + (^= (temperingMaskB . temperingShiftS)) . + (^= (temperingShiftU)) -- parameter to rand MUST be a list of (_N) words! rand :: [W] -> [W] -rand init = map tempering r2 where - r = init ++ r2 - r2 = zipWith xor (map f r3) (drop parm_M r) - r3 = zipWith (\x y -> (x .&. uPPER_MASK) .|. (y .&. lOWER_MASK)) r (tail r) - f y = (y .>>. 1) `xor` (mag01 (y .&. 1)) - +rand seed = map tempering r2 where + r = seed ++ r2 + r2 = zipWith xor (map f r3) (drop parmM r) + r3 = zipWith (\x y -> (x .&. upperMask) .|. (y .&. lowerMask)) r (tail r) + f y = (y .>>. 1) `xor` (mag01 (odd y)) + genrand :: W -> [W] genrand = rand . sgenrand -test = sequence $ map print $ take 1000 $ genrand 4357 +test :: IO () +test = sequence_ $ map print $ take 1000 $ genrand 4357 diff --git a/Numeric/Random/Spectrum/Brown.hs b/Numeric/Random/Spectrum/Brown.hs --- a/Numeric/Random/Spectrum/Brown.hs +++ b/Numeric/Random/Spectrum/Brown.hs @@ -14,7 +14,7 @@ module Numeric.Random.Spectrum.Brown (brown) where -brown :: [Double] -- ^ noise +brown :: [Double] -- ^ noise -> [Double] -- ^ brown noise brown = scanl1 (+) diff --git a/Numeric/Random/Spectrum/Pink.hs b/Numeric/Random/Spectrum/Pink.hs --- a/Numeric/Random/Spectrum/Pink.hs +++ b/Numeric/Random/Spectrum/Pink.hs @@ -17,40 +17,43 @@ module Numeric.Random.Spectrum.Pink (kellet, voss) where +import DSP.Basic (downsample, upsampleAndHold) +import Data.List (tails) + ------------------------------------------------------------------------------- -- rb-j filter --- pole zero --- ---- ---- --- 0.99572754 0.98443604 --- 0.94790649 0.83392334 --- 0.53567505 0.07568359 +-- pole zero +-- ---- ---- +-- 0.99572754 0.98443604 +-- 0.94790649 0.83392334 +-- 0.53567505 0.07568359 ------------------------------------------------------------------------------- -- | Kellet's filter --- b0 = 0.99886 * b0 + white * 0.0555179; --- b1 = 0.99332 * b1 + white * 0.0750759; --- b2 = 0.96900 * b2 + white * 0.1538520; --- b3 = 0.86650 * b3 + white * 0.3104856; --- b4 = 0.55000 * b4 + white * 0.5329522; --- b5 = -0.7616 * b5 - white * 0.0168980; --- pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362; --- b6 = white * 0.115926; +-- b0 = 0.99886 * b0 + white * 0.0555179; +-- b1 = 0.99332 * b1 + white * 0.0750759; +-- b2 = 0.96900 * b2 + white * 0.1538520; +-- b3 = 0.86650 * b3 + white * 0.3104856; +-- b4 = 0.55000 * b4 + white * 0.5329522; +-- b5 = -0.7616 * b5 - white * 0.0168980; +-- pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362; +-- b6 = white * 0.115926; -kellet :: [Double] -- ^ noise +kellet :: [Double] -- ^ noise -> [Double] -- ^ pinked noise kellet w = kellet' w 0 0 0 0 0 0 0 where kellet' [] _ _ _ _ _ _ _ = [] kellet' (white:ws) b0 b1 b2 b3 b4 b5 b6 = pink : kellet' ws b0' b1' b2' b3' b4' b5' b6' - where b0' = 0.99886 * b0 + white * 0.0555179 - b1' = 0.99332 * b1 + white * 0.0750759 - b2' = 0.96900 * b2 + white * 0.1538520 - b3' = 0.86650 * b3 + white * 0.3104856 - b4' = 0.55000 * b4 + white * 0.5329522 + where b0' = 0.99886 * b0 + white * 0.0555179 + b1' = 0.99332 * b1 + white * 0.0750759 + b2' = 0.96900 * b2 + white * 0.1538520 + b3' = 0.86650 * b3 + white * 0.3104856 + b4' = 0.55000 * b4 + white * 0.5329522 b5' = -0.7616 * b5 - white * 0.0168980 pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362 b6' = white * 0.115926 @@ -60,33 +63,25 @@ -- voss algorithm add :: Num a => [[a]] -> [a] -add xs | any (== []) xs = [] - | otherwise = foldl1 (+) (map head xs) : add (map tail xs) - -hold :: Int -> [a] -> [a] -hold n xs = hold' n n xs - where hold' _ _ [] = [] - hold' n 1 (x:xs) = x : hold' n n xs - hold' n i (x:xs) = x : hold' n (i-1) (x:xs) +add = foldl1 (zipWith (+)) split :: Int -> [a] -> [[a]] -split n xs = split' n n xs - where split' _ 0 _ = [] - split' n i (x:xs) = split'' n (x:xs) : split' n (i-1) xs - split'' _ [] = [] - split'' n (x:xs) = x : split'' n (drop n (x:xs)) +split n = take n . map (downsample n) . (++repeat []) . tails mkOctaves :: [[a]] -> [[a]] -mkOctaves xss = mkOctaves' 1 xss +mkOctaves = mkOctaves' 1 where mkOctaves' _ [] = [] - mkOctaves' n (xs:xss) = hold n xs : mkOctaves' (2*n) xss + mkOctaves' n (xs:xss) = upsampleAndHold n xs : mkOctaves' (2*n) xss -- | Voss's algorithm -- -- UNTESTED, but the algorithm looks like it is working based on my hand -- tests. +-- TODO: Since the input noise is consumed in different speed for the octaves +-- the function will certainly leak memory. + voss :: Int -- ^ number of octaves to sum -> [Double] -- ^ noise -> [Double] -- ^ pinked noise @@ -99,4 +94,4 @@ ------------------------------------------------------------------------------- --- vm w = +-- vm w = diff --git a/Numeric/Random/Spectrum/Purple.hs b/Numeric/Random/Spectrum/Purple.hs --- a/Numeric/Random/Spectrum/Purple.hs +++ b/Numeric/Random/Spectrum/Purple.hs @@ -18,7 +18,7 @@ module Numeric.Random.Spectrum.Purple (purple) where -purple :: [Double] -- ^ noise +purple :: [Double] -- ^ noise -> [Double] -- ^ purple noise purple xs = zipWith (-) xs (0:xs) diff --git a/Numeric/Random/Spectrum/White.hs b/Numeric/Random/Spectrum/White.hs --- a/Numeric/Random/Spectrum/White.hs +++ b/Numeric/Random/Spectrum/White.hs @@ -16,7 +16,7 @@ module Numeric.Random.Spectrum.White (white) where -white :: [Double] -- ^ noise +white :: [Double] -- ^ noise -> [Double] -- ^ white noise white = id diff --git a/Numeric/Special/Trigonometric.hs b/Numeric/Special/Trigonometric.hs --- a/Numeric/Special/Trigonometric.hs +++ b/Numeric/Special/Trigonometric.hs @@ -1,6 +1,6 @@ -module Numeric.Special.Trigonometric (csc, sec, cot, +module Numeric.Special.Trigonometric (csc, sec, cot, acsc, asec, acot, - csch, sech, coth, + csch, sech, coth, acsch, asech, acoth ) where diff --git a/Numeric/Statistics/Covariance.hs b/Numeric/Statistics/Covariance.hs --- a/Numeric/Statistics/Covariance.hs +++ b/Numeric/Statistics/Covariance.hs @@ -29,5 +29,5 @@ where mu1 = mean x1 mu2 = mean x2 n = fromIntegral $ length $ x1 - f1 = \x -> (x - mu1)^2 - f2 = \x -> (x - mu2)^2 + f1 = \x -> (x - mu1)^(2::Int) + f2 = \x -> (x - mu2)^(2::Int) diff --git a/Numeric/Statistics/Median.hs b/Numeric/Statistics/Median.hs --- a/Numeric/Statistics/Median.hs +++ b/Numeric/Statistics/Median.hs @@ -14,13 +14,31 @@ -- ----------------------------------------------------------------------------- -module Numeric.Statistics.Median (median) where +module Numeric.Statistics.Median (median, medianFast) where import Data.List -- | Compute the median of a list median :: (Ord a, Fractional a) => [a] -> a -median x | odd n = sort x !! (n `div` 2) - | even n = ((sort x !! (n `div` 2 - 1)) + (sort x !! (n `div` 2))) / 2 +median x = + if odd n + then sort x !! (n `div` 2) + else ((sort x !! (n `div` 2 - 1)) + (sort x !! (n `div` 2))) / 2 where n = length x + + +{- | +Compute the center of the list in a more lazy manner +and thus halves memory requirement. +-} + +medianFast :: (Ord a, Fractional a) => [a] -> a +medianFast [] = error "medianFast: empty list has no median" +medianFast zs = + let recurse (x0:_) (_:[]) = x0 + recurse (x0:x1:_) (_:_:[]) = (x0+x1)/2 + recurse (_:xs) (_:_:ys) = recurse xs ys + recurse _ _ = + error "median: this error cannot occur in the way 'recurse' is called" + in recurse zs zs diff --git a/Numeric/Statistics/Moment.hs b/Numeric/Statistics/Moment.hs --- a/Numeric/Statistics/Moment.hs +++ b/Numeric/Statistics/Moment.hs @@ -49,7 +49,7 @@ -- where mu = mean x var :: (Fractional a) => [a] -> a -var xs = Prelude.sum (map (\x -> (x - mu)^2) xs) / (n - 1) +var xs = Prelude.sum (map (\x -> (x - mu)^(2::Int)) xs) / (n - 1) where mu = mean xs n = fromIntegral $ length $ xs @@ -74,7 +74,7 @@ -- @ Skew(X) = 1\/N sum(i=1..N) ((x_i-mu)\/sigma)^3 @ skew :: (RealFloat a) => [a] -> a -skew xs = Prelude.sum (map (\x -> ((x - mu) / sigma)^3) xs) / n +skew xs = Prelude.sum (map (\x -> ((x - mu) / sigma)^(3::Int)) xs) / n where mu = mean xs sigma = stddev xs n = fromIntegral $ length $ xs @@ -84,7 +84,7 @@ -- @ Kurt(X) = ( 1\/N sum(i=1..N) ((x_i-mu)\/sigma)^4 ) - 3@ kurtosis :: (RealFloat a) => [a] -> a -kurtosis xs = Prelude.sum (map (\x -> ((x - mu) / sigma)^4) xs) / n - 3 +kurtosis xs = Prelude.sum (map (\x -> ((x - mu) / sigma)^(4::Int)) xs) / n - 3 where mu = mean xs sigma = stddev xs n = fromIntegral $ length $ xs diff --git a/Numeric/Statistics/TTest.hs b/Numeric/Statistics/TTest.hs --- a/Numeric/Statistics/TTest.hs +++ b/Numeric/Statistics/TTest.hs @@ -31,8 +31,8 @@ where t = (mu1 - mu2) / s_d mu1 = Prelude.sum x1 / n1 mu2 = Prelude.sum x2 / n2 - v1 = Prelude.sum (map (\x -> (x - mu1)^2) x1) - v2 = Prelude.sum (map (\x -> (x - mu2)^2) x2) + v1 = Prelude.sum (map (\x -> (x - mu1)^(2::Int)) x1) + v2 = Prelude.sum (map (\x -> (x - mu2)^(2::Int)) x2) n1 = fromIntegral $ length $ x1 n2 = fromIntegral $ length $ x2 s_d = sqrt (((v1 + v2) / (n1+n2-2)) * (1/n1 + 1/n2)) @@ -62,5 +62,3 @@ var2 = var x2 s_d = sqrt ((var1 + var2 - 2 * cov x1 x2) / n) n = fromIntegral $ length $ x1 - - diff --git a/Numeric/Transform/Fourier/DFT.hs b/Numeric/Transform/Fourier/DFT.hs --- a/Numeric/Transform/Fourier/DFT.hs +++ b/Numeric/Transform/Fourier/DFT.hs @@ -47,8 +47,8 @@ {-# specialize dft' :: Array Int (Complex Double) -> Array Int (Complex Double) -> Int -> Array Int (Complex Double) #-} dft' :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Array a (Complex b) -dft' a w 1 = a +dft' a _ 1 = a dft' a w n = listArray (0,n-1) (sum [ a!k | k <- [0..(n-1)] ] : [ a!0 + sum [ a!k * wik i k | k <- [1..(n-1)] ] | i <- [1..(n-1)] ]) - where wik 0 k = 1 - wik i 0 = 1 + where wik 0 _ = 1 + wik _ 0 = 1 wik i k = w!(i*k `mod` n) 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 @@ -22,13 +22,16 @@ import Numeric.Transform.Fourier.FFTHard import Numeric.Transform.Fourier.R2DIF -import Numeric.Transform.Fourier.R2DIT +-- import Numeric.Transform.Fourier.R2DIT import Numeric.Transform.Fourier.R4DIF -import Numeric.Transform.Fourier.SRDIF +-- import Numeric.Transform.Fourier.SRDIF import Numeric.Transform.Fourier.CT import Numeric.Transform.Fourier.PFA import Numeric.Transform.Fourier.Rader +import DSP.Basic (uninterleave) + + ------------------------------------------------------------------------------- -- | This is the driver routine for calculating FFT's. All of the @@ -139,8 +142,7 @@ xa2 m = (xpr!m - cos w * xpi!m + sin w * xmr!m) :+ (xmi!m + sin w * xpi!m + cos w * xmr!m) where w = pi * fromIntegral m / fromIntegral n2 - rfft_unzip [] = [] - rfft_unzip (x1:x2:xs) = (x1:+x2) : rfft_unzip xs + rfft_unzip = uncurry (zipWith (:+)) . uninterleave n = (snd (bounds a) + 1) n2 = n `div` 2 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 @@ -12,9 +12,11 @@ -- ----------------------------------------------------------------------------- -module Numeric.Transform.Fourier.FFTUtils (fft_mag, fft_db, fft_phase, fft_grd, fft_info, - rfft_mag, rfft_db, rfft_phase, rfft_grd, rfft_info, - write_fft_info, write_rfft_info) where +module Numeric.Transform.Fourier.FFTUtils ( + fft_mag, fft_db, fft_phase, fft_grd, fft_info, + rfft_mag, rfft_db, rfft_phase, rfft_grd, rfft_info, + write_fft_info, write_rfft_info, + ) where import System.IO import Data.Array @@ -23,83 +25,116 @@ import Numeric.Transform.Fourier.FFT import DSP.Unwrap +magsq :: RealFloat a => Complex a -> a magsq (x:+y) = x*x + y*y +log10 :: Floating a => a -> a log10 0 = -1.0e9 log10 x = logBase 10 x +dot :: RealFloat a => Complex a -> Complex a -> a dot a b = realPart a * realPart b + imagPart a * imagPart b +eps :: Double eps = 1.0e-1 :: Double -- General functions +fft_mag :: (RealFloat b, Integral a, Ix a) => + Array a (Complex b) -> Array a b fft_mag x = fmap magnitude $ fft $ x +fft_db :: (RealFloat b, Integral a, Ix a) => + Array a (Complex b) -> Array a b fft_db x = fmap (10 *) $ fmap log10 $ fmap magsq $ fft $ x +fft_phase :: (Integral a, Ix a) => + Array a (Complex Double) -> Array a Double fft_phase x = unwrap eps $ fmap phase $ fft $ x +fft_grd :: (Integral i, RealFloat a, Ix i) => + Array i (Complex a) -> Array i a fft_grd x = listArray (bounds x') [ dot (x'!i) (dx'!i) / magsq (x'!i) | i <- indices x' ] where x' = fft x dx' = fft $ listArray (bounds x) [ fromIntegral i * x!i | i <- indices x ] +fft_info :: (Integral i, Ix i) => + Array i (Complex Double) + -> (Array i Double, Array i Double, Array i Double, Array i Double) fft_info x = (mag,db,arg,grd) where x' = fft x dx' = fft $ listArray (bounds x) [ fromIntegral i * x!i | i <- indices x ] mag = fmap magnitude $ x' - db = fmap (10 *) $ fmap log10 $ fmap magsq $ x' - arg = unwrap eps $ fmap phase $ x' - grd = listArray (bounds x') [ dot (x'!i) (dx'!i) / magsq (x'!i) | i <- indices x' ] + db = fmap (10 *) $ fmap log10 $ fmap magsq $ x' + arg = unwrap eps $ fmap phase $ x' + grd = listArray (bounds x') [ dot (x'!i) (dx'!i) / magsq (x'!i) | i <- indices x' ] +rfft_mag :: (RealFloat b, Integral a, Ix a) => + Array a b -> Array a b rfft_mag x = fmap magnitude $ rfft $ x +rfft_db :: (RealFloat b, Integral a, Ix a) => + Array a b -> Array a b rfft_db x = fmap (10 *) $ fmap log10 $ fmap magsq $ rfft $ x +rfft_phase :: (Integral a, Ix a) => + Array a Double -> Array a Double rfft_phase x = unwrap eps $ fmap phase $ rfft $ x +rfft_grd :: (Integral i, Ix i, RealFloat a) => + Array i a -> Array i a rfft_grd x = listArray (bounds x') [ dot (x'!i) (dx'!i) / magsq (x'!i) | i <- indices x' ] where x' = rfft x dx' = rfft $ listArray (bounds x) [ fromIntegral i * x!i | i <- indices x ] - dot a b = realPart a * realPart b + imagPart a * imagPart b -- I/O +rfft_info :: (Integral i, Ix i) => + Array i Double + -> (Array i Double, Array i Double, Array i Double, Array i Double) rfft_info x = (mag,db,arg,grd) where x' = rfft x dx' = rfft $ listArray (bounds x) [ fromIntegral i * x!i | i <- indices x ] mag = fmap magnitude $ x' - db = fmap (10 *) $ fmap log10 $ fmap magsq $ x' - arg = unwrap eps $ fmap phase $ x' - grd = listArray (bounds x') [ dot (x'!i) (dx'!i) / magsq (x'!i) | i <- indices x' ] + db = fmap (10 *) $ fmap log10 $ fmap magsq $ x' + arg = unwrap eps $ fmap phase $ x' + grd = listArray (bounds x') [ dot (x'!i) (dx'!i) / magsq (x'!i) | i <- indices x' ] +hPrintIndex :: (Integral a, Integral i, Show b) => + Handle -> i -> (a, b) -> IO () hPrintIndex h n (i,x) = do - hPutStr h $ show (fromIntegral i / fromIntegral n) - hPutStr h $ " " - hPutStrLn h $ show x + hPutStr h $ show (fromIntegral i / fromIntegral n :: Double) + hPutStr h $ " " + hPutStrLn h $ show x +write_cvector :: (Show e, Integral i, Ix i) => + 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 + let n = (snd $ bounds x) + 1 + h <- openFile f WriteMode + sequence $ map (hPrintIndex h n) $ assocs $ x + hClose h +write_fft_info :: (Ix i, Integral i) => + [Char] -> 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 - write_cvector (b ++ "_db.out") mag - write_cvector (b ++ "_arg.out") mag - write_cvector (b ++ "_grd.out") mag + let (mag,db,arg,grd) = fft_info x + write_cvector (b ++ "_mag.out") mag + write_cvector (b ++ "_db.out") db + write_cvector (b ++ "_arg.out") arg + write_cvector (b ++ "_grd.out") grd +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 + let n = (snd $ bounds x) + 1 + h <- openFile f WriteMode + sequence $ map (hPrintIndex h n) $ take (n `div` 2) $ assocs $ x + hClose h +write_rfft_info :: [Char] -> Array Int Double -> IO () write_rfft_info b x = do - let (mag,db,arg,grd) = rfft_info x - write_rvector (b ++ "_mag.out") mag - write_rvector (b ++ "_db.out") db - write_rvector (b ++ "_arg.out") arg - write_rvector (b ++ "_grd.out") grd + let (mag,db,arg,grd) = rfft_info x + write_rvector (b ++ "_mag.out") mag + write_rvector (b ++ "_db.out") db + write_rvector (b ++ "_arg.out") arg + write_rvector (b ++ "_grd.out") grd diff --git a/Numeric/Transform/Fourier/Goertzel.hs b/Numeric/Transform/Fourier/Goertzel.hs --- a/Numeric/Transform/Fourier/Goertzel.hs +++ b/Numeric/Transform/Fourier/Goertzel.hs @@ -8,7 +8,7 @@ -- Stability : experimental -- Portability : portable -- --- This is an implementation of Goertzel's algorithm, which computes on +-- This is an implementation of Goertzel's algorithm, which computes one -- bin of a DFT. A description can be found in Oppenheim and Schafer's -- /Discrete Time Signal Processing/, pp 585-587. -- @@ -30,12 +30,12 @@ -> b -- ^ k -> Complex a -- ^ X[k] -cgoertzel x k = g (elems x) 0 0 +cgoertzel x0 k = g (elems x0) 0 0 where w = 2 * pi * fromIntegral k / fromIntegral n a = 2 * cos w g [] x1 x2 = x1 * cis w - x2 g (x:xs) x1@(x1r:+x1i) x2 = g xs (x + (a*x1r:+a*x1i) - x2) x1 - n = (snd $ bounds x) - 1 + n = (snd $ bounds x0) - 1 -- | Power via Goertzel's algorithm for complex inputs @@ -43,7 +43,7 @@ -> b -- ^ k -> a -- ^ |X[k]|^2 -cgoertzel_power x k = (magnitude $ cgoertzel x k)^2 +cgoertzel_power x k = (magnitude $ cgoertzel x k)^(2::Int) -- | Goertzel's algorithm for real inputs @@ -51,12 +51,12 @@ -> b -- ^ k -> Complex a -- ^ X[k] -rgoertzel x k = g (elems x) 0 0 +rgoertzel x0 k = g (elems x0) 0 0 where w = 2 * pi * fromIntegral k / fromIntegral n a = 2 * cos w g [] x1 x2 = ((x1 - cos w * x2) :+ x2 * sin w) g (x:xs) x1 x2 = g xs (x + a * x1 - x2) x1 - n = (snd $ bounds x) - 1 + n = (snd $ bounds x0) - 1 -- | Power via Goertzel's algorithm for real inputs @@ -64,9 +64,9 @@ -> b -- ^ k -> a -- ^ |X[k]|^2 -rgoertzel_power x k = g (elems x) 0 0 +rgoertzel_power x0 k = g (elems x0) 0 0 where w = 2 * pi * fromIntegral k / fromIntegral n a = 2 * cos w - g [] x1 x2 = x1^2 + x2^2 - a * x1 * x2 + g [] x1 x2 = x1^(2::Int) + x2^(2::Int) - a * x1 * x2 g (x:xs) x1 x2 = g xs (x + a * x1 - x2) x1 - n = (snd $ bounds x) - 1 + n = (snd $ bounds x0) - 1 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 @@ -49,7 +49,7 @@ {-# specialize find_inverse :: Int -> Int -> Int #-} find_inverse :: (Integral a) => a -> a -> a -find_inverse a n = find_inverse' a n 1 +find_inverse a0 n0 = find_inverse' a0 n0 1 where find_inverse' a n a' | (a*a') `mod` n == 1 = a' | otherwise = find_inverse' a n (a'+1) 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 @@ -14,6 +14,7 @@ module Numeric.Transform.Fourier.R2DIF (fft_r2dif) where +import DSP.Basic (interleave) import Data.List import Data.Array import Data.Complex @@ -38,6 +39,4 @@ ye = fft ae yo = fft ao y = listArray (0,n-1) (interleave (elems ye) (elems yo)) - interleave [] [] = [] - interleave (e:es) (o:os) = e : o : interleave es os n2 = n `div` 2 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 @@ -14,6 +14,7 @@ module Numeric.Transform.Fourier.R4DIF (fft_r4dif) where +import DSP.Basic (interleave) import Data.List import Data.Array import Data.Complex @@ -43,8 +44,6 @@ j = 0 :+ 1 wn = cis (-2 * pi / fromIntegral n) w = listArray (0,n-1) $ iterate (* wn) 1 - interleave [] [] = [] - interleave (e:es) (o:os) = e : o : interleave es os n2 = n `div` 2 n4 = n `div` 4 n34 = 3 * n4 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 @@ -38,7 +38,7 @@ f' = array (0,n-1) ((0, sum [ f!i | i <- [0..(n-1)] ]) : [ (a ^* i, f!0 + hg!i) | i <- [0..(n-2)] ]) wn = cis (-2 * pi / fromIntegral n) w = listArray (0,n-1) $ iterate (* wn) 1 - i ^* 0 = 1 + _ ^* 0 = 1 i ^* j = (i * (i ^* (j-1))) `mod` n a = generator n @@ -62,10 +62,10 @@ f' = array (0,n-1) ((0, sum [ f!i | i <- [0..(n-1)] ]) : [ (a ^* i, f!0 + hg!i) | i <- [0..(n-2)] ]) wn = cis (-2 * pi / fromIntegral n) w = listArray (0,n-1) $ iterate (* wn) 1 - i ^* 0 = 1 + _ ^* 0 = 1 i ^* j = (i * (i ^* (j-1))) `mod` n a = generator n - ifft a = fmap (/ fromIntegral (n-1)) $ fmap swap $ fft $ fmap swap a + ifft b = fmap (/ fromIntegral (n-1)) $ fmap swap $ fft $ fmap swap b swap (x:+y) = (y:+x) -- Haskell translation of find_generator from FFTW @@ -74,8 +74,8 @@ generator :: (Integral a) => a -> a generator p = findgen 1 - where findgen 0 = error "rader: generator: no primative root?" + where findgen 0 = error "rader: generator: no primitive root?" findgen x | (period x x) == (p - 1) = x | otherwise = findgen ((x + 1) `mod` p) - period x 1 = 1 + period _ 1 = 1 period x prod = 1 + (period x (prod * x `mod` p)) 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 @@ -14,7 +14,8 @@ module Numeric.Transform.Fourier.SRDIF (fft_srdif) where -import Data.List +import DSP.Basic (interleave) +-- import Data.List import Data.Array import Data.Complex @@ -41,8 +42,6 @@ j = 0 :+ 1 wn = cis (-2 * pi / fromIntegral n) w = listArray (0,n-1) $ iterate (* wn) 1 - interleave [] [] = [] - interleave (e:es) (o:os) = e : o : interleave es os n2 = n `div` 2 n4 = n `div` 4 n34 = 3 * n4 diff --git a/Numeric/Transform/Fourier/SlidingFFT.hs b/Numeric/Transform/Fourier/SlidingFFT.hs --- a/Numeric/Transform/Fourier/SlidingFFT.hs +++ b/Numeric/Transform/Fourier/SlidingFFT.hs @@ -40,23 +40,42 @@ -> [Complex a] -- ^ x[n] -> [Array Int (Complex a)] -- ^ [X[k]] +sfft _ [] = error "sfft: input must have at least on value" sfft n (x:xs) = x' : sfft' n x xs x' where x' = fft $ listArray (0,n-1) $ reverse $ take n (x:xs) {-# specialize sfft' :: Int -> Complex Float -> [Complex Float] -> Array Int (Complex Float) -> [Array Int (Complex Float)] #-} {-# specialize sfft' :: Int -> Complex Double -> [Complex Double] -> Array Int (Complex Double) -> [Array Int (Complex Double)] #-} -sfft' :: RealFloat a => Int -> Complex a -> [Complex a] -> Array Int (Complex a) -> [Array Int (Complex a)] +sfft' :: RealFloat a => Int + -> Complex a + -> [Complex a] + -> Array Int (Complex a) + -> [Array Int (Complex a)] sfft' n xn (x:xs) x' | enough n (x:xs) = x'' : sfft' n x xs x'' | otherwise = [] where x'' = listArray (0,n-1) [ x0 - xn + x'!i * w i | i <- [0..(n-1)] ] x0 = xs !! (n-2) w i = cis $ -2 * pi * fromIntegral i / fromIntegral n +sfft' _ _ [] _ = error "sfft': input must have at least on value" -- We can't use Prelude.length because we may be operating on infinite, -- or ginormous lists. So enough will return True is there is enough -- data to perform the next FFT update, or False if there is not enough. +enough :: Int -> [a] -> Bool enough _ [] = False -enough 1 (x:_) = True -enough n (x:xs) = enough (n-1) xs +enough 1 (_:_) = True +enough n (_:xs) = enough (n-1) xs + +{- +Lemming: Me seems that the right implementation is + +enough 0 _ = True +enough _ [] = False +enough n (_:xs) = enough (n-1) xs + +or + +enough n xs = n<=0 || not (null (drop (n-1) xs)) +-} diff --git a/Polynomial/Basic.hs b/Polynomial/Basic.hs --- a/Polynomial/Basic.hs +++ b/Polynomial/Basic.hs @@ -29,21 +29,23 @@ -- | Evaluate a polynomial using Horner's method. polyeval :: Num a => [a] -> a -> a -polyeval [] x = 0 +polyeval [] _ = 0 polyeval (p:ps) x = p + x * polyeval ps x -- | Add two polynomials polyadd :: Num a => [a] -> [a] -> [a] -polyadd [] [] = [] polyadd [] ys = ys polyadd xs [] = xs polyadd (x:xs) (y:ys) = (x+y) : polyadd xs ys +polyAddScalar :: Num a => a -> [a] -> [a] +polyAddScalar c [] = [c] +polyAddScalar c (x:xs) = (c+x):xs + -- | Subtract two polynomials polysub :: Num a => [a] -> [a] -> [a] -polysub [] [] = [] polysub [] ys = map negate ys polysub xs [] = xs polysub (x:xs) (y:ys) = (x-y) : polysub xs ys @@ -56,58 +58,77 @@ -- | Multiply two polynomials polymult :: Num a => [a] -> [a] -> [a] -polymult (x:[]) ys = map (x*) ys -polymult (x:xs) ys = polyadd (map (x*) ys) (polymult xs (0:ys)) +polymult ys = + foldr (\x acc -> polyadd (polyscale x ys) (0 : acc)) [] +polymultAlt :: Num a => [a] -> [a] -> [a] +polymultAlt _ [] = [] +polymultAlt ys (x:xs) = polyadd (polyscale x ys) (0 : polymult ys xs) + -- | Divide two polynomials polydiv :: Fractional a => [a] -> [a] -> [a] -polydiv x y = reverse $ polydiv' (reverse x) (reverse y) - where polydiv' (x:xs) y | length (x:xs) < length y = [] - | otherwise = z : (polydiv' (tail (polysub (x:xs) (polymult [z] y))) y) - where z = x / head y +polydiv x0 y0 = reverse $ polydiv' (reverse x0) (reverse y0) + where polydiv' (x:xs) y + | length (x:xs) < length y = [] + | otherwise = z : (polydiv' (tail (polysub (x:xs) (polymult [z] y))) y) + where z = x / head y + polydiv' [] _ = [] -- | Modulus of two polynomials (remainder of division) polymod :: Fractional a => [a] -> [a] -> [a] -polymod x y = reverse $ polymod' (reverse x) (reverse y) - where polymod' (x:xs) y | length (x:xs) < length y = (x:xs) - | otherwise = polymod' (tail (polysub (x:xs) (polymult [z] y))) y - where z = x / head y +polymod x0 y0 = reverse $ polymod' (reverse x0) (reverse y0) + where polymod' (x:xs) y + | length (x:xs) < length y = (x:xs) + | otherwise = polymod' (tail (polysub (x:xs) (polymult [z] y))) y + where z = x / head y + polymod' [] _ = [] -- | Raise a polynomial to a non-negative integer power polypow :: (Num a, Integral b) => [a] -> b -> [a] -polypow x 0 = [ 1 ] +polypow _ 0 = [ 1 ] polypow x 1 = x -polypow x 2 = polymult x x -polypow x n | even n = polymult x2 x2 - | odd n = polymult x (polymult x2 x2) - where x2 = polypow x (n `div` 2) +polypow x n | even n = xSqr + | otherwise = polymult x xSqr + where xSqr = polymult x2 x2 + x2 = polypow x (n `div` 2) -- | Polynomial substitution y(n) = x(w(n)) polysubst :: Num a => [a] -> [a] -> [a] -polysubst w x = foldr polyadd [0] (polysubst' 0 w x ) +polysubst ws = foldr (\x -> polyAddScalar x . polymult ws) [] + +polysubstAlt :: Num a => [a] -> [a] -> [a] +polysubstAlt w0 x0 = foldr polyadd [0] (polysubst' 0 w0 x0) where polysubst' _ _ [] = [] - polysubst' n w (x:xs) = map (x*) (polypow w n) : polysubst' (n+1) w xs + polysubst' n w (x:xs) = polyscale x (polypow w (n::Int)) : polysubst' (n+1) w xs +{- | +Polynomial substitution @y(n) = x(w(n))@ +where the coefficients of @x@ are also polynomials. +-} +polyPolySubst :: Num a => [a] -> [[a]] -> [a] +polyPolySubst ws = foldr (\x -> polyadd x . polymult ws) [] + -- | Polynomial derivative polyderiv :: Num a => [a] -> [a] -polyderiv (x:xs) = polyderiv' 1 xs +polyderiv [] = [] +polyderiv (_:xs0) = polyderiv' 1 xs0 where polyderiv' _ [] = [] polyderiv' n (x:xs) = n * x : polyderiv' (n+1) xs -- | Polynomial integration polyinteg :: Fractional a => [a] -> a -> [a] -polyinteg x c = c : polyinteg' 1 x +polyinteg x0 c = c : polyinteg' 1 x0 where polyinteg' _ [] = [] polyinteg' n (x:xs) = x / n : polyinteg' (n+1) xs -- | Convert roots to a polynomial roots2poly :: Num a => [a] -> [a] -roots2poly (r:[]) = [-r, 1] +roots2poly [] = [1] roots2poly (r:rs) = polymult [-r, 1] (roots2poly rs) diff --git a/Polynomial/Maclaurin.hs b/Polynomial/Maclaurin.hs --- a/Polynomial/Maclaurin.hs +++ b/Polynomial/Maclaurin.hs @@ -26,8 +26,6 @@ polycos, polysin, polyatan, polycosh, polysinh, polyatanh) where -import Polynomial.Basic - -- A few utility lists ifacs :: [Double] diff --git a/Polynomial/Roots.hs b/Polynomial/Roots.hs --- a/Polynomial/Roots.hs +++ b/Polynomial/Roots.hs @@ -81,7 +81,7 @@ -> Int -- ^ iteration limit -> [Complex a] -- ^ the polynomial -> [Complex a] -- ^ the roots -roots eps count as = +roots eps0 count0 as0 = -- -- List of complex roots of a polynomial -- a0 + a1*x + a2*x^2... @@ -91,7 +91,7 @@ -- count is a maximum count of iterations allowed -- Require: list 'as' must have at least two elements -- and the last element must not be zero - roots' eps count as [] + roots' eps0 count0 as0 [] where roots' eps count as xs | length as <= 2 = x:xs @@ -100,14 +100,14 @@ where x = laguerre eps count as 0 bs = drop 1 (reverse (drop 1 as)) - deflate z bs cs - | bs == [] = cs + deflate z bs' cs + | bs' == [] = cs | otherwise = - deflate z (tail bs) (((head bs)+z*(head cs)):cs) + deflate z (tail bs') (((head bs')+z*(head cs)):cs) laguerre :: RealFloat a => a -> Int -> [Complex a] -> Complex a -> Complex a -laguerre eps count as x +laguerre eps0 count as0 x0 -- -- One of the roots of the polynomial 'as', -- where @@ -116,12 +116,12 @@ -- x is initial guess of the root -- This method is due to Laguerre. -- - | count <= 0 = x - | magnitude (x - x') < eps = x' - | otherwise = laguerre eps (count - 1) as x' - where x' = laguerre2 eps as as' as'' x - as' = polyderiv as - as'' = polyderiv as' + | count <= 0 = x0 + | magnitude (x0 - x0') < eps0 = x0' + | otherwise = laguerre eps0 (count - 1) as0 x0' + where x0' = laguerre2 eps0 as0 as0' as0'' x0 + as0' = polyderiv as0 + as0'' = polyderiv as0' laguerre2 eps as as' as'' x -- One iteration step | magnitude b < eps = x @@ -137,7 +137,7 @@ b = polyeval as x d = polyeval as' x f = polyeval as'' x - g2 = g^2 + g2 = g^(2::Int) n = fromIntegral (length as) -- Original Copyright Notice diff --git a/dsp.cabal b/dsp.cabal --- a/dsp.cabal +++ b/dsp.cabal @@ -1,18 +1,18 @@ Name: dsp -Version: 0.1 +Version: 0.2 License: GPL -Copyright: Matt Donadio, 2003 -Author: Matt Donadio +License-File: LICENSE +Copyright: Matthew Donadio, 2003 +Author: Matthew Donadio Maintainer: Henning Thielemann Stability: Experimental Homepage: http://haskelldsp.sourceforge.net/ -Synopsis: Haskell Digital Signal Processing -Description: Digital Signal Processing, Fourier Transform, Linear Algebra, Interpolation +Synopsis: Digital Signal Processing +Description: Digital Signal Processing, Fourier Transform, Filter design, Frequency estimation, Interpolation, Linear Algebra, Polynomials Category: Sound Tested-With: GHC Build-Depends: base -GHC-Options: -O2 --- -Wall +GHC-Options: -O2 -Wall Exposed-modules: DSP.Basic DSP.Convolution