----------------------------------------------------------------------------- -- | -- Module : DSP.Filter.Analog.Prototype -- Copyright : (c) Matthew Donadio 2003 -- License : GPL -- -- Maintainer : m.p.donadio@ieee.org -- Stability : experimental -- Portability : portable -- -- Module for generating analog filter prototypes -- ----------------------------------------------------------------------------- -- Notes (mainly for self): -- The gain of an analog filter is -- gain = abs \$ realPart \$ product zeros / product poles -- = abs \$ b_m / a_n -- For a Butterworth filter, the product of the poles is one, so we don't -- have to worry about any gain. -- For a Chebyshev 1 filter, the product of the poles is a_n, which is -- the head of the polynomial. We make this b_0 to set the gain in the -- passband. -- For a Chebyshev 2 filter, we use the full gain formula because we want -- to set the gain to unity at DC. -- TODO: Do we want to include Bessel filters? module DSP.Filter.Analog.Prototype where import Data.Complex (Complex((:+)), realPart) import Polynomial.Basic (roots2poly) -- | Generates Butterworth filter prototype butterworth :: Int -- ^ N -> ([Double],[Double]) -- ^ (b,a) butterworth n = (num, den) 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 ] den = map realPart \$ roots2poly \$ poles -- | Generates Chebyshev filter prototype chebyshev1 :: Double -- ^ epsilon -> Int -- ^ N -> ([Double],[Double]) -- ^ (b,a) chebyshev1 eps n = (num, den) where poles = [ (-u k) :+ (w k) | k <- [0..(n-1)] ] u k = sinh v0 * sin (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) w k = cosh v0 * cos (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) num = [ gain ] den = map realPart \$ roots2poly \$ poles v0 = asinh (1/eps) / fromIntegral n gain = if even n then abs \$ head den / sqrt (1 + eps^(2::Int)) else abs \$ head den -- | Generates Inverse Chebyshev filter prototype chebyshev2 :: Double -- ^ epsilon -> Int -- ^ N -> ([Double],[Double]) -- ^ (b,a) chebyshev2 eps n = (num, den) where zeros = [ 0 :+ 1 / wz k | k <- [0..(n-1)], 2*k+1 /= n ] poles = [ 1 / ((-u k) :+ (w k)) | k <- [0..(n-1)] ] wz k = cos (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) u k = sinh v0 * sin (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) w k = cosh v0 * cos (fromIntegral (2*k+1) * pi / fromIntegral (2*n)) num = map (*gain) \$ map realPart \$ roots2poly \$ zeros den = map realPart \$ roots2poly \$ poles v0 = asinh (1/eps) / fromIntegral n gain = abs \$ realPart \$ product poles / product zeros