-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.IIR.Transform
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Digital IIR filter transforms
--
-- Reference: R&G, pg 260; O&S, pg 434; P&M, pg 699
--
-- Notation follows O&S
--
-----------------------------------------------------------------------------

-- TODO: These need more testing.  I checked the lp2hp case against O&S
-- which verifies substitute and lp2hp,nd I triple checked the parameters
-- for the others.  I need to find test vectors for the other cases for
-- proper testing, though.

module DSP.Filter.IIR.Transform (d_lp2lp, d_lp2hp, d_lp2bp, d_lp2bs) where

import DSP.Filter.Analog.Transform (substitute)
import Numeric.Special.Trigonometric (cot)


-- | 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')

d_lp2lp :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
d_lp2lp Double
tp Double
wp ([Double]
num,[Double]
den) = ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den)
    where nsub :: [Double]
nsub = [-Double
a, Double
1]
          dsub :: [Double]
dsub = [Double
1, -Double
a]
          a :: Double
a = forall a. Floating a => a -> a
sin ((Double
tpforall a. Num a => a -> a -> a
-Double
wp)forall a. Fractional a => a -> a -> a
/Double
2) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sin ((Double
tpforall a. Num a => a -> a -> a
+Double
wp)forall a. Fractional a => a -> a -> a
/Double
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')

d_lp2hp :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
d_lp2hp Double
tp Double
wp ([Double]
num,[Double]
den) = ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den)
    where nsub :: [Double]
nsub = [Double
a, Double
1]
          dsub :: [Double]
dsub = [-Double
1, -Double
a]
          a :: Double
a = -forall a. Floating a => a -> a
cos ((Double
tpforall a. Num a => a -> a -> a
+Double
wp)forall a. Fractional a => a -> a -> a
/Double
2) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
cos ((Double
tpforall a. Num a => a -> a -> a
-Double
wp)forall a. Fractional a => a -> a -> a
/Double
2)

-- | Lowpass to Bandpass: z^-1 -->

d_lp2bp :: Double -- ^ theta_p
        -> Double -- ^ omega_p1
        -> Double -- ^ omega_p2
        -> ([Double], [Double]) -- ^ (b,a)
        -> ([Double], [Double]) -- ^ (b',a')

d_lp2bp :: Double
-> Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
d_lp2bp Double
tp Double
wp1 Double
wp2 ([Double]
num,[Double]
den) = ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den)
    where nsub :: [Double]
nsub = [ (Double
kforall a. Num a => a -> a -> a
-Double
1)forall a. Fractional a => a -> a -> a
/(Double
kforall a. Num a => a -> a -> a
+Double
1), -Double
2forall a. Num a => a -> a -> a
*Double
aforall a. Num a => a -> a -> a
*Double
kforall a. Fractional a => a -> a -> a
/(Double
kforall a. Num a => a -> a -> a
+Double
1), Double
1 ]
          dsub :: [Double]
dsub = [ Double
1, -Double
2forall a. Num a => a -> a -> a
*Double
aforall a. Num a => a -> a -> a
*Double
kforall a. Fractional a => a -> a -> a
/(Double
kforall a. Num a => a -> a -> a
+Double
1), (Double
kforall a. Num a => a -> a -> a
-Double
1)forall a. Fractional a => a -> a -> a
/(Double
kforall a. Num a => a -> a -> a
+Double
1) ]
          a :: Double
a = forall a. Floating a => a -> a
cos ((Double
wp2forall a. Num a => a -> a -> a
+Double
wp1)forall a. Fractional a => a -> a -> a
/Double
2) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
cos ((Double
wp2forall a. Num a => a -> a -> a
-Double
wp1)forall a. Fractional a => a -> a -> a
/Double
2)
          k :: Double
k = forall a. Floating a => a -> a
cot ((Double
wp2forall a. Num a => a -> a -> a
-Double
wp1)forall a. Fractional a => a -> a -> a
/Double
2) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
tan (Double
tpforall a. Fractional a => a -> a -> a
/Double
2)

-- | Lowpass to Bandstop: z^-1 -->

d_lp2bs :: Double -- ^ theta_p
        -> Double -- ^ omega_p1
        -> Double -- ^ omega_p2
        -> ([Double], [Double]) -- ^ (b,a)
        -> ([Double], [Double]) -- ^ (b',a')

d_lp2bs :: Double
-> Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
d_lp2bs Double
tp Double
wp1 Double
wp2 ([Double]
num,[Double]
den) = ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den)
    where nsub :: [Double]
nsub = [ (Double
1forall a. Num a => a -> a -> a
-Double
k)forall a. Fractional a => a -> a -> a
/(Double
1forall a. Num a => a -> a -> a
+Double
k), -Double
2forall a. Num a => a -> a -> a
*Double
aforall a. Fractional a => a -> a -> a
/(Double
1forall a. Num a => a -> a -> a
+Double
k), Double
1 ]
          dsub :: [Double]
dsub = [ Double
1, -Double
2forall a. Num a => a -> a -> a
*Double
aforall a. Fractional a => a -> a -> a
/(Double
1forall a. Num a => a -> a -> a
+Double
k), (Double
1forall a. Num a => a -> a -> a
-Double
k)forall a. Fractional a => a -> a -> a
/(Double
1forall a. Num a => a -> a -> a
+Double
k) ]
          a :: Double
a = forall a. Floating a => a -> a
cos ((Double
wp2forall a. Num a => a -> a -> a
+Double
wp1)forall a. Fractional a => a -> a -> a
/Double
2) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
cos ((Double
wp2forall a. Num a => a -> a -> a
-Double
wp1)forall a. Fractional a => a -> a -> a
/Double
2)
          k :: Double
k = forall a. Floating a => a -> a
cot ((Double
wp2forall a. Num a => a -> a -> a
-Double
wp1)forall a. Fractional a => a -> a -> a
/Double
2) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
tan (Double
tpforall a. Fractional a => a -> a -> a
/Double
2)

{-

Test vectors

O&S, pg 435

 num = polypow  [ 0.001836, 0.001836 ] 4
 den = polymult [ 0.6493, -1.5548, 1 ] [ 0.8482, -1.4996, 1 ]

-}