module DSP.Filter.IIR.Design where
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
m = length b 1
n = length a 1
mkButterworth :: (Double, Double)
-> (Double, Double)
-> (Array Int Double, Array Int Double)
mkButterworth (wp,dp) (ws,ds) = poly2iir $
bilinear 1 $
a_lp2lp wc $
butterworth n
where n = ceiling $ log (((1/ds)^!21) / ((1/(1dp))^!21)) / 2 / log (ws' / wp')
wc = ws' / ((1/ds)^!21)**(1/2/fromIntegral n)
wp' = prewarp wp 1
ws' = prewarp ws 1
mkChebyshev1 :: (Double, Double)
-> (Double, Double)
-> (Array Int Double, Array Int Double)
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)
k = wp' / ws'
n = ceiling $ acosh (1/k1) / log ((1 + sqrt (1 k^!2)) / k)
mkChebyshev2 :: (Double, Double)
-> (Double, Double)
-> (Array Int Double, Array Int Double)
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)
g = 1 dp
n = ceiling $ acosh (g / eps / sqrt (1 g^!2)) / acosh (ws' / wp')