dsp 0.1 → 0.2
raw patch · 71 files changed
+1956/−953 lines, 71 files
Files
- DSP/Basic.hs +67/−25
- DSP/Convolution.hs +11/−5
- DSP/Correlation.hs +25/−18
- DSP/Covariance.hs +15/−14
- DSP/Estimation/Frequency/FCI.hs +8/−6
- DSP/Estimation/Frequency/PerMax.hs +21/−18
- DSP/Estimation/Frequency/Pisarenko.hs +7/−1
- DSP/Estimation/Frequency/QuinnFernandes.hs +5/−1
- DSP/Estimation/Frequency/WLP.hs +31/−14
- DSP/Estimation/Spectral/AR.hs +6/−5
- DSP/Estimation/Spectral/ARMA.hs +7/−3
- DSP/Estimation/Spectral/KayData.hs +68/−66
- DSP/Estimation/Spectral/MA.hs +0/−1
- DSP/Filter/Analog/Prototype.hs +7/−5
- DSP/Filter/Analog/Response.hs +5/−4
- DSP/Filter/Analog/Transform.hs +59/−38
- DSP/Filter/FIR/FIR.hs +31/−24
- DSP/Filter/FIR/Kaiser.hs +10/−3
- DSP/Filter/FIR/PolyInterp.hs +144/−144
- DSP/Filter/FIR/Sharpen.hs +4/−5
- DSP/Filter/FIR/Smooth.hs +8/−4
- DSP/Filter/FIR/Taps.hs +8/−2
- DSP/Filter/FIR/Window.hs +66/−42
- DSP/Filter/IIR/Bilinear.hs +17/−11
- DSP/Filter/IIR/Cookbook.lhs +12/−10
- DSP/Filter/IIR/Design.hs +20/−17
- DSP/Filter/IIR/IIR.hs +30/−18
- DSP/Filter/IIR/Prony.hs +3/−3
- DSP/Filter/IIR/Transform.hs +32/−54
- DSP/Multirate/CIC.hs +4/−4
- DSP/Multirate/Halfband.hs +7/−19
- DSP/Multirate/Polyphase.hs +3/−4
- DSP/Source/Basic.hs +3/−3
- DSP/Source/Oscillator.hs +6/−5
- DSP/Unwrap.hs +1/−1
- LICENSE +674/−0
- Matrix/Cholesky.hs +1/−1
- Matrix/LU.hs +38/−36
- Matrix/Levinson.hs +2/−2
- Matrix/Matrix.hs +3/−3
- Matrix/Simplex.hs +73/−56
- Numeric/Random/Distribution/Binomial.hs +1/−1
- Numeric/Random/Distribution/Gamma.hs +1/−1
- Numeric/Random/Distribution/Geometric.hs +1/−1
- Numeric/Random/Distribution/Normal.hs +44/−37
- Numeric/Random/Distribution/Poisson.hs +64/−13
- Numeric/Random/Distribution/Uniform.hs +18/−16
- Numeric/Random/Generator/MT19937.hs +42/−29
- Numeric/Random/Spectrum/Brown.hs +1/−1
- Numeric/Random/Spectrum/Pink.hs +30/−35
- Numeric/Random/Spectrum/Purple.hs +1/−1
- Numeric/Random/Spectrum/White.hs +1/−1
- Numeric/Special/Trigonometric.hs +2/−2
- Numeric/Statistics/Covariance.hs +2/−2
- Numeric/Statistics/Median.hs +21/−3
- Numeric/Statistics/Moment.hs +3/−3
- Numeric/Statistics/TTest.hs +2/−4
- Numeric/Transform/Fourier/DFT.hs +3/−3
- Numeric/Transform/Fourier/FFT.hs +6/−4
- Numeric/Transform/Fourier/FFTUtils.hs +66/−31
- Numeric/Transform/Fourier/Goertzel.hs +9/−9
- Numeric/Transform/Fourier/PFA.hs +1/−1
- Numeric/Transform/Fourier/R2DIF.hs +1/−2
- Numeric/Transform/Fourier/R4DIF.hs +1/−2
- Numeric/Transform/Fourier/Rader.hs +5/−5
- Numeric/Transform/Fourier/SRDIF.hs +2/−3
- Numeric/Transform/Fourier/SlidingFFT.hs +22/−3
- Polynomial/Basic.hs +44/−23
- Polynomial/Maclaurin.hs +0/−2
- Polynomial/Roots.hs +13/−13
- dsp.cabal +7/−7
DSP/Basic.hs view
@@ -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
DSP/Convolution.hs view
@@ -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
DSP/Correlation.hs view
@@ -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) ]
DSP/Covariance.hs view
@@ -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
DSP/Estimation/Frequency/FCI.hs view
@@ -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
DSP/Estimation/Frequency/PerMax.hs view
@@ -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
DSP/Estimation/Frequency/Pisarenko.hs view
@@ -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
DSP/Estimation/Frequency/QuinnFernandes.hs view
@@ -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
DSP/Estimation/Frequency/WLP.hs view
@@ -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
DSP/Estimation/Spectral/AR.hs view
@@ -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
DSP/Estimation/Spectral/ARMA.hs view
@@ -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] ]
DSP/Estimation/Spectral/KayData.hs view
@@ -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) ]
DSP/Estimation/Spectral/MA.hs view
@@ -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
DSP/Filter/Analog/Prototype.hs view
@@ -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
DSP/Filter/Analog/Response.hs view
@@ -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)
DSP/Filter/Analog/Transform.hs view
@@ -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
DSP/Filter/FIR/FIR.hs view
@@ -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' ]
DSP/Filter/FIR/Kaiser.hs view
@@ -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
DSP/Filter/FIR/PolyInterp.hs view
@@ -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
DSP/Filter/FIR/Sharpen.hs view
@@ -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
DSP/Filter/FIR/Smooth.hs view
@@ -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
DSP/Filter/FIR/Taps.hs view
@@ -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
DSP/Filter/FIR/Window.hs view
@@ -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) ..]
DSP/Filter/IIR/Bilinear.hs view
@@ -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
DSP/Filter/IIR/Cookbook.lhs view
@@ -27,6 +27,8 @@ > import DSP.Filter.IIR.IIR +> import DSP.Basic((^!))+ Cookbook formulae for audio EQ biquad filter coefficients ----------------------------------------------------------------------------- by Robert Bristow-Johnson <robert@wavemechanics.com>@@ -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.)
DSP/Filter/IIR/Design.hs view
@@ -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')
DSP/Filter/IIR/IIR.hs view
@@ -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
DSP/Filter/IIR/Prony.hs view
@@ -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)
DSP/Filter/IIR/Transform.hs view
@@ -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) {-
DSP/Multirate/CIC.hs view
@@ -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) {-
DSP/Multirate/Halfband.hs view
@@ -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
DSP/Multirate/Polyphase.hs view
@@ -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 {-
DSP/Source/Basic.hs view
@@ -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
DSP/Source/Oscillator.hs view
@@ -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
DSP/Unwrap.hs view
@@ -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
+ LICENSE view
@@ -0,0 +1,674 @@+ GNU GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ 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.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ 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 <http://www.gnu.org/licenses/>.++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:++ <program> Copyright (C) <year> <name of author>+ 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+<http://www.gnu.org/licenses/>.++ 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+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
Matrix/Cholesky.hs view
@@ -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
Matrix/LU.hs view
@@ -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'' -}
Matrix/Levinson.hs view
@@ -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)
Matrix/Matrix.hs view
@@ -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] ]
Matrix/Simplex.hs view
@@ -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
Numeric/Random/Distribution/Binomial.hs view
@@ -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)
Numeric/Random/Distribution/Gamma.hs view
@@ -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
Numeric/Random/Distribution/Geometric.hs view
@@ -29,6 +29,6 @@ geometric :: Double -- ^ p -> [Double] -- ^ U -> [Double] -- ^ X- + geometric p us = map (\u -> 1 + log u / log (1 - p)) us
Numeric/Random/Distribution/Normal.hs view
@@ -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"
Numeric/Random/Distribution/Poisson.hs view
@@ -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]
Numeric/Random/Distribution/Uniform.hs view
@@ -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
Numeric/Random/Generator/MT19937.hs view
@@ -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
Numeric/Random/Spectrum/Brown.hs view
@@ -14,7 +14,7 @@ module Numeric.Random.Spectrum.Brown (brown) where -brown :: [Double] -- ^ noise +brown :: [Double] -- ^ noise -> [Double] -- ^ brown noise brown = scanl1 (+)
Numeric/Random/Spectrum/Pink.hs view
@@ -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 =
Numeric/Random/Spectrum/Purple.hs view
@@ -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)
Numeric/Random/Spectrum/White.hs view
@@ -16,7 +16,7 @@ module Numeric.Random.Spectrum.White (white) where -white :: [Double] -- ^ noise +white :: [Double] -- ^ noise -> [Double] -- ^ white noise white = id
Numeric/Special/Trigonometric.hs view
@@ -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
Numeric/Statistics/Covariance.hs view
@@ -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)
Numeric/Statistics/Median.hs view
@@ -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
Numeric/Statistics/Moment.hs view
@@ -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
Numeric/Statistics/TTest.hs view
@@ -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--
Numeric/Transform/Fourier/DFT.hs view
@@ -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)
Numeric/Transform/Fourier/FFT.hs view
@@ -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
Numeric/Transform/Fourier/FFTUtils.hs view
@@ -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
Numeric/Transform/Fourier/Goertzel.hs view
@@ -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
Numeric/Transform/Fourier/PFA.hs view
@@ -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)
Numeric/Transform/Fourier/R2DIF.hs view
@@ -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
Numeric/Transform/Fourier/R4DIF.hs view
@@ -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
Numeric/Transform/Fourier/Rader.hs view
@@ -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))
Numeric/Transform/Fourier/SRDIF.hs view
@@ -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
Numeric/Transform/Fourier/SlidingFFT.hs view
@@ -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))+-}
Polynomial/Basic.hs view
@@ -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)
Polynomial/Maclaurin.hs view
@@ -26,8 +26,6 @@ polycos, polysin, polyatan, polycosh, polysinh, polyatanh) where -import Polynomial.Basic- -- A few utility lists ifacs :: [Double]
Polynomial/Roots.hs view
@@ -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
dsp.cabal view
@@ -1,18 +1,18 @@ Name: dsp-Version: 0.1+Version: 0.2 License: GPL-Copyright: Matt Donadio, 2003-Author: Matt Donadio <m.p.donadio@ieee.org>+License-File: LICENSE+Copyright: Matthew Donadio, 2003+Author: Matthew Donadio <m.p.donadio@ieee.org> Maintainer: Henning Thielemann <haskell@henning-thielemann.de> 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