-- Hoogle documentation, generated by Haddock -- See Hoogle, http://www.haskell.org/hoogle/ -- | Haskell Digital Signal Processing -- -- Digital Signal Processing, Fourier Transform, Linear Algebra, -- Interpolation @package dsp @version 0.2.5.1 -- | Module to perform the linear convolution of two sequences module DSP.Convolution -- | conv convolves two finite sequences conv :: (Ix a, Integral a, Num b) => Array a b -> Array a b -> Array a b test :: Bool -- | 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 -- | raw cross-correllation rxy :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Complex b -- | biased cross-correllation rxy_b :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Complex b -- | unbiased cross-correllation rxy_u :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Complex b -- | raw auto-correllation rxx :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Complex b -- | biased auto-correllation rxx_b :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Complex b -- | unbiased auto-correllation rxx_u :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Complex b test :: Bool -- | 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 best -- here, but it also requires three (x,f(x)) pairs module DSP.Estimation.Frequency.PerMax -- | Discrete frequency periodigram maximizer permax :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | This is an implementation of the Quinn-Fernandes algorithm for -- estimating the frequency of a real sinusoid in noise. module DSP.Estimation.Frequency.QuinnFernandes -- | The Quinn-Fernandes algorithm qf :: (Ix a, Integral a, RealFloat b) => Array a b -> b -> b -- | This module contains a few algorithms for ARMA parameter estimation. -- Algorithms are taken from Steven M. Kay, _Modern Spectral Estimation: -- Theory and Application_, which is one of the standard texts on the -- subject. When possible, variable conventions are the same in the code -- as they are found in the text. -- -- BROKEN: DO NOT USE module DSP.Estimation.Spectral.ARMA -- | THIS DOES NOT WORK arma_mywe :: (RealFloat b, Integral i, Ix i) => Array i (Complex b) -> i -> i -> Array i (Complex b) -- | Test vectors from Kay, Modern Spectral Estimation module DSP.Estimation.Spectral.KayData -- | Complex test data xc :: Array Int (Complex Double) -- | Real test data xr :: Array Int Double -- | Finite Impuse Response filtering functions module DSP.Filter.FIR.FIR -- | 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 -- the type of the filter. In the following functions, we use the O&S -- convention that m is the order of the filter, which is equal to the -- number of taps minus one. fir :: (Num a, Eq a) => Array Int a -> [a] -> [a] test :: Bool -- | Functions for creating rectangular windowed FIR filters module DSP.Filter.FIR.Taps -- | Lowpass filter lpf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -> a -> Array a b -- | Highpass filter hpf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -> a -> Array a b -- | Bandpass filter bpf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -> b -> a -> Array a b -- | Bandstop filter bsf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -> b -> a -> Array a b -- | Multiband filter mbf :: (Ix a, Integral a, Enum b, Floating b, Eq b) => [b] -> [b] -> a -> Array a b -- | Raised-cosine filter rc :: (Ix a, Integral a, Enum b, Floating b, Eq b) => b -> b -> a -> Array a b -- | IIR functions -- -- IMPORTANT NOTE: -- -- 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) = ------------------------
--   
-- --
--   1 + sum(k=1..N) a_k*z^-1
--   
module DSP.Filter.IIR.IIR -- | This is an integrator when a==1, and a leaky integrator when 0 -- < a < 1. -- --
--   y[n] = a * y[n-1] + x[n]
--   
integrator :: Num a => a -> [a] -> [a] -- | First order section, DF1 -- --
--   v[n] = b0 * x[n] + b1 * x[n-1]
--   
-- --
--   y[n] = v[n] - a1 * y[n-1]
--   
fos_df1 :: Num a => a -> a -> a -> [a] -> [a] -- | First order section, DF2 -- --
--   w[n] = -a1 * w[n-1] + x[n]
--   
-- --
--   y[n] = b0 * w[n] + b1 * w[n-1]
--   
fos_df2 :: Num a => a -> a -> a -> [a] -> [a] -- | First order section, DF2T -- --
--   v0[n] = b0 * x[n] + v1[n-1]
--   
-- --
--   y[n] = v0[n]
--   
-- --
--   v1[n] = -a1 * y[n] + b1 * x[n]
--   
fos_df2t :: Num a => a -> a -> a -> [a] -> [a] -- | Direct Form I for a second order section -- --
--   v[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2]
--   
-- --
--   y[n] = v[n] - a1 * y[n-1] - a2 * y[n-2]
--   
biquad_df1 :: Num a => a -> a -> a -> a -> a -> [a] -> [a] -- | Direct Form II for a second order section (biquad) -- --
--   w[n] = -a1 * w[n-1] - a2 * w[n-2] + x[n]
--   
-- --
--   y[n] = b0 * w[n] + b1 * w[n-1] + b2 * w[n-2]
--   
biquad_df2 :: Num a => a -> a -> a -> a -> a -> [a] -> [a] -- | Transposed Direct Form II for a second order section -- --
--   v0[n] = b0 * x[n] + v1[n-1]
--   
-- --
--   y[n] = v0[n]
--   
-- --
--   v1[n] = -a1 * y[n] + b1 * x[n] + v2[n-1]
--   
-- --
--   v2[n] = -a2 * y[n] + b2 * x[n]
--   
biquad_df2t :: Num a => a -> a -> a -> a -> a -> [a] -> [a] -- | Direct Form I IIR -- --
--   v[n] = sum(k=0..M) b_k*x[n-k]
--   
-- --
--   y[n] = v[n] - sum(k=1..N) a_k*y[n-k]
--   
-- -- v[n] is calculated with fir iir_df1 :: (Num a, Eq a) => (Array Int a, Array Int a) -> [a] -> [a] -- | Direct Form II IIR -- --
--   w[n] = x[n] - sum(k=1..N) a_k*w[n-k]
--   
-- --
--   y[n] = sum(k=0..M) b_k*w[n-k]
--   
iir_df2 :: Num a => (Array Int a, Array Int a) -> [a] -> [a] xt :: [Double] yt :: [Double] f1 :: Fractional a => [a] -> [a] f2 :: Fractional a => [a] -> [a] f3 :: Fractional a => [a] -> [a] f4 :: [Double] -> [Double] f5 :: [Double] -> [Double] -- | Flowgraph functions -- -- DO NOT USE YET module DSP.Flowgraph -- | Cascade of functions, eg -- --
--   cascade [ f1, f2, f3 ] x == (f3 . f2 . f1) x
--   
cascade :: Num a => [[a] -> [a]] -> [a] -> [a] -- | Gain node -- --
--   y[n] = a * x[n]
--   
gain :: Num a => a -> [a] -> [a] -- | Bias node -- --
--   y[n] = x[n] + a
--   
bias :: Num a => a -> [a] -> [a] -- | Adder node -- --
--   z[n] = x[n] + y[n]
--   
adder :: Num a => [a] -> [a] -> [a] -- | Polyphase interpolators and decimators -- -- Reference: C&R module DSP.Multirate.Polyphase -- | Polyphase interpolator poly_interp :: (Num a, Eq a) => Int -> Array Int a -> [a] -> [a] -- | Basic signals module DSP.Source.Basic -- | all zeros zeros :: Num a => [a] -- | single impulse impulse :: Num a => [a] -- | unit step step :: Num a => [a] -- | ramp ramp :: Num a => [a] -- | Basic functions for manipulating signals module DSP.Basic -- | linspace generates a list of values linearly spaced between -- specified start and end values (array will include both start and end -- values). -- --
--   linspace 0.0 1.0 5 == [ 0.0, 0.25, 0.5, 0.75 1.0 ]
--   
linspace :: Double -> Double -> Int -> [Double] -- | logspace generates a list of values logarithmically spaced -- between the values 10 ** start and 10 ** end (array will include both -- start and end values). -- --
--   logspace 0.0 1.0 4 == [ 1.0, 2.1544, 4.6416, 10.0 ]
--   
logspace :: Double -> Double -> Int -> [Double] -- | delay is the unit delay function, eg, -- --
--   delay1 [ 1, 2, 3 ] == [ 0, 1, 2, 3 ]
--   
delay1 :: Num a => [a] -> [a] -- | delay is the n sample delay function, eg, -- --
--   delay 3 [ 1, 2, 3 ] == [ 0, 0, 0, 1, 2, 3 ]
--   
delay :: Num a => Int -> [a] -> [a] -- | downsample throws away every n'th sample, eg, -- --
--   downsample 2 [ 1, 2, 3, 4, 5, 6 ] == [ 1, 3, 5 ]
--   
downsample :: Int -> [a] -> [a] downsampleRec :: Int -> [a] -> [a] -- | 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] upsampleRec :: Num a => Int -> [a] -> [a] -- | upsampleAndHold replicates each sample n times, eg, -- --
--   upsampleAndHold 3 [ 1, 2, 3 ] == [ 1, 1, 1, 2, 2, 2, 3, 3, 3 ]
--   
upsampleAndHold :: Int -> [a] -> [a] -- | 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] -- | 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 split. uninterleave :: [a] -> ([a], [a]) -- | 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 -- | generates a Just if the given condition holds toMaybe :: Bool -> a -> Maybe a -- | Computes the square of the Euclidean norm of a 2D point norm2sqr :: Num a => (a, a) -> a -- | Power with fixed exponent type. This eliminates warnings about using -- default types. (^!) :: Num a => a -> Int -> a infixr 8 ^! -- | Halfband interpolators and decimators -- -- Reference: C&R module DSP.Multirate.Halfband -- | Halfband interpolator hb_interp :: (Num a, Eq a) => Array Int a -> [a] -> [a] -- | Halfband decimator hb_decim :: (Num a, Eq a) => Array Int a -> [a] -> [a] -- | CIC filters -- -- R = rate change -- -- M = differential delay in combs -- -- N = number of stages module DSP.Multirate.CIC -- | CIC interpolator cic_interpolate :: Num a => Int -> Int -> Int -> [a] -> [a] -- | CIC interpolator cic_decimate :: Num a => Int -> Int -> Int -> [a] -> [a] -- | Cookbook formulae for audio EQ biquad filter coefficients by Robert -- Bristow-Johnson robert@wavemechanics.com -- -- -- http://www.harmony-central.com/Computer/Programming/Audio-EQ-Cookbook.txt module DSP.Filter.IIR.Cookbook lpf :: Floating a => a -> a -> [a] -> [a] hpf :: Floating a => a -> a -> [a] -> [a] bpf_csg :: Floating a => a -> a -> [a] -> [a] bpf_cpg :: Floating a => a -> a -> [a] -> [a] notch :: Floating a => a -> a -> [a] -> [a] apf :: Floating a => a -> a -> [a] -> [a] peakingEQ :: Floating a => a -> a -> a -> [a] -> [a] lowShelf :: Floating a => a -> a -> a -> [a] -> [a] highShelf :: Floating a => a -> a -> a -> [a] -> [a] -- | 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)) -- -- Procedure: -- --
    --
  1. Filter the signal once with H(z)
  2. --
  3. Double this
  4. --
  5. Subtract this from 3x
  6. --
  7. Filter this twice by H(z) or once by H(z)^2
  8. --
module DSP.Filter.FIR.Sharpen -- | Filter shaprening routine sharpen :: (Num a, Eq a) => Array Int a -> [a] -> [a] -- | This module contains a few algorithms for weighted linear predictors -- for estimating the frequency of a complex sinusoid in noise. module DSP.Estimation.Frequency.WLP -- | The weighted linear predictor form of the frequency estimator wlp :: (Ix a, Integral a, RealFloat b) => Array a b -> Array a (Complex b) -> b -- | WLP using Lank, Reed, and Pollon's window lrp :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> b -- | WLP using kay's window kay :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> b -- | WLP using Lovell and Williamson's window lw :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> b -- | WLP using Clarkson, Kootsookos, and Quinn's window ckq :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> b -> b -> b -- | This module contains an implementation of Pisarenko Harmonic -- Decomposition for a single real sinusoid. For this case, eigenvalues -- do not need to be computed. module DSP.Estimation.Frequency.Pisarenko -- | Pisarenko's method for a single sinusoid pisarenko :: (Ix a, Integral a, Floating b) => Array a b -> b -- | This module contains a few simple algorithms for interpolating the -- peak location of a DFT/FFT. module DSP.Estimation.Frequency.FCI -- | Quinn's First Estimator (FCI1) quinn1 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | Quinn's Second Estimator (FCI2) quinn2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | Quinn's Third Estimator (FCI3) quinn3 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | Eric Jacobsen's Estimator jacobsen :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | MacLeod's Three Point Estimator macleod3 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | MacLeod's Three Point Estimator macleod5 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | Rife and Vincent's Estimator rv :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> b -- | NCO and NCOM functions module DSP.Source.Oscillator -- | 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, 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. nco :: RealFloat a => a -> a -> [a] -- | ncom mixes (multiplies) x by a real sine wave with normalized -- frequency wn. This is usually called an NCOM: Numerically Controlled -- Oscillator and Modulator. ncom :: RealFloat a => a -> a -> [a] -> [a] -- | quadrature_nco returns an infinite list representing a complex -- phasor with a phase step of wn radians, ie a quadrature nco with -- normalized frequency wn radians/sample. Since Haskell uses lazy -- evaluation, rotate will only be computed once, so this NCO uses only -- one sin and one cos for the entire list, at the expense of 4 mults, 1 -- add, and 1 subtract per point. quadrature_nco :: RealFloat a => a -> a -> [Complex a] -- | complex_ncom mixes the complex input x with a quardatue nco -- with normalized frequency wn radians/sample using complex multiplies -- (perform a complex spectral shift) complex_ncom :: RealFloat a => a -> a -> [Complex a] -> [Complex a] -- | quadrature_ncom mixes the complex input x with a quadrature nco -- with normalized frequency wn radians/sample in quadrature (I/Q -- modulation) quadrature_ncom :: RealFloat a => a -> a -> [Complex a] -> [a] agc :: RealFloat a => Complex a -> Complex a -- | Simple phase unwrapping algorithm module DSP.Unwrap -- | This is the simple phase unwrapping algorithm from Oppenheim and -- Schafer. unwrap :: (Ix a, Integral a, Ord b, Floating b) => b -> Array a b -> Array a b -- | Commonly used window functions. Except for the Parzen window, the -- results of all of these look right, but I have to check them -- against either Matlab or my C code. -- -- More windowing functions exist, but I have to dig through my papers to -- find the equations. module DSP.Window -- | Applys a window, w, to a sequence x window :: Array Int Double -> Array Int Double -> Array Int Double -- | rectangular window rectangular :: Int -> Array Int Double -- | Bartlett window bartlett :: Int -> Array Int Double -- | Hanning window hanning :: Int -> Array Int Double -- | Hamming window hamming :: Int -> Array Int Double -- | Blackman window blackman :: Int -> Array Int Double -- | rectangular window kaiser :: Double -> Int -> Array Int Double -- | Generalized Hamming window gen_hamming :: Double -> Int -> Array Int Double -- | rectangular window parzen :: Int -> Array Int Double -- | Deprecated: Use DSP.Window instead module DSP.Filter.FIR.Window -- | This module implements the Kaiser Window Method for designing FIR -- filters. module DSP.Filter.FIR.Kaiser -- | Designs a lowpass Kaiser filter kaiser_lpf :: Double -> Double -> Double -> Double -> Array Int Double -- | Designs a highpass Kaiser filter kaiser_hpf :: Double -> Double -> Double -> Double -> Array Int Double -- | This module contains a routine that solves the system Ax=b, where A is -- positive definite, using Cholesky decomposition. module Matrix.Cholesky cholesky :: (Ix a, Integral a, RealFloat b) => Array (a, a) (Complex b) -> Array a (Complex b) -> Array a (Complex b) -- | This module contains an implementation of Levinson-Durbin recursion. module Matrix.Levinson -- | levinson takes an array, r, of autocorrelation values, and a model -- order, p, and returns an array, a, of the model estimate and rho, the -- noise power. levinson :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b), b) -- | This module contains a few algorithms for AR parameter estimation. -- Algorithms are taken from Steven M. Kay, /Modern Spectral Estimation: -- Theory and Application/, which is one of the standard texts on the -- subject. When possible, variable conventions are the same in the code -- as they are found in the text. module DSP.Estimation.Spectral.AR -- | Computes an AR(p) model estimate from x using the Yule-Walker method ar_yw :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b), b) -- | Computes an AR(p) model estimate from x using the covariance method ar_cov :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b), b) -- | Computes an AR(p) model estimate from x using the modified covariance -- method ar_mcov :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b), b) -- | Computes an AR(p) model estimate from x using the Burg' method ar_burg :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b), b) -- | This module contains one algorithm for MA parameter estimation. It is -- taken from Steven M. Kay, _Modern Spectral Estimation: Theory and -- Application_, which is one of the standard texts on the subject. When -- possible, variable conventions are the same in the code as they are -- found in the text. module DSP.Estimation.Spectral.MA -- | Computes an MA(q) model estimate from x using the Durbin's method -- where l is the order of the AR process used in the algorithm ma_durbin :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> a -> (Array a (Complex b), b) -- | Two-step simplex algorithm -- -- I only guarantee that this module wastes inodes module Matrix.Simplex -- | Type for results of the simplex algorithm data Simplex a Unbounded :: Simplex a Infeasible :: Simplex a Optimal :: a -> Simplex a -- | The simplex algorithm for standard form: -- -- min c'x -- -- where Ax = b, x >= 0 -- -- a!(0,0) = -z -- -- a!(0,j) = c' -- -- a!(i,0) = b -- -- a!(i,j) = A_ij simplex :: Array (Int, Int) Double -> Simplex (Array (Int, Int) Double) -- | The two-phase simplex algorithm twophase :: Array (Int, Int) Double -> Simplex (Array (Int, Int) Double) instance GHC.Show.Show a => GHC.Show.Show (Matrix.Simplex.Simplex a) instance GHC.Read.Read a => GHC.Read.Read (Matrix.Simplex.Simplex a) module Matrix.Vector generate :: Ix i => (i, i) -> (i -> a) -> Array i a fromList :: [a] -> Array Int a toList :: Array Int a -> [a] norm :: (Ix i, Floating a) => Array i a -> a scale :: (Ix i, Num a) => a -> Array i a -> Array i a lift2 :: Ix i => (a -> b -> c) -> Array i a -> Array i b -> Array i c add :: (Ix i, Num a) => Array i a -> Array i a -> Array i a sub :: (Ix i, Num a) => Array i a -> Array i a -> Array i a module Matrix.Sparse data Matrix i j a bounds :: Matrix i j a -> ((i, j), (i, j)) fromMap :: (Ord i, Ord j) => ((i, j), (i, j)) -> Map (i, j) a -> Matrix i j a fromRows :: (Ord i, Ord j) => ((i, j), (i, j)) -> Map i (Map j a) -> Matrix i j a fromColumns :: (Ord i, Ord j) => ((i, j), (i, j)) -> Map j (Map i a) -> Matrix i j a fromDense :: (Ix i, Ix j) => Array (i, j) a -> Matrix i j a toRows :: (Ord i, Ord j) => Matrix i j a -> Map i (Map j a) toColumns :: (Ord i, Ord j) => Matrix i j a -> Map j (Map i a) toDense :: (Ix i, Ix j, Num a) => Matrix i j a -> Array (i, j) a getRow :: (Ord i, Ord j) => i -> Matrix i j a -> Map j a getColumn :: (Ord i, Ord j) => j -> Matrix i j a -> Map i a mulVector :: (Ix i, Ix j, Num a) => Matrix i j a -> Array j a -> Array i a instance (GHC.Show.Show i, GHC.Show.Show j, GHC.Show.Show a) => GHC.Show.Show (Matrix.Sparse.Matrix i j a) instance GHC.Base.Functor (Matrix.Sparse.Matrix i j) module Matrix.QR.Givens -- | Solve a sparse overconstrained linear problem, i.e. minimize -- ||Ax-b||. A must have dimensions m x n with -- m>=n and it must have full-rank. None of these conditions -- is checked. leastSquares :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Matrix i j a -> Array i a -> Array j a -- | The decomposition routine is pretty simple. It does not try to -- minimize fill-up by a clever ordering of rotations. However, for -- banded matrices it will work as expected. decompose :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Matrix i j a -> ([Rotation i a], Upper i j a) solve :: (Ix i, Ix j, Fractional a) => ([Rotation i a], Upper i j a) -> Array i a -> Array j a -- | Only sensible for square matrices, but the function does not check -- whether the matrix is square. det :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Matrix i j a -> a -- | Absolute value of the determinant. This is also sound for non-square -- matrices. detAbs :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Matrix i j a -> a data Rotation i a rotateVector :: (Ix i, Num a) => Rotation i a -> Array i a -> Array i a data Upper i j a -- | Assumes that Upper matrix is at least as high as wide and that -- it has full rank. solveUpper :: (Ix i, Ix j, Fractional a) => Upper i j a -> Array i a -> Array j a detUpper :: (Ix i, Ix j, Fractional a) => Upper i j a -> a instance (GHC.Show.Show i, GHC.Show.Show j, GHC.Show.Show a) => GHC.Show.Show (Matrix.QR.Givens.Upper i j a) instance (GHC.Show.Show i, GHC.Show.Show a) => GHC.Show.Show (Matrix.QR.Givens.Rotation i a) -- | Basic matrix routines module Matrix.Matrix -- | Matrix-matrix multiplication: A x B = C mm_mult :: (Ix i, Ix j, Ix k, Num a) => Array (i, j) a -> Array (j, k) a -> Array (i, k) a -- | Matrix-vector multiplication: A x b = c mv_mult :: (Ix i, Ix j, Num a) => Array (i, j) a -> Array j a -> Array i a -- | Transpose of a matrix m_trans :: (Ix i, Ix j, Num a) => Array (i, j) a -> Array (j, i) a -- | Hermitian transpose (conjugate transpose) of a matrix m_hermit :: (Ix i, Ix j, RealFloat a) => Array (i, j) (Complex a) -> Array (j, i) (Complex a) columnBounds :: (Ix i, Ix j) => Array (i, j) a -> (i, i) rowBounds :: (Ix i, Ix j) => Array (i, j) a -> (j, j) getColumn :: (Ix i, Ix j) => j -> Array (i, j) e -> Array i e getRow :: (Ix i, Ix j) => i -> Array (i, j) e -> Array j e toColumns :: (Ix i, Ix j) => Array (i, j) a -> [Array i a] toRows :: (Ix i, Ix j) => Array (i, j) a -> [Array j a] -- | We need the bounds of the row indices for empty input lists. fromColumns :: Ix i => (i, i) -> [Array i a] -> Array (i, Int) a fromRows :: Ix j => (j, j) -> [Array j a] -> Array (Int, j) a outer :: (Ix i, Ix j, Num a) => Array i a -> Array j a -> Array (i, j) a inner :: (Ix i, Num a) => Array i a -> Array i a -> a module Matrix.QR.Householder -- | Solve an overconstrained linear problem, i.e. minimize -- ||Ax-b||. A must have dimensions m x n with -- m>=n and it must have full-rank. None of these conditions -- is checked. leastSquares :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Array (i, j) a -> Array i a -> Array j a decompose :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Array (i, j) a -> ([Reflection i a], Upper i j a) solve :: (Ix i, Ix j, Fractional a) => ([Reflection i a], Upper i j a) -> Array i a -> Array j a -- | Only sensible for square matrices, but the function does not check -- whether the matrix is square. -- -- For non-square matrices the sign is not of much use. It depends on the -- implementation. It is not uniquely defined by requiring that the -- determinant of the orthogonal transformation has positive sign. det :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Array (i, j) a -> a -- | Absolute value of the determinant. This is also sound for non-square -- matrices. detAbs :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Array (i, j) a -> a data Reflection i a reflectMatrix :: (Ix i, Ix j, Num a) => Reflection i a -> Array (i, j) a -> Array (i, j) a reflectVector :: (Ix i, Num a) => Reflection i a -> Array i a -> Array i a data Upper i j a matrixFromUpper :: (Ix i, Ix j, Num a) => Upper i j a -> Array (i, j) a -- | Assumes that Upper matrix is at least as high as wide and that -- it has full rank. solveUpper :: (Ix i, Ix j, Fractional a) => Upper i j a -> Array i a -> Array j a detUpper :: (Ix i, Ix j, Fractional a) => Upper i j a -> a -- | Module implementing LU decomposition and related functions module Matrix.LU -- | LU decomposition via Crout's Algorithm lu :: Array (Int, Int) Double -> Array (Int, Int) Double -- | Solution to Ax=b via LU decomposition lu_solve :: Array (Int, Int) Double -> Array Int Double -> Array Int Double -- | Improve a solution to Ax=b via LU decomposition improve :: Array (Int, Int) Double -> Array (Int, Int) Double -> Array Int Double -> Array Int Double -> Array Int Double -- | Matrix inversion via LU decomposition inverse :: Array (Int, Int) Double -> Array (Int, Int) Double -- | Determinant of a matrix via LU decomposition lu_det :: Array (Int, Int) Double -> Double -- | LU solver using original matrix solve :: Array (Int, Int) Double -> Array Int Double -> Array Int Double -- | Determinant computation by implicit LU decomposition with -- permutations. det :: Array (Int, Int) Double -> Double -- | 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 module DSP.Filter.IIR.Prony -- | Implementation of Prony's method prony :: Int -> Int -> Array Int Double -> (Array Int Double, Array Int Double) -- | Function approximation using Chebyshev polynomials -- --
--   f(x) = ( sum (k=0..N-1) c_k * T_k(x) ) - 0.5 * c_0
--   
-- -- over the interval [a,b] -- -- Reference: NRiC module Numeric.Approximation.Chebyshev -- | Calculates the Chebyshev approximation to f(x) over -- [a,b] cheby_approx :: (Double -> Double) -> Double -> Double -> Int -> [Double] -- | Evaluates the Chebyshev approximation to f(x) over -- [a,b] at x cheby_eval :: [Double] -> Double -> Double -> Double -> Double -- | UNTESTED -- -- Module for transforming a list of uniform random variables into a list -- of binomial random variables. -- -- Reference: Ross module Numeric.Random.Distribution.Binomial -- | Generates a list of binomial random variables from a list of uniforms binomial :: Int -> Double -> [Double] -> [Double] -- | UNTESTED -- -- Module for transforming a list of uniform random variables into a list -- of exponential random variables. -- --
--   f(x) = lambda * exp(-lambda*x)
--   
-- --
--   F(x) = 1 - exp(-lambda*x)
--   
-- --
--   lambda = 1 / mu
--   
-- -- Reference: Ross module Numeric.Random.Distribution.Exponential -- | Generates a list of exponential random variables from a list of -- uniforms via the inverse transformation method -- --
--   F(x) = 1 - exp(-lambda*x)
--   
-- --
--   F^-1(x) = -log(1 - x) / lambda
--   
exponential_inv :: Double -> [Double] -> [Double] -- | UNTESTED -- -- Module for transforming a list of uniform random variables into a list -- of gamma random variables. -- --
--   f(x) = lambda * exp(-lambda*x) * (lambda * x)^(t-1) / Gamma(t)
--   
-- -- Reference: Ross module Numeric.Random.Distribution.Gamma -- | Generates a list of gamma random variables from a list of uniforms via -- the inverse transformation method gamma :: Int -> Double -> [Double] -> [Double] -- | UNTESTED -- -- Module for transforming a list of uniform random variables into a list -- of geometric random variables. -- --
--   P{X=n} = (1-p)^(n-1)*p
--   
-- -- Reference: Ross module Numeric.Random.Distribution.Geometric -- | Generates a list of geometric random variables from a list of uniforms geometric :: Double -> [Double] -> [Double] -- | Module for transforming a list of uniform random variables into a list -- of normal random variables. module Numeric.Random.Distribution.Normal -- | Normal random variables via the Central Limit Theorm (not explicity -- given, but see Ross) -- -- If mu=0 and sigma=1, then this will generate numbers in the range -- [-n2,n2] normal_clt :: Int -> (Double, Double) -> [Double] -> [Double] -- | 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) -> [Double] -> [Double] -- | 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) -> [Double] -> [Double] -- | 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) -> [Double] -> [Double] -- | Functions for turning a list of random integers (as Word32) in -- a list of Uniform RV's module Numeric.Random.Distribution.Uniform -- | 32 bits in [0,1] uniform32cc :: [Word32] -> [Double] -- | 32 bits in [0,1) uniform32co :: [Word32] -> [Double] -- | 32 bits in (0,1] uniform32oc :: [Word32] -> [Double] -- | 32 bits in (0,1) uniform32oo :: [Word32] -> [Double] -- | 53 bits in [0,1], ie 64-bit IEEE 754 in [0,1] uniform53cc :: [Word32] -> [Double] -- | 53 bits in [0,1), ie 64-bit IEEE 754 in [0,1) uniform53co :: [Word32] -> [Double] -- | 53 bits in (0,1] uniform53oc :: [Word32] -> [Double] -- | 53 bits in (0,1) uniform53oo :: [Word32] -> [Double] -- | transforms uniform [0,1] to [a,b] uniform :: Double -> Double -> [Double] -> [Double] -- | A Haskell program for MT19937 pseudorandom number generator module Numeric.Random.Generator.MT19937 type W = Word32 genrand :: W -> [W] test :: IO () -- | Function for brown noise, which is integrated white noise module Numeric.Random.Spectrum.Brown brown :: [Double] -> [Double] -- | Functions for pinking noise -- -- http://www.firstpr.com.au/dsp/pink-noise/ module Numeric.Random.Spectrum.Pink -- | Kellet's filter kellet :: [Double] -> [Double] -- | Voss's algorithm -- -- UNTESTED, but the algorithm looks like it is working based on my hand -- tests. voss :: Int -> [Double] -> [Double] -- | Function for purple noise, which is differentiated white noise -- -- This currently just does a simple first-order difference. This is -- equivalent to filtering the white noise with h[n] = [1,-1] -- A better solution would be to use a proper FIR differentiator. module Numeric.Random.Spectrum.Purple purple :: [Double] -> [Double] -- | Function for white noise -- -- This is pretty useless, but it is here to be comprehensive module Numeric.Random.Spectrum.White white :: [Double] -> [Double] module Numeric.Special.Trigonometric csc :: Floating a => a -> a sec :: Floating a => a -> a cot :: Floating a => a -> a acsc :: Floating a => a -> a asec :: Floating a => a -> a acot :: Floating a => a -> a csch :: Floating a => a -> a sech :: Floating a => a -> a coth :: Floating a => a -> a acsch :: Floating a => a -> a asech :: Floating a => a -> a acoth :: Floating a => a -> a -- | Simple module for computing the median on a list -- -- Reference: Ross, NRiC module Numeric.Statistics.Median -- | Compute the median of a list median :: (Ord a, Fractional a) => [a] -> a -- | Compute the center of the list in a more lazy manner and thus halves -- memory requirement. medianFast :: (Ord a, Fractional a) => [a] -> a -- | Simple module for computing the various moments of a list -- -- Reference: Ross, NRiC module Numeric.Statistics.Moment -- | Compute the mean of a list -- --
--   Mean(X) = 1/N sum(i=1..N) x_i
--   
mean :: Fractional a => [a] -> a -- | Compute the variance of a list -- --
--   Var(X) = sigma^2
--   
-- --
--   = 1/N-1 sum(i=1..N) (x_i-mu)^2
--   
var :: Fractional a => [a] -> a -- | Compute the standard deviation of a list -- --
--   StdDev(X) = sigma = sqrt (Var(X))
--   
stddev :: RealFloat a => [a] -> a -- | Compute the average deviation of a list -- --
--   AvgDev(X) = 1/N sum(i=1..N) |x_i-mu|
--   
avgdev :: RealFloat a => [a] -> a -- | Compute the skew of a list -- --
--   Skew(X) = 1/N sum(i=1..N) ((x_i-mu)/sigma)^3
--   
skew :: RealFloat a => [a] -> a -- | Compute the kurtosis of a list -- --
--   Kurt(X) = ( 1/N sum(i=1..N) ((x_i-mu)/sigma)^4 ) - 3
--   
kurtosis :: RealFloat a => [a] -> a -- | UNTESTED -- -- Simple module for computing the covariance of two lists -- --
--   Cov(X1,X2) = 1/(N-1) * sum (i=1..N) ((x1_i - mu1)(x2_i - mu2))
--   
-- -- Reference: Ross, NRiC module Numeric.Statistics.Covariance cov :: Fractional a => [a] -> [a] -> a -- | UNTESTED -- -- 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 -- | Generates a list of poisson random variables from a list of uniforms. poisson :: Double -> [Double] -> [Int] test :: Int -> Double -> Double testHead :: Int -> Double -> Double -- | 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. module DSP.Covariance -- | raw cross-covariance -- -- We define covariance in terms of correlation. -- -- Cxy(X,Y) = E[(X - E[X])(Y - E[Y])] = E[XY] - E[X]E[Y] = Rxy(X,Y) - -- E[X]E[Y] cxy :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Complex b -- | biased cross-covariance cxy_b :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Complex b -- | unbiased cross-covariance cxy_u :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -> a -> Complex b -- | raw auto-covariance -- -- Cxx(X,X) = E[(X - E[X])(X - E[X])] = E[XX] - E[X]E[X] = Rxy(X,X) - -- E[X]^2 cxx :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Complex b -- | biased auto-covariance cxx_b :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Complex b -- | unbiased auto-covariance cxx_u :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Complex b -- | UNTESTED: DO NOT USE -- -- Student's t-test functions -- -- Reference: NRiC module Numeric.Statistics.TTest ttest :: [Double] -> [Double] -> Double tutest :: [Double] -> [Double] -> Double tptest :: [Double] -> [Double] -> Double -- | Cooley-Tukey algorithm for computing the FFT module Numeric.Transform.Fourier.CT -- | Cooley-Tukey algorithm doing row FFT's then column FFT's fft_ct1 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Cooley-Tukey algorithm doing column FFT's then row FFT's fft_ct2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Not so naive implementation of a Discrete Fourier Transform. module Numeric.Transform.Fourier.DFT dft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -- | Hard-coded FFT transforms module Numeric.Transform.Fourier.FFTHard -- | Length 2 FFT fft'2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -- | Length 3 FFT fft'3 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -- | Length 4 FFT fft'4 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -- | 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. module Numeric.Transform.Fourier.Goertzel -- | Goertzel's algorithm for complex inputs cgoertzel :: (RealFloat a, Ix b, Integral b) => Array b (Complex a) -> b -> Complex a -- | Power via Goertzel's algorithm for complex inputs cgoertzel_power :: (RealFloat a, Ix b, Integral b) => Array b (Complex a) -> b -> a -- | Goertzel's algorithm for real inputs rgoertzel :: (RealFloat a, Ix b, Integral b) => Array b a -> b -> Complex a -- | Power via Goertzel's algorithm for real inputs rgoertzel_power :: (RealFloat a, Ix b, Integral b) => Array b a -> b -> a -- | Prime Factor Algorithm module Numeric.Transform.Fourier.PFA -- | Prime Factor Algorithm doing row FFT's then column FFT's fft_pfa :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Radix-2 Decimation in Frequency FFT module Numeric.Transform.Fourier.R2DIF -- | Radix-2 Decimation in Frequency FFT fft_r2dif :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Radix-2 Decimation in Time FFT module Numeric.Transform.Fourier.R2DIT -- | Radix-2 Decimation in Time FFT fft_r2dit :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Radix-4 Decimation in Frequency FFT module Numeric.Transform.Fourier.R4DIF -- | Radix-4 Decimation in Frequency FFT fft_r4dif :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Rader's Algorithm for computing prime length FFT's module Numeric.Transform.Fourier.Rader -- | Rader's Algorithm using direct convolution fft_rader1 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> Array a (Complex b) -- | Rader's Algorithm using FFT convolution fft_rader2 :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | FFT driver functions module Numeric.Transform.Fourier.FFT -- | This is the driver routine for calculating FFT's. All of the recursion -- in the various algorithms are defined in terms of fft. fft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -- | Inverse FFT, including scaling factor, defined in terms of fft ifft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a (Complex b) -- | This is the algorithm for computing 2N-point real FFT with an N-point -- complex FFT, defined in terms of fft rfft :: (Ix a, Integral a, RealFloat b) => Array a b -> Array a (Complex b) -- | This is the algorithm for computing a 2N-point real inverse FFT with -- an N-point complex FFT, defined in terms of ifft irfft :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> Array a b -- | Algorithm for 2 N-point real FFT's computed with N-point complex FFT, -- defined in terms of fft r2fft :: (Ix a, Integral a, RealFloat b) => Array a b -> Array a b -> (Array a (Complex b), Array a (Complex b)) -- | Utility functions based on the FFT module Numeric.Transform.Fourier.FFTUtils fft_mag :: (RealFloat b, Integral a, Ix a) => Array a (Complex b) -> Array a b fft_db :: (RealFloat b, Integral a, Ix a) => Array a (Complex b) -> Array a b fft_phase :: (Integral a, Ix a) => Array a (Complex Double) -> Array a Double fft_grd :: (Integral i, RealFloat a, Ix i) => Array i (Complex a) -> Array i a fft_info :: (Integral i, Ix i) => Array i (Complex Double) -> (Array i Double, Array i Double, Array i Double, Array i Double) rfft_mag :: (RealFloat b, Integral a, Ix a) => Array a b -> Array a b rfft_db :: (RealFloat b, Integral a, Ix a) => Array a b -> Array a b rfft_phase :: (Integral a, Ix a) => Array a Double -> Array a Double rfft_grd :: (Integral i, Ix i, RealFloat a) => Array i a -> Array i a rfft_info :: (Integral i, Ix i) => Array i Double -> (Array i Double, Array i Double, Array i Double, Array i Double) write_fft_info :: (Ix i, Integral i) => String -> Array i (Complex Double) -> IO () write_rfft_info :: String -> Array Int Double -> IO () -- | Module to perform fast linear convolution of two sequences module DSP.FastConvolution -- | fast_conv convolves two finite sequences using DFT -- relationships fast_conv :: RealFloat b => Array Int (Complex b) -> Array Int (Complex b) -> Array Int (Complex b) -- | Split-Radix Decimation in Frequency FFT module Numeric.Transform.Fourier.SRDIF -- | Split-Radix Decimation in Frequency FFT fft_srdif :: (Ix a, Integral a, RealFloat b) => Array a (Complex b) -> a -> (Array a (Complex b) -> Array a (Complex b)) -> Array a (Complex b) -- | Sliding FFT Algorithm module Numeric.Transform.Fourier.SlidingFFT -- | Sliding FFT sfft :: RealFloat a => Int -> [Complex a] -> [Array Int (Complex a)] -- | Simple module for handling polynomials. module Polynomial.Basic -- | Evaluate a polynomial using Horner's method. polyeval :: Num a => [a] -> a -> a -- | Add two polynomials polyadd :: Num a => [a] -> [a] -> [a] polyAddScalar :: Num a => a -> [a] -> [a] -- | Subtract two polynomials polysub :: Num a => [a] -> [a] -> [a] -- | Scale a polynomial polyscale :: Num a => a -> [a] -> [a] -- | Multiply two polynomials polymult :: Num a => [a] -> [a] -> [a] polymultAlt :: Num a => [a] -> [a] -> [a] -- | Divide two polynomials polydiv :: Fractional a => [a] -> [a] -> [a] -- | Modulus of two polynomials (remainder of division) polymod :: Fractional a => [a] -> [a] -> [a] -- | Raise a polynomial to a non-negative integer power polypow :: (Num a, Integral b) => [a] -> b -> [a] -- | Polynomial substitution y(n) = x(w(n)) polysubst :: Num a => [a] -> [a] -> [a] polysubstAlt :: Num a => [a] -> [a] -> [a] -- | Polynomial substitution y(n) = x(w(n)) where the coefficients -- of x are also polynomials. polyPolySubst :: Num a => [a] -> [[a]] -> [a] -- | Polynomial derivative polyderiv :: Num a => [a] -> [a] -- | Polynomial integration polyinteg :: Fractional a => [a] -> a -> [a] -- | Convert roots to a polynomial roots2poly :: Num a => [a] -> [a] -- | 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
--   
-- --
--   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. --
  3. Multiply the results by (1+z^-1)^n == [ 1, 1 ]^n
  4. --
  5. Add up all of the common terms
  6. --
  7. Normalize all of the coeficients by a0
  8. --
-- -- where n is the maximum order of the numerator and denominator module DSP.Filter.IIR.Bilinear zm :: (Integral b, Fractional a) => a -> b -> [a] zp :: (Integral b, Num a) => b -> [a] step1 :: Fractional a => a -> [a] -> [[a]] step2 :: (Num a, Integral b) => b -> [[a]] -> [[a]] step3 :: Num a => [[a]] -> [a] step4 :: Fractional a => a -> [a] -> [a] -- | Performs the bilinear transform bilinear :: Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Function for frequency prewarping prewarp :: Double -> Double -> Double -- | 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) module DSP.Filter.FIR.Smooth -- | designs smooth FIR filters smoothfir :: (Ix a, Integral a, Fractional b) => a -> a -> Array a b -- | 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. module DSP.Filter.FIR.PolyInterp -- | mkcoef takes the continuous impluse response function (one of -- the functions below, f) and number of points in the -- interpolation, p, time shifts it by x, samples it, -- and creates an array with the interpolation coeficients that can be -- used as a FIR filter. mkcoef :: (Num a, Ix b, Integral b) => (a -> a) -> b -> a -> Array b a bspline_1p0o :: (Ord a, Fractional a) => a -> a bspline_2p1o :: (Ord a, Fractional a) => a -> a bspline_4p3o :: (Ord a, Fractional a) => a -> a bspline_6p5o :: (Ord a, Fractional a) => a -> a lagrange_4p3o :: (Ord a, Fractional a) => a -> a lagrange_6p5o :: (Ord a, Fractional a) => a -> a hermite_4p3o :: (Ord a, Fractional a) => a -> a hermite_6p3o :: (Ord a, Fractional a) => a -> a hermite_6p5o :: (Ord a, Fractional a) => a -> a sndosc_4p5o :: (Ord a, Fractional a) => a -> a sndosc_6p5o :: (Ord a, Fractional a) => a -> a watte_4p2o :: (Ord a, Fractional a) => a -> a parabolic2x_4p2o :: (Ord a, Fractional a) => a -> a optimal_2p3o2x :: (Ord a, Fractional a) => a -> a optimal_2p3o4x :: (Ord a, Fractional a) => a -> a optimal_2p3o8x :: (Ord a, Fractional a) => a -> a optimal_2p3o16x :: (Ord a, Fractional a) => a -> a optimal_2p3o32x :: (Ord a, Fractional a) => a -> a optimal_4p2o2x :: (Ord a, Fractional a) => a -> a optimal_4p2o4x :: (Ord a, Fractional a) => a -> a optimal_4p2o8x :: (Ord a, Fractional a) => a -> a optimal_4p2o16x :: (Ord a, Fractional a) => a -> a optimal_4p2o32x :: (Ord a, Fractional a) => a -> a optimal_4p3o2x :: (Ord a, Fractional a) => a -> a optimal_4p3o4x :: (Ord a, Fractional a) => a -> a optimal_4p3o8x :: (Ord a, Fractional a) => a -> a optimal_4p3o16x :: (Ord a, Fractional a) => a -> a optimal_4p3o32x :: (Ord a, Fractional a) => a -> a optimal_4p4o2x :: (Ord a, Fractional a) => a -> a optimal_4p4o4x :: (Ord a, Fractional a) => a -> a optimal_4p4o8x :: (Ord a, Fractional a) => a -> a optimal_4p4o16x :: (Ord a, Fractional a) => a -> a optimal_4p4o32x :: (Ord a, Fractional a) => a -> a optimal_6p4o2x :: (Ord a, Fractional a) => a -> a optimal_6p4o4x :: (Ord a, Fractional a) => a -> a optimal_6p4o8x :: (Ord a, Fractional a) => a -> a optimal_6p4o16x :: (Ord a, Fractional a) => a -> a optimal_6p4o32x :: (Ord a, Fractional a) => a -> a optimal_6p5o2x :: (Ord a, Fractional a) => a -> a optimal_6p5o4x :: (Ord a, Fractional a) => a -> a optimal_6p5o8x :: (Ord a, Fractional a) => a -> a optimal_6p5o16x :: (Ord a, Fractional a) => a -> a optimal_6p5o32x :: (Ord a, Fractional a) => a -> a -- | Analog prototype filter transforms module DSP.Filter.Analog.Transform -- | Lowpass to lowpass: s --> s/wc a_lp2lp :: Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Lowpass to highpass: s --> wc/s a_lp2hp :: Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Lowpass to bandpass: s --> (s^2 + wl*wu) / (s(wu-wl)) a_lp2bp :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Lowpass to bandstop: s --> (s(wu-wl)) / (s^2 + wl*wu) a_lp2bs :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double]) substitute :: ([Double], [Double]) -> ([Double], [Double]) -> ([Double], [Double]) propSubstituteRecip :: ([Double], [Double]) -> ([Double], [Double]) -> Bool propSubstituteAlt :: ([Double], [Double]) -> ([Double], [Double]) -> Bool -- | Digital IIR filter transforms -- -- Reference: R&G, pg 260; O&S, pg 434; P&M, pg 699 -- -- Notation follows O&S module DSP.Filter.IIR.Transform -- | Lowpass to lowpass: z^-1 --> (z^-1 - a)/(1 - a*z^-1) d_lp2lp :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Lowpass to Highpass: z^-1 --> -(z^-1 + a)/(1 + a*z^-1) d_lp2hp :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Lowpass to Bandpass: z^-1 --> d_lp2bp :: Double -> Double -> Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Lowpass to Bandstop: z^-1 --> d_lp2bs :: Double -> Double -> Double -> ([Double], [Double]) -> ([Double], [Double]) -- | Module for generating analog filter prototypes module DSP.Filter.Analog.Prototype -- | Generates Butterworth filter prototype butterworth :: Int -> ([Double], [Double]) -- | Generates Chebyshev filter prototype chebyshev1 :: Double -> Int -> ([Double], [Double]) -- | Generates Inverse Chebyshev filter prototype chebyshev2 :: Double -> Int -> ([Double], [Double]) -- | Lowpass, Highpass, Bandpass IIR design functions -- -- Method: -- --
    --
  1. Design analog prototype
  2. --
  3. Perform analog-to-analog frequency transformation
  4. --
  5. Perform bilinear transform
  6. --
module DSP.Filter.IIR.Design poly2iir :: ([a], [b]) -> (Array Int a, Array Int b) -- | Generates lowpass Butterworth IIR filters butterworthLowpass :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) butterworthHighpass :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) butterworthBandpass :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) -- | Generates lowpass Chebyshev IIR filters chebyshev1Lowpass :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) -- | Generates lowpass Inverse Chebyshev IIR filters chebyshev2Lowpass :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) -- | Generates lowpass Butterworth IIR filters -- | Deprecated: Use butterworthLowpass instead mkButterworth :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) -- | Generates lowpass Chebyshev IIR filters -- | Deprecated: Use chebyshev1Lowpass instead mkChebyshev1 :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) -- | Generates lowpass Inverse Chebyshev IIR filters -- | Deprecated: Use chebyshev2Lowpass instead mkChebyshev2 :: (Double, Double) -> (Double, Double) -> (Array Int Double, Array Int Double) -- | Simple module for generating Chebyshev polynomials -- --
--   T_0(x) = 1
--   
-- --
--   T_1(x) = x
--   
-- --
--   T_N+1(x) = 2x T_N(x) - T_N-1(x)
--   
module Polynomial.Chebyshev -- | generates Chebyshev polynomials cheby :: (Integral a, Num b) => a -> [b] -- | Module for generating analog filter responses -- -- Formulas are from Oppenheim and Schafer, Appendix B module DSP.Filter.Analog.Response -- | Butterworth filter response function butterworth_H :: Int -> Double -> Double -> Double -- | Chebyshev filter response function chebyshev1_H :: Int -> Double -> Double -> Double -> Double -- | Inverse Chebyshev filter response function -- -- Note that w_c is a property of the stopband for this filter chebyshev2_H :: Int -> Double -> Double -> Double -> Double -- | Simple module for generating Maclaurin series representation of a few -- functions: -- --
--   f(x) = sum [ a_i * x^i | i <- [0..] ]
--   
-- -- The Int parameter for all functions is the order of -- the polynomial, eg: -- --
--   [ a_i | i <- [0..N] ]
--   
-- -- and not the number of non-zero terms module Polynomial.Maclaurin -- | e^x polyexp :: Int -> [Double] -- | ln (1+x), 0 <= x <= 1 polyln1 :: Int -> [Double] -- | cos x polycos :: Int -> [Double] -- | sin x polysin :: Int -> [Double] -- | atan x, -1 < x < 1 polyatan :: Int -> [Double] -- | cosh x polycosh :: Int -> [Double] -- | sinh x polysinh :: Int -> [Double] -- | atanh x polyatanh :: Int -> [Double] -- | Root finder using Laguerre's method module Polynomial.Roots -- | Root finder using Laguerre's method roots :: RealFloat a => a -> Int -> [Complex a] -> [Complex a] -- | Matched-z transform -- -- References: Proakis and Manolakis, Rabiner and Gold module DSP.Filter.IIR.Matchedz -- | Performs the matched-z transform matchedz :: Double -> ([Double], [Double]) -> ([Double], [Double])