-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.IIR.Bilinear
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- 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.  Multiply the results by (1+z^-1)^n == [ 1, 1 ]^n
--
-- 3.  Add up all of the common terms
--
-- 4.  Normalize all of the coeficients by a0
--
-- where n is the maximum order of the numerator and denominator
--
-----------------------------------------------------------------------------

-- TODO: Rework to replace roots2poly using the fact that most poles
-- and\/or zeros are either complex conjugate pairs, or real only.

-- TODO: Do we want to include prewarping?

module DSP.Filter.IIR.Bilinear {- (bilinear, prewarp) -} where

import Polynomial.Basic

-- Computes (2\/ts * (1-z^-1))^n == [ -2\/ts, 2\/ts ]^n

zm :: (Integral b, Fractional a) => a -> b -> [a]
zm :: forall b a. (Integral b, Fractional a) => a -> b -> [a]
zm a
ts b
n = forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [ -a
2forall a. Fractional a => a -> a -> a
/a
ts, a
2forall a. Fractional a => a -> a -> a
/a
ts ] b
n

-- Computes (1+z^-1)^n == [ 1, 1 ]^n

zp :: (Integral b, Num a) => b -> [a]
zp :: forall b a. (Integral b, Num a) => b -> [a]
zp b
n = forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [ a
1, a
1 ] b
n

-- Step 1: Substitute s^n with (2\/ts * (1-z^-1))^n == [ -2\/ts, 2\/ts ]^n
-- in num and den

step1 :: Fractional a => a -> [a] -> [[a]]
step1 :: forall a. Fractional a => a -> [a] -> [[a]]
step1 a
ts = forall {t}. Integral t => t -> [a] -> [[a]]
step1' (Int
0::Int)
    where step1' :: t -> [a] -> [[a]]
step1' t
_ []     = []
          step1' t
n (a
x:[a]
xs) = forall a. Num a => a -> [a] -> [a]
polyscale a
x (forall b a. (Integral b, Fractional a) => a -> b -> [a]
zm a
ts t
n) forall a. a -> [a] -> [a]
: t -> [a] -> [[a]]
step1' (t
nforall a. Num a => a -> a -> a
+t
1) [a]
xs

-- Step 2: Multiply the num and den by (1+z^-1)^n == [ 1, 1 ]^n

step2 :: (Num a, Integral b) => b -> [[a]] -> [[a]]
step2 :: forall a b. (Num a, Integral b) => b -> [[a]] -> [[a]]
step2 b
_ []     = []
step2 b
n ([a]
x:[[a]]
xs) = forall a. Num a => [a] -> [a] -> [a]
polymult (forall b a. (Integral b, Num a) => b -> [a]
zp b
n) [a]
x forall a. a -> [a] -> [a]
: forall a b. (Num a, Integral b) => b -> [[a]] -> [[a]]
step2 (b
nforall a. Num a => a -> a -> a
-b
1) [[a]]
xs

-- Step 3: Add up all of the common terms

step3 :: Num a => [[a]] -> [a]
step3 :: forall a. Num a => [[a]] -> [a]
step3 [[a]]
x = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr forall a. Num a => [a] -> [a] -> [a]
polyadd [a
0] [[a]]
x

-- Step 4: Normalize all of the coeficients by a0

step4 :: Fractional a => a -> [a] -> [a]
step4 :: forall a. Fractional a => a -> [a] -> [a]
step4 a
a0 [a]
x = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/a
a0) [a]
x

-- Glue it all together

-- | Performs the bilinear transform

bilinear :: Double -- ^ T_s
	 -> ([Double],[Double]) -- ^ (b,a)
	 -> ([Double],[Double]) -- ^ (b',a')

bilinear :: Double -> ([Double], [Double]) -> ([Double], [Double])
bilinear Double
ts ([Double]
num,[Double]
den) = ([Double]
num'', [Double]
den'')
    where n :: Int
n = forall a. Ord a => a -> a -> a
max (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
num forall a. Num a => a -> a -> a
- Int
1) (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
den forall a. Num a => a -> a -> a
- Int
1)
	  num' :: [Double]
num' = forall a. Num a => [[a]] -> [a]
step3 forall a b. (a -> b) -> a -> b
$ forall a b. (Num a, Integral b) => b -> [[a]] -> [[a]]
step2 Int
n forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => a -> [a] -> [[a]]
step1 Double
ts forall a b. (a -> b) -> a -> b
$ [Double]
num
          den' :: [Double]
den' = forall a. Num a => [[a]] -> [a]
step3 forall a b. (a -> b) -> a -> b
$ forall a b. (Num a, Integral b) => b -> [[a]] -> [[a]]
step2 Int
n forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => a -> [a] -> [[a]]
step1 Double
ts forall a b. (a -> b) -> a -> b
$ [Double]
den
          a0 :: Double
a0 = forall a. [a] -> a
last [Double]
den'
	  num'' :: [Double]
num'' = forall a. Fractional a => a -> [a] -> [a]
step4 Double
a0 [Double]
num'
	  den'' :: [Double]
den'' = forall a. Fractional a => a -> [a] -> [a]
step4 Double
a0 [Double]
den'

-- | Function for frequency prewarping

prewarp :: Double -- ^ w_c
	-> Double -- ^ T_s
	-> Double -- ^ W_c

prewarp :: Double -> Double -> Double
prewarp Double
wc Double
ts = Double
2forall a. Fractional a => a -> a -> a
/Double
ts forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
tan (Double
wc forall a. Fractional a => a -> a -> a
/ Double
2)

{-

-- Test, section 6.5.1 from Lyon's book

num1 = [ 17410.145 ]
den1 = [ 17410.145, 137.94536, 1 ]

(num1',den1') = bilinear 0.01 (num1,den1)

-- Test, from O&S, p 421

num2 = [ 0.202238 ]
den2 = polymult (polymult [ 0.5871, 0.3996, 1 ] [ 0.5871, 1.0836, 1 ] ) [ 0.5871, 1.4802, 1 ]

(num2',den2') = bilinear 1 (num2,den2)

bilinear ([0, 0, 0, 0, 1], reverse [ 1, 158881.5000000000000000000000, 6734684542.320000000000000000, 33433292062222.63200000000000, 26749649944094120.95199999999, 5301498365227355432.219999999, 308666240537082938598.7999999 ]) 48000

> num3 = [ 0, 0, 0, 0, 72687672654.5 ]
> den3 = reverse [ 1, 158881.5000000000000000000000, 6734684542.320000000000000000, 33433292062222.63200000000000, 26749649944094120.95199999999, 5301498365227355432.219999999, 308666240537082938598.7999999 ]

num31 = [ 0.0, 519.2365 ]
den31 = polypow [ 129.4, 1.0 ] 2

num32 = [ 0.0, 519.2365 ]
den32 = [ 676.7, 1.0 ]

num33 = [ 0.0, 519.2365 ]
den33 = [ 4636.0, 1.0 ]

num34 = [ 0.0, 519.2365 ]
den34 = polypow [ 76655.0, 1.0 ] 2

-}