dsp 0.2.2 → 0.2.3
raw patch · 15 files changed
+389/−215 lines, 15 filesdep +dspdep ~arraydep ~basenew-component:exe:dsp-demo-articlenew-component:exe:dsp-demo-fft-benchnew-component:exe:dsp-demo-fft-testnew-component:exe:dsp-demo-freqnew-component:exe:dsp-demo-iirnew-component:exe:dsp-demo-noise
Dependencies added: dsp
Dependency ranges changed: array, base
Files
- DSP/Filter/IIR/Design.hs +77/−32
- Numeric/Statistics/Moment.hs +2/−2
- Numeric/Transform/Fourier/CT.hs +2/−2
- Numeric/Transform/Fourier/FFT.hs +2/−2
- Numeric/Transform/Fourier/FFTHard.hs +4/−4
- Numeric/Transform/Fourier/FFTUtils.hs +2/−2
- Numeric/Transform/Fourier/PFA.hs +2/−2
- Polynomial/Roots.hs +9/−9
- demo/Article.hs +16/−10
- demo/FFTBench.hs +34/−22
- demo/FFTTest.hs +42/−18
- demo/FreqDemo.hs +50/−38
- demo/IIRDemo.hs +29/−16
- demo/NoiseDemo.hs +38/−36
- dsp.cabal +80/−20
DSP/Filter/IIR/Design.hs view
@@ -8,7 +8,7 @@ -- Stability : experimental -- Portability : portable ----- Lowpass IIR design functions+-- Lowpass, Highpass, Bandpass IIR design functions -- -- Method: --@@ -20,48 +20,85 @@ -- ----------------------------------------------------------------------------- -module DSP.Filter.IIR.Design where+module DSP.Filter.IIR.Design (+ poly2iir,+ butterworthLowpass, butterworthHighpass, butterworthBandpass,+ chebyshev1Lowpass, chebyshev2Lowpass,+ mkButterworth, mkChebyshev1, mkChebyshev2,+ ) where -import DSP.Filter.Analog.Prototype-import DSP.Filter.Analog.Transform-import DSP.Filter.IIR.Bilinear+import qualified DSP.Filter.Analog.Prototype as Analog+import DSP.Filter.Analog.Transform (a_lp2lp, a_lp2hp, a_lp2bp)+import DSP.Filter.IIR.Bilinear (bilinear, prewarp) import DSP.Basic ((^!)) -import Data.Array+import Data.Array (Array, listArray) + 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- m = length b - 1- n = length a - 1+poly2iir (a,b) =+ let toArray x = listArray (0, length x - 1) $ reverse x+ in (toArray a, toArray b) + -- | Generates lowpass Butterworth IIR filters -mkButterworth :: (Double, Double) -- ^ (wp,dp)- -> (Double, Double) -- ^ (ws,ds)- -> (Array Int Double, Array Int Double) -- ^ (b,a)+butterworthLowpass, mkButterworth ::+ (Double, Double) -- ^ (wp,dp)+ -> (Double, Double) -- ^ (ws,ds)+ -> (Array Int Double, Array Int Double) -- ^ (b,a)+butterworthLowpass p s =+ let (n, s') = butterworthParams p s+ in butterworth (a_lp2lp $ wc n s') n -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)- wp' = prewarp wp 1- ws' = prewarp ws 1+butterworthHighpass, butterworthBandpass ::+ (Double, Double) ->+ (Double, Double) ->+ (Array Int Double, Array Int Double)+butterworthHighpass p s =+ let (n, s') = butterworthParams p s+ in butterworth (a_lp2hp $ wc n s') n +butterworthBandpass p@(wp, _dp) s@(ws, _ds) =+ let (n, _s') = butterworthParams p s+ in butterworth (a_lp2bp wp ws) n++butterworth ::+ (([Double], [Double]) -> ([Double], [Double])) ->+ Int -> (Array Int Double, Array Int Double)+butterworth analogToAnalog n =+ poly2iir $ bilinear 1 $ analogToAnalog $ Analog.butterworth n++butterworthParams ::+ (Double, Double) ->+ (Double, Double) ->+ (Int, (Double, Double))+butterworthParams (wp, dp) (ws, ds) =+ let n = ceiling $ log (((1/ds)^!2-1) / ((1/(1-dp))^!2-1)) / 2 / log (ws' / wp')+ wp' = prewarp wp 1+ ws' = prewarp ws 1+ in (n, (ws', ds))++wc :: Floating a => Int -> (a, a) -> a+wc n (ws', ds) = ws' / ((1/ds)^!2 - 1) ** (1/2/fromIntegral n)+++{-# DEPRECATED mkButterworth "Use butterworthLowpass instead" #-}+mkButterworth = butterworthLowpass++ -- | Generates lowpass Chebyshev IIR filters -mkChebyshev1 :: (Double, Double) -- ^ (wp,dp)- -> (Double, Double) -- ^ (ws,ds)- -> (Array Int Double, Array Int Double) -- ^ (b,a)+chebyshev1Lowpass, mkChebyshev1 ::+ (Double, Double) -- ^ (wp,dp)+ -> (Double, Double) -- ^ (ws,ds)+ -> (Array Int Double, Array Int Double) -- ^ (b,a) -mkChebyshev1 (wp,dp) (ws,ds) = poly2iir $+chebyshev1Lowpass (wp,dp) (ws,ds) = poly2iir $ bilinear 1 $ a_lp2lp wp' $- chebyshev1 eps n+ Analog.chebyshev1 eps n where wp' = prewarp wp 1 ws' = prewarp ws 1 eps = sqrt ((2 - dp)*dp) / (1 - dp)@@ -70,18 +107,26 @@ k = wp' / ws' n = ceiling $ acosh (1/k1) / log ((1 + sqrt (1 - k^!2)) / k) +{-# DEPRECATED mkChebyshev1 "Use chebyshev1Lowpass instead" #-}+mkChebyshev1 = chebyshev1Lowpass++ -- | Generates lowpass Inverse Chebyshev IIR filters -mkChebyshev2 :: (Double, Double) -- ^ (wp,dp)- -> (Double, Double) -- ^ (ws,ds)- -> (Array Int Double, Array Int Double) -- ^ (b,a)+chebyshev2Lowpass, mkChebyshev2 ::+ (Double, Double) -- ^ (wp,dp)+ -> (Double, Double) -- ^ (ws,ds)+ -> (Array Int Double, Array Int Double) -- ^ (b,a) -mkChebyshev2 (wp,dp) (ws,ds) = poly2iir $+chebyshev2Lowpass (wp,dp) (ws,ds) = poly2iir $ bilinear 1 $ a_lp2lp ws' $- chebyshev2 eps n+ Analog.chebyshev2 eps n where wp' = prewarp wp 1 ws' = prewarp ws 1 eps = ds / sqrt (1 - ds^!2) g = 1 - dp n = ceiling $ acosh (g / eps / sqrt (1 - g^!2)) / acosh (ws' / wp')++{-# DEPRECATED mkChebyshev2 "Use chebyshev2Lowpass instead" #-}+mkChebyshev2 = chebyshev2Lowpass
Numeric/Statistics/Moment.hs view
@@ -14,8 +14,8 @@ -- ----------------------------------------------------------------------------- -module Numeric.Statistics.Moment (mean, var, - stddev, avgdev, +module Numeric.Statistics.Moment (mean, var,+ stddev, avgdev, skew, kurtosis) where -- TODO: does mean pass though the list twice? once to compute the sum,
Numeric/Transform/Fourier/CT.hs view
@@ -81,14 +81,14 @@ {-# specialize rows :: Array (Int,Int) (Complex Float) -> [Array Int (Complex Float)] #-} {-# specialize rows :: Array (Int,Int) (Complex Double) -> [Array Int (Complex Double)] #-} -rows :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] +rows :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] rows x = [ listArray (0,m) [ x!(i,j) | j <- [0..m] ] | i <- [0..l] ] where ((_,_),(l,m)) = bounds x {-# specialize cols :: Array (Int,Int) (Complex Float) -> [Array Int (Complex Float)] #-} {-# specialize cols :: Array (Int,Int) (Complex Double) -> [Array Int (Complex Double)] #-} -cols :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] +cols :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] cols x = [ listArray (0,l) [ x!(i,j) | i <- [0..l] ] | j <- [0..m] ] where ((_,_),(l,m)) = bounds x
Numeric/Transform/Fourier/FFT.hs view
@@ -135,10 +135,10 @@ xmi = listArray (0,n2-1) (0 : [ (xi!m - xi!(n2-m)) / 2 | m <- [1..(n2-1)] ]) xr = fmap realPart x xi = fmap imagPart x- xa1 m = (xpr!m + cos w * xpi!m - sin w * xmr!m) :+ + xa1 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- xa2 m = (xpr!m - cos w * xpi!m + sin w * xmr!m) :+ + 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 = uncurry (zipWith (:+)) . uninterleave
Numeric/Transform/Fourier/FFTHard.hs view
@@ -27,7 +27,7 @@ fft'2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n] -> Array a (Complex b) -- ^ X[k] -fft'2 a = array (0,1) [ (0, ((tmp1 + tmp2) :+ (tmp3 + tmp4))), +fft'2 a = array (0,1) [ (0, ((tmp1 + tmp2) :+ (tmp3 + tmp4))), (1, ((tmp1 - tmp2) :+ (tmp3 - tmp4) )) ] where tmp1 = realPart (a!0) tmp3 = imagPart (a!0)@@ -68,9 +68,9 @@ fft'4 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -- ^ x[n] -> Array a (Complex b) -- ^ X[k] -fft'4 a = array (0,3) [ (0, (tmp3 + tmp6) :+ (tmp15 + tmp16)), - (1, (tmp11 + tmp14) :+ (tmp9 - tmp10)), - (2, (tmp3 - tmp6) :+ (tmp15 - tmp16)), +fft'4 a = array (0,3) [ (0, (tmp3 + tmp6) :+ (tmp15 + tmp16)),+ (1, (tmp11 + tmp14) :+ (tmp9 - tmp10)),+ (2, (tmp3 - tmp6) :+ (tmp15 - tmp16)), (3, (tmp11 - tmp14) :+ (tmp10 + tmp9)) ] where tmp1 = realPart (a!0) tmp7 = imagPart (a!0)
Numeric/Transform/Fourier/FFTUtils.hs view
@@ -61,7 +61,7 @@ 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) +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'@@ -92,7 +92,7 @@ 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) +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'
Numeric/Transform/Fourier/PFA.hs view
@@ -56,14 +56,14 @@ {-# specialize rows :: Array (Int,Int) (Complex Float) -> [Array Int (Complex Float)] #-} {-# specialize rows :: Array (Int,Int) (Complex Double) -> [Array Int (Complex Double)] #-} -rows :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] +rows :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] rows x = [ listArray (0,m) [ x!(i,j) | j <- [0..m] ] | i <- [0..l] ] where ((_,_),(l,m)) = bounds x {-# specialize cols :: Array (Int,Int) (Complex Float) -> [Array Int (Complex Float)] #-} {-# specialize cols :: Array (Int,Int) (Complex Double) -> [Array Int (Complex Double)] #-} -cols :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] +cols :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -> [Array a (Complex b)] cols x = [ listArray (0,l) [ x!(i,j) | i <- [0..l] ] | j <- [0..m] ] where ((_,_),(l,m)) = bounds x
Polynomial/Roots.hs view
@@ -59,7 +59,7 @@ -- 300 is a count of permitted iterations -- (You set it to small number just in case you -- do not trust the algorithm. But if you do,- -- then set it to something big, say 300) + -- then set it to something big, say 300) The answer is [2.0 :+ 0.0, 1.0 :+ 0.0]; that is, both roots are real and equal to 2 and 1:@@ -67,7 +67,7 @@ x^2 - 3 x + 2 = (x - 2) (x - 1) = 0 -} -module Polynomial.Roots (roots) where +module Polynomial.Roots (roots) where import Data.Complex @@ -90,19 +90,19 @@ -- eps is a desired accuracy -- 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 + -- and the last element must not be zero roots' eps0 count0 as0 [] where- roots' eps count as xs + roots' eps count as xs | length as <= 2 = x:xs- | otherwise = + | otherwise = roots' eps count (deflate x bs [last as]) (x:xs) where x = laguerre eps count as 0 bs = drop 1 (reverse (drop 1 as)) deflate z bs' cs | bs' == [] = cs- | otherwise = + | otherwise = deflate z (tail bs') (((head bs')+z*(head cs)):cs) @@ -125,9 +125,9 @@ laguerre2 eps as as' as'' x -- One iteration step | magnitude b < eps = x- | magnitude gp < magnitude gm = + | magnitude gp < magnitude gm = if gm == 0 then x - 1 else x - n/gm- | otherwise = + | otherwise = if gp == 0 then x - 1 else x - n/gp where gp = g + delta gm = g - delta@@ -155,5 +155,5 @@ -- License: -- -- GNU General Public License, GPL--- +-- -----------------------------------------------------------------------------
demo/Article.hs view
@@ -1,18 +1,19 @@ -- This program was used to generate the data for ----- Matthew Donadio, "Lost Knowledge Refound: Sharpened FIR Filters," +-- Matthew Donadio, "Lost Knowledge Refound: Sharpened FIR Filters," -- IEEE Signal Processing Magazine, to appear -module Main where+module Main (main) where -import Data.Array+import DSP.Filter.FIR.Sharpen (sharpen)+import DSP.Filter.FIR.FIR (fir)+import DSP.Source.Basic (impulse) -import DSP.Filter.FIR.FIR-import DSP.Filter.FIR.Sharpen-import DSP.Source.Basic+import Numeric.Transform.Fourier.FFTUtils (write_rfft_info) -import Numeric.Transform.Fourier.FFTUtils+import Data.Array (Array, listArray) + n :: Int n = 1000 @@ -22,11 +23,16 @@ 0.021409, -0.090271, -0.013137, 0.047422, 0.015799, -0.022174, -0.016674 ] +y1, y2, y3 :: [Double] y1 = fir h $ impulse y2 = fir h $ fir h $ impulse y3 = sharpen h $ impulse +example :: String -> [Double] -> IO ()+example name y = write_rfft_info name $ listArray (0,n-1) $ y++main :: IO () main = do- write_rfft_info "y1" $ listArray (0,999) $ y1- write_rfft_info "y2" $ listArray (0,999) $ y2- write_rfft_info "y3" $ listArray (0,999) $ y3+ example "y1" y1+ example "y2" y2+ example "y3" y3
demo/FFTBench.hs view
@@ -1,38 +1,47 @@-module Main where+module Main (main) where -import Data.Array-import Data.Complex+import Numeric.Transform.Fourier.FFT (fft)+import Numeric.Transform.Fourier.FFTHard (fft'2, fft'4)+import Numeric.Transform.Fourier.R2DIF (fft_r2dif)+import Numeric.Transform.Fourier.R2DIT (fft_r2dit)+import Numeric.Transform.Fourier.R4DIF (fft_r4dif)+import Numeric.Transform.Fourier.SRDIF (fft_srdif)+import Numeric.Transform.Fourier.CT (fft_ct1, fft_ct2)+import Numeric.Transform.Fourier.Rader (fft_rader1, fft_rader2) -import Numeric.Transform.Fourier.FFT-import Numeric.Transform.Fourier.FFTHard-import Numeric.Transform.Fourier.R2DIF-import Numeric.Transform.Fourier.R2DIT-import Numeric.Transform.Fourier.R4DIF-import Numeric.Transform.Fourier.SRDIF-import Numeric.Transform.Fourier.CT-import Numeric.Transform.Fourier.PFA-import Numeric.Transform.Fourier.Rader-import Numeric.Transform.Fourier.DFT+import Numeric.Random.Generator.MT19937 (genrand)+import Numeric.Random.Distribution.Uniform (uniform53cc) -import Numeric.Random.Generator.MT19937-import Numeric.Random.Distribution.Uniform+import Control.Monad (when)+import Data.Complex (Complex((:+)), magnitude)+import Data.Array (Array, listArray, elems, bounds) -len = 2048 :: Int-iter = 100 :: Int +len, iter :: Int+len = 2048+iter = 100++m1 :: Double -> Double m1 x = x - 1 +real, imag :: [Double] real = map m1 $ map (2*) $ uniform53cc $ genrand 42 imag = map m1 $ map (2*) $ uniform53cc $ genrand 43 -x = zipWith (:+) real imag+xl :: [Complex Double]+xl = zipWith (:+) real imag gendata :: [Complex Double] -> Int -> [Array Int (Complex Double)]-gendata xs n = map (listArray (0,n-1)) $ gendata' xs n- where gendata' xs n = take n xs : gendata' (drop n xs) n+gendata xs n = map (listArray (0,n-1)) $ slice xs n -calc f xs iter = magnitude $ sum $ map sum $ map elems $ map f $ take iter xs+slice :: [a] -> Int -> [[a]]+slice xs n = take n xs : slice (drop n xs) n +calc :: (array -> Array Int (Complex Double)) -> [array] -> Int -> Double+calc f xs n = magnitude $ sum $ map (sum . elems . f) $ take n xs++f1, f2, f3, f4, f5, f6, f7, f8 ::+ Array Int (Complex Double) -> Array Int (Complex Double) f1 xs | n == 2 = fft'2 xs | n == 4 = fft'4 xs | otherwise = fft_r2dit xs n f1@@ -90,11 +99,14 @@ f8 xs = fft_rader2 xs n fft where n = (snd $ bounds xs) + 1 +main :: IO () main = do- let xs = (gendata x len)+ let xs = gendata xl len print $ calc f1 xs iter print $ calc f2 xs iter print $ calc f3 xs iter print $ calc f4 xs iter print $ calc f5 xs iter print $ calc f6 xs iter+ when False $ print $ calc f7 xs iter+ when False $ print $ calc f8 xs iter
demo/FFTTest.hs view
@@ -8,17 +8,21 @@ -- generator bottleneck, Proc. 27th ACM Symposium on the Theory of -- Computing, 407-416 (1995). -module Main where--import System.Environment-import Data.Array-import Data.Complex+module Main (main) where import Numeric.Random.Generator.MT19937 import Numeric.Random.Distribution.Uniform -import Numeric.Transform.Fourier.FFT+import Numeric.Transform.Fourier.FFT (fft) +import DSP.Basic ((^!))++import System.Environment (getArgs)++import Data.Complex (Complex((:+)), cis)+import Data.Array (Array, Ix, listArray, elems, bounds, range, (!))++ -- Generates random test vectors gendata :: Int -> W -> Array Int (Complex Double)@@ -26,21 +30,34 @@ -- A few functions over arrays -aadd x y = listArray (0,n) [ x!i + y!i | i <- [0..n] ]- where n = snd $ bounds x+aadd, asub ::+ (Ix i, Num e) =>+ Array i e -> Array i e -> Array i e -asub x y = listArray (0,n) [ x!i - y!i | i <- [0..n] ]- where n = snd $ bounds x+aadd x y = listArray bnds [ x!i + y!i | i <- range bnds ]+ where bnds = bounds x -arot x = listArray (0,n) $ xs' ++ [x']- where xs' = tail $ elems x- x' = head $ elems x- n = snd $ bounds x+asub x y = listArray bnds [ x!i - y!i | i <- range bnds ]+ where bnds = bounds x +arot ::+ (Ix i, Num e) =>+ Array i e -> Array i e++arot xa =+ listArray (bounds xa) $+ case elems xa of+ [] -> []+ x:xs -> xs ++ [x]++ascale ::+ (Ix i, Num e) =>+ e -> Array i e -> Array i e ascale a x = fmap (a*) x -- linearity test: aFFT(x) + bFFT(y) == FFT(ax+by) +lin_test :: Int -> Double lin_test n = acomp z1 z2 where x = gendata n 42 y = gendata n 44@@ -54,6 +71,7 @@ -- impulse response test: rect == FFT(x) + FFT(impulse - x) +imp_test :: Int -> Double imp_test n = acomp a' (aadd b' c') where zeros = 0 : zeros a = listArray (0,n-1) $ (1 :+ 0) : zeros@@ -65,6 +83,7 @@ -- shift test: x[n-m] <-> W_N^km X[k] +shift_test :: Int -> Double shift_test n = acomp a' c' where a = gendata n 42 b = arot a@@ -74,11 +93,15 @@ -- determines peak error (from FFTW) +acomp ::+ (Ix i, RealFloat a) =>+ Array i (Complex a) -> Array i (Complex a) -> a+ acomp x y = (maximum $ zipWith (/) a mag) where a = zipWith calc_a (elems x) (elems y) mag = zipWith calc_mag (elems x) (elems y)- calc_a (xr:+xi) (yr:+yi) = sqrt $ (xr - yr)^2 + (xi - yi)^2- calc_mag (xr:+xi) (yr:+yi) = 0.5 * (sqrt (xr^2+xi^2) + sqrt (yr^2+yi^2)) + tol+ calc_a (xr:+xi) (yr:+yi) = sqrt $ (xr - yr)^!2 + (xi - yi)^!2+ calc_mag (xr:+xi) (yr:+yi) = 0.5 * (sqrt (xr^!2+xi^!2) + sqrt (yr^!2+yi^!2)) + tol tol = 1.0e-6 @@ -90,8 +113,9 @@ where ok = lin_test n < tol && imp_test n < tol && shift_test n < tol tol = 1.0e-6 -testfft :: Int -> Int -> IO [()]-testfft n1 n2 = sequence $ map test1fft [n1..n2]+testfft :: Int -> Int -> IO ()+testfft n1 n2 = mapM_ test1fft [n1..n2] +main :: IO () main = do args <- getArgs testfft (read $ args !! 0) (read $ args !! 1)
demo/FreqDemo.hs view
@@ -14,27 +14,31 @@ -- along with this program; if not, write to the Free Software -- Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -module Main where+module Main (main) where -import Data.Array-import Data.Complex+import DSP.Estimation.Frequency.Pisarenko (pisarenko)+import DSP.Estimation.Frequency.PerMax (permax)+import DSP.Estimation.Frequency.FCI+ (quinn1, quinn2, quinn3, jacobsen, macleod3, macleod5, rv)+import DSP.Estimation.Frequency.QuinnFernandes (qf)+import DSP.Estimation.Frequency.WLP (lrp, kay, lw) -import Numeric+import DSP.Source.Oscillator (nco, quadrature_nco) -import Numeric.Random.Generator.MT19937-import Numeric.Random.Distribution.Normal-import Numeric.Random.Distribution.Uniform+import DSP.Basic ((^!)) -import DSP.Source.Oscillator+import Numeric.Random.Generator.MT19937 (genrand)+import Numeric.Random.Distribution.Uniform (uniform53oc)+import Numeric.Random.Distribution.Normal (normal_ar) -import Numeric.Transform.Fourier.FFT+import Numeric.Transform.Fourier.FFT (rfft) -import DSP.Estimation.Frequency.Pisarenko-import DSP.Estimation.Frequency.PerMax-import DSP.Estimation.Frequency.FCI-import DSP.Estimation.Frequency.QuinnFernandes-import DSP.Estimation.Frequency.WLP+import Numeric (showFFloat) +import Data.Complex (Complex((:+)))+import Data.Array (Array, listArray)++ -- Parameters rho :: Double@@ -54,21 +58,23 @@ -- Vectors -y :: Array Int Double-y = listArray (0,n-1) $ zipWith (+) noise $ map (rho *) $ nco w phi+ya :: Array Int Double+ya = listArray (0,n-1) $ zipWith (+) noise $ map (rho *) $ nco w phi where noise = normal_ar (0, sig2) $ uniform53oc $ genrand 42- sig2 = (rho^2 / 2) / (10 ** (snr / 10))+ sig2 = (rho^!2 / 2) / (10 ** (snr / 10)) -z :: Array Int (Complex Double)-z = listArray (0,n-1) $ zipWith (+) noise $ map ((rho :+ 0) *) $ quadrature_nco w phi+za :: Array Int (Complex Double)+za = listArray (0,n-1) $ zipWith (+) noise $ map ((rho :+ 0) *) $ quadrature_nco w phi where noise = zipWith (:+) (normal_ar (0, sig2) $ uniform53oc $ genrand 42) (normal_ar (0, sig2) $ uniform53oc $ genrand 43)- sig2 = (rho^2 / 2) / (10 ** (snr / 10))+ sig2 = (rho^!2 / 2) / (10 ** (snr / 10)) -- The tests +dfp :: Array Int (Complex Double) -> [(String, Double)] dfp z = [ ("Periodigram Maximizer\t\t\t", permax z k) ] where k = round $ w / 2 / pi * fromIntegral n +fci :: Array Int Double -> [(String, Double)] fci y = [ ("Quinn's First Estimator\t\t\t", quinn1 y' k / 2), ("Quinn's Second Estimator\t\t", quinn2 y' k / 2), ("Quinn's Third Estimator\t\t\t", quinn3 y' k / 2),@@ -79,12 +85,15 @@ where y' = rfft y k = round $ w / 2 / pi * fromIntegral n +scm :: Array Int Double -> [(String, Double)] scm y = [ ("Pisarenko's Method\t\t\t", pisarenko y) ] +offline :: Array Int Double -> [(String, Double)] offline y = [ ("Quinn-Fernandes\t\t\t\t", qf y w') ] where k = round $ w / 2 / pi * fromIntegral n- w' = 2 * pi * fromIntegral k / fromIntegral n+ w' = 2 * pi * fromInteger k / fromIntegral n +fastblock :: Array Int (Complex Double) -> [(String, Double)] fastblock z = [ ("Lank, Reed, and Pollon\t\t\t", lrp z), ("Kay\t\t\t\t\t", kay z), ("Lovell and Williamson\t\t\t", lw z) ]@@ -93,22 +102,25 @@ -- Glue it all together -showone (s,w') = putStrLn $ s ++ ": w=" ++ (showFFloat (Just 6) w' $ " err=" ++ showFFloat (Just 6) (abs (w-w')) "")+showone :: (String, Double) -> IO ()+showone (s,w') =+ putStrLn $ s ++ ": w=" ++ (showFFloat (Just 6) w' $ " err=" ++ showFFloat (Just 6) (abs (w-w')) "") +main :: IO () main = do- putStrLn "==> Parameters"- putStrLn $ "rho=\t" ++ show rho- putStrLn $ "w=\t" ++ show w- putStrLn $ "phi=\t" ++ show phi- putStrLn $ "snr=\t" ++ show snr- putStrLn $ "n=\t" ++ show n- putStrLn "==> Periodigram Techniques"- sequence $ map showone $ dfp z- putStrLn "==> Fourier Coefficient Interpolation Techniques"- sequence $ map showone $ fci y- putStrLn "==> Sample Covariance Methods"- sequence $ map showone $ scm y- putStrLn "==> Offline Filtering Techniques"- sequence $ map showone $ offline y- putStrLn "==> Fast Block Techniques"- sequence $ map showone $ fastblock z+ putStrLn "==> Parameters"+ putStrLn $ "rho=\t" ++ show rho+ putStrLn $ "w=\t" ++ show w+ putStrLn $ "phi=\t" ++ show phi+ putStrLn $ "snr=\t" ++ show snr+ putStrLn $ "n=\t" ++ show n+ putStrLn "==> Periodigram Techniques"+ mapM_ showone $ dfp za+ putStrLn "==> Fourier Coefficient Interpolation Techniques"+ mapM_ showone $ fci ya+ putStrLn "==> Sample Covariance Methods"+ mapM_ showone $ scm ya+ putStrLn "==> Offline Filtering Techniques"+ mapM_ showone $ offline ya+ putStrLn "==> Fast Block Techniques"+ mapM_ showone $ fastblock za
demo/IIRDemo.hs view
@@ -14,29 +14,42 @@ -- along with this program; if not, write to the Free Software -- Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -module Main where+module Main (main) where -import Data.Array+import qualified DSP.Filter.IIR.Design as IIR+import DSP.Filter.IIR.IIR (iir_df1)+import DSP.Source.Basic (impulse) -import DSP.Filter.IIR.IIR-import DSP.Filter.IIR.Design+import Numeric.Transform.Fourier.FFTUtils (write_rfft_info) -import Numeric.Transform.Fourier.FFTUtils+import Data.Array (Array, listArray) -import DSP.Source.Basic -- Examples from Oppenheim and Schafer -ex7'3 = mkButterworth (0.2 * pi, 1 - 0.89125) (0.3 * pi, 0.17783)-ex7'8 = mkChebyshev1 (0.2 * pi, 1 - 0.89125) (0.3 * pi, 0.17783)+ex7'3lp, ex7'3hp, ex7'3bp, ex7'8, ex7'5, ex7'6a, ex7'6b ::+ (Array Int Double, Array Int Double) -ex7'5 = mkButterworth (0.4 * pi, 0.01) (0.6 * pi, 0.001)-ex7'6a = mkChebyshev1 (0.4 * pi, 0.01) (0.6 * pi, 0.001)-ex7'6b = mkChebyshev2 (0.4 * pi, 0.01) (0.6 * pi, 0.001)+ex7'3lp = IIR.butterworthLowpass (0.2 * pi, 1 - 0.89125) (0.3 * pi, 0.17783)+ex7'3hp = IIR.butterworthHighpass (0.2 * pi, 1 - 0.89125) (0.3 * pi, 0.17783)+ex7'3bp = IIR.butterworthBandpass (0.2 * pi, 1 - 0.89125) (0.3 * pi, 0.17783) +ex7'8 = IIR.chebyshev1Lowpass (0.2 * pi, 1 - 0.89125) (0.3 * pi, 0.17783)++ex7'5 = IIR.butterworthLowpass (0.4 * pi, 0.01) (0.6 * pi, 0.001)+ex7'6a = IIR.chebyshev1Lowpass (0.4 * pi, 0.01) (0.6 * pi, 0.001)+ex7'6b = IIR.chebyshev2Lowpass (0.4 * pi, 0.01) (0.6 * pi, 0.001)++example :: String -> (Array Int Double, Array Int Double) -> IO ()+example name coeffs =+ write_rfft_info name $ listArray (0,999) $ iir_df1 coeffs impulse++main :: IO () main = do- write_rfft_info "ex-7.3" $ listArray (0,999) $ iir_df1 ex7'3 $ impulse- write_rfft_info "ex-7.8" $ listArray (0,999) $ iir_df1 ex7'8 $ impulse- write_rfft_info "ex-7.5" $ listArray (0,999) $ iir_df1 ex7'5 $ impulse- write_rfft_info "ex-7.6a" $ listArray (0,999) $ iir_df1 ex7'6a $ impulse- write_rfft_info "ex-7.6b" $ listArray (0,999) $ iir_df1 ex7'6b $ impulse+ example "ex-7.3lp" ex7'3lp+ example "ex-7.3hp" ex7'3hp+ example "ex-7.3bp" ex7'3bp+ example "ex-7.8" ex7'8+ example "ex-7.5" ex7'5+ example "ex-7.6a" ex7'6a+ example "ex-7.6b" ex7'6b
demo/NoiseDemo.hs view
@@ -2,31 +2,32 @@ module Main (main) where --- Import the System functions that we need+-- Import a portion of the Numeric.Random library -import System.Environment-import System.IO-import System.Exit+import Numeric.Random.Generator.MT19937 (genrand)+import Numeric.Random.Distribution.Uniform (uniform53oc)+import Numeric.Random.Distribution.Normal (normal_ar)+import Numeric.Random.Spectrum.White (white)+import Numeric.Random.Spectrum.Pink (kellet)+import Numeric.Random.Spectrum.Purple (purple)+import Numeric.Random.Spectrum.Brown (brown) --- We need support for complex numbers and arrays+-- We do some simple FFT analysis -import Data.Complex-import Data.Array+import Numeric.Transform.Fourier.FFT (rfft) --- Import a portion of the Numeric.Random library+-- Import the System functions that we need -import Numeric.Random.Generator.MT19937-import Numeric.Random.Distribution.Uniform-import Numeric.Random.Distribution.Normal-import Numeric.Random.Spectrum.White-import Numeric.Random.Spectrum.Pink-import Numeric.Random.Spectrum.Purple-import Numeric.Random.Spectrum.Brown+import System.Environment (getProgName, getArgs)+import System.IO (IOMode(WriteMode), withFile, hPutStrLn, hPutStr)+import System.Exit (exitFailure) --- We do some simple FFT analysis+-- We need support for complex numbers and arrays -import Numeric.Transform.Fourier.FFT+import Data.Complex (Complex((:+)))+import Data.Array (Array, listArray, elems, bounds, assocs) + -- Noise parameters mu :: Double@@ -82,11 +83,15 @@ add as bs = listArray (bounds as) $ zipWith (+) (elems as) (elems bs) n = fromIntegral $ length xs --- chunk creates n3 sublists from xs of n1 elemets, and overlapping --- n2 points--chunk :: Int -> Int -> Int -> [Double] -> [[Double]]-chunk n1 n2 n3 xs = take n1 xs : chunk n1 n2 n3 (drop (n1-n2) xs)+{- |+'chunk' creates sublists from xs of n1 elements,+and overlapping n2 points+-}+chunk :: Int -> Int -> [a] -> [[a]]+chunk n1 n2 =+ let m = n1-n2+ go xs = take n1 xs : go (drop m xs)+ in go -- avg calculates an averaged RFFT using a rectangular window -- n1 is the length of each FFT@@ -94,14 +99,14 @@ -- n3 is the number of FFTs to average avgrfft :: Int -> Int -> Int -> [Double] -> Array Int Double-avgrfft n1 n2 n3 xs = avg $ take n3 $ map dbrfft $ map (listArray (0,n1-1)) $ chunk n1 n2 n3 xs+avgrfft n1 n2 n3 xs =+ avg $ take n3 $ map (dbrfft . listArray (0,n1-1)) $ chunk n1 n2 xs -- simple function to write out an array to a file dump :: String -> Array Int Double -> IO ()-dump filename xs = do h <- openFile filename WriteMode- sequence $ map (dump' h) $ assocs $ xs- hClose h+dump filename xs =+ withFile filename WriteMode $ \h -> mapM_ (dump' h) $ assocs xs where dump' h (f,m) = do hPutStr h $ show f hPutStr h $ " " hPutStrLn h $ show m@@ -118,21 +123,18 @@ -- simple function to parse the command line -parseargs :: IO (Int,Int,Int)-parseargs = do args <- getArgs- if length args == 3- then do let n1 = read $ args !! 0- n2 = read $ args !! 1- n3 = read $ args !! 2- return (n1,n2,n3)- else usage+parseArgs :: IO (Int,Int,Int)+parseArgs = do+ args <- getArgs+ case map read args of+ [n1,n2,n3] -> return (n1,n2,n3)+ _ -> usage -- glue it all together main :: IO ()-main = do (n1,n2,n3) <- parseargs+main = do (n1,n2,n3) <- parseArgs dump "white.out" $ avgrfft n1 n2 n3 $ white_gn dump "pink.out" $ avgrfft n1 n2 n3 $ pink_gn dump "brown.out" $ avgrfft n1 n2 n3 $ brown_gn dump "purple.out" $ avgrfft n1 n2 n3 $ purple_gn- return ()
dsp.cabal view
@@ -1,5 +1,5 @@ Name: dsp-Version: 0.2.2+Version: 0.2.3 License: GPL License-File: COPYING Copyright: Matt Donadio, 2003@@ -14,19 +14,11 @@ Tested-With: GHC==6.4.1, GHC==6.8.2 Tested-With: GHC==7.4.2, GHC==7.6.3-Cabal-Version: >=1.6+Cabal-Version: >=1.8 Build-Type: Simple -Extra-Source-Files:- demo/Article.hs- demo/FFTBench.hs- demo/FFTTest.hs- demo/FreqDemo.hs- demo/IIRDemo.hs- demo/NoiseDemo.hs- Source-Repository this- Tag: 0.2.2+ Tag: 0.2.3 Type: darcs Location: http://code.haskell.org/~thielema/dsp/ @@ -35,12 +27,16 @@ Location: http://code.haskell.org/~thielema/dsp/ Flag splitBase- description: Choose the new smaller, split-up base package.+ Description: Choose the new smaller, split-up base package. +Flag buildExamples+ Description: Build demo executables+ Default: True+ Library If flag(splitBase) Build-Depends:- array >=0.1 && <0.5,+ array >=0.1 && <0.6, random >=1.0 && <1.1, base >= 2 && <5 Else@@ -131,10 +127,74 @@ Other-Modules: Numeric.Transform.Fourier.Eigensystem --- Executable:--- Article.hs--- FFTBench.hs--- FFTTest.hs--- FreqDemo.hs--- IIRDemo.hs--- NoiseDemo.hs+Executable dsp-demo-article+ Main-Is: Article.hs+ Hs-Source-Dirs: demo+ GHC-Options: -Wall+ If flag(buildExamples)+ Build-Depends:+ dsp,+ array,+ base+ Else+ Buildable: False++Executable dsp-demo-fft-bench+ Main-Is: FFTBench.hs+ Hs-Source-Dirs: demo+ GHC-Options: -Wall+ If flag(buildExamples)+ Build-Depends:+ dsp,+ array,+ base+ Else+ Buildable: False++Executable dsp-demo-fft-test+ Main-Is: FFTTest.hs+ Hs-Source-Dirs: demo+ GHC-Options: -Wall+ If flag(buildExamples)+ Build-Depends:+ dsp,+ array,+ base+ Else+ Buildable: False++Executable dsp-demo-freq+ Main-Is: FreqDemo.hs+ Hs-Source-Dirs: demo+ GHC-Options: -Wall+ If flag(buildExamples)+ Build-Depends:+ dsp,+ array,+ base+ Else+ Buildable: False++Executable dsp-demo-iir+ Main-Is: IIRDemo.hs+ Hs-Source-Dirs: demo+ GHC-Options: -Wall+ If flag(buildExamples)+ Build-Depends:+ dsp,+ array,+ base+ Else+ Buildable: False++Executable dsp-demo-noise+ Main-Is: NoiseDemo.hs+ Hs-Source-Dirs: demo+ GHC-Options: -Wall+ If flag(buildExamples)+ Build-Depends:+ dsp,+ array,+ base+ Else+ Buildable: False