-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.IIR.IIR
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- 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@
--
-----------------------------------------------------------------------------

-- TODO: Should these use Arrays for a and b?  Tuples?

{-

Reference:

@Book{dsp,
  author = 	 "Alan V. Oppenheim and Ronald W. Schafer",
  title = 	 "Discrete-Time Signal Processing",
  publisher = 	 "Prentice-Hall",
  year = 	 1989,
  address =	 "Englewood Cliffs",
  series =       {Prentice-Hall Signal Processing Series}
}

However, we differ in the convention of the sign of the poles, as
noted in the module header.

-}

module DSP.Filter.IIR.IIR (integrator,
    fos_df1, fos_df2, fos_df2t,
    biquad_df1, biquad_df2, biquad_df2t,
    iir_df1, iir_df2,
    -- for testing
    xt, yt, f1, f2, f3, f4, f5,
    ) where

import Data.Array

import DSP.Filter.FIR.FIR

-- | This is an integrator when a==1, and a leaky integrator when @0 \< a \< 1@.
--
--  @y[n] = a * y[n-1] + x[n]@

{-# specialize integrator :: Float -> [Float] -> [Float] #-}
{-# specialize integrator :: Double -> [Double] -> [Double] #-}

integrator :: Num a => a -- ^ a
	   -> [a] -- ^ x[n]
	   -> [a] -- ^ y[n]

integrator :: forall a. Num a => a -> [a] -> [a]
integrator a
a [a]
x = forall a. Num a => a -> a -> [a] -> [a]
integrator' a
a a
0 [a]
x


integrator' :: Num a => a -> a-> [a] -> [a]
integrator' :: forall a. Num a => a -> a -> [a] -> [a]
integrator' a
_ a
_  []     = []
integrator' a
a a
y1 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => a -> a -> [a] -> [a]
integrator' a
a a
y [a]
xs
    where y :: a
y = a
a forall a. Num a => a -> a -> a
* a
y1 forall a. Num a => a -> a -> a
+ a
x

-- | First order section, DF1
--
--	@v[n] = b0 * x[n] + b1 * x[n-1]@
--
--	@y[n] = v[n] - a1 * y[n-1]@

{-# specialize fos_df1 :: Float -> Float -> Float -> [Float] -> [Float] #-}
{-# specialize fos_df1 :: Double -> Double -> Double -> [Double] -> [Double] #-}

fos_df1 :: Num a => a -- ^ a_1
	-> a -- ^ b_0
	-> a -- ^ b_1
	-> [a] -- ^ x[n]
	-> [a] -- ^ y[n]

fos_df1 :: forall a. Num a => a -> a -> a -> [a] -> [a]
fos_df1 a
a1 a
b0 a
b1 [a]
x = forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
fos_df1' a
a1 a
b0 a
b1 a
0 a
0 [a]
x

fos_df1' :: Num a => a -> a -> a -> a -> a -> [a] -> [a]
fos_df1' :: forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
fos_df1' a
_  a
_  a
_  a
_  a
_  []     = []
fos_df1' a
a1 a
b0 a
b1 a
x1 a
y1 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
fos_df1' a
a1 a
b0 a
b1 a
x a
y [a]
xs
    where v :: a
v = a
b0 forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ a
b1 forall a. Num a => a -> a -> a
* a
x1
	  y :: a
y = a
v      forall a. Num a => a -> a -> a
- a
a1 forall a. Num a => a -> a -> a
* a
y1


-- | First order section, DF2
--
--	@w[n] = -a1 * w[n-1] + x[n]@
--
--	@y[n] = b0 * w[n] + b1 * w[n-1]@

{-# specialize fos_df2 :: Float -> Float -> Float -> [Float] -> [Float] #-}
{-# specialize fos_df2 :: Double -> Double -> Double -> [Double] -> [Double] #-}

fos_df2 :: Num a => a -- ^ a_1
	-> a -- ^ b_0
	-> a -- ^ b_1
	-> [a] -- ^ x[n]
	-> [a] -- ^ y[n]

fos_df2 :: forall a. Num a => a -> a -> a -> [a] -> [a]
fos_df2 a
a1 a
b0 a
b1 [a]
x = forall a. Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2' a
a1 a
b0 a
b1 a
0 [a]
x

fos_df2' :: Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2' :: forall a. Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2' a
_  a
_  a
_  a
_  []     = []
fos_df2' a
a1 a
b0 a
b1 a
w1 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2' a
a1 a
b0 a
b1 a
w [a]
xs
    where w :: a
w = a
x forall a. Num a => a -> a -> a
- a
a1 forall a. Num a => a -> a -> a
* a
w1
          y :: a
y = a
b0 forall a. Num a => a -> a -> a
* a
w forall a. Num a => a -> a -> a
+ a
b1 forall a. Num a => a -> a -> a
* a
w1

-- | First order section, DF2T
--
--	@v0[n] = b0 * x[n] + v1[n-1]@
--
--	@y[n] = v0[n]@
--
--	@v1[n] = -a1 * y[n] + b1 * x[n]@

{-# specialize fos_df2t :: Float -> Float -> Float -> [Float] -> [Float] #-}
{-# specialize fos_df2t :: Double -> Double -> Double -> [Double] -> [Double] #-}

fos_df2t :: Num a => a -- ^ a_1
	    -> a -- ^ b_0
	    -> a -- ^ b_1
	    -> [a] -- ^ x[n]
	    -> [a] -- ^ y[n]

fos_df2t :: forall a. Num a => a -> a -> a -> [a] -> [a]
fos_df2t a
a1 a
b0 a
b1 [a]
x = forall a. Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2t' a
a1 a
b0 a
b1 a
0 [a]
x

fos_df2t' :: Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2t' :: forall a. Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2t' a
_  a
_  a
_  a
_   []     = []
fos_df2t' a
a1 a
b0 a
b1 a
v11 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => a -> a -> a -> a -> [a] -> [a]
fos_df2t' a
a1 a
b0 a
b1 a
v1 [a]
xs
    where v0 :: a
v0 = a
b0 forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ a
v11
          y :: a
y  = a
v0
	  v1 :: a
v1 = -a
a1 forall a. Num a => a -> a -> a
* a
y forall a. Num a => a -> a -> a
+ a
b1 forall a. Num a => a -> a -> a
* a
x

-- | 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]@

{-# specialize biquad_df1 :: Float -> Float -> Float -> Float -> Float -> [Float] -> [Float] #-}
{-# specialize biquad_df1 :: Double -> Double -> Double -> Double -> Double -> [Double] -> [Double] #-}

biquad_df1 :: Num a => a -- ^ a_1
	   -> a -- ^ a_2
	   -> a -- ^ b_0
	   -> a -- ^ b_1
	   -> a -- ^ b_2
	   -> [a] -- ^ x[n]
	   -> [a] -- ^ y[n]

biquad_df1 :: forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
biquad_df1 a
a1 a
a2 a
b0 a
b1 a
b2 [a]
x = forall a.
Num a =>
a -> a -> a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df1 a
a1 a
a2 a
b0 a
b1 a
b2 a
0 a
0 a
0 a
0 [a]
x

df1 :: Num a => a -> a -> a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df1 :: forall a.
Num a =>
a -> a -> a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df1 a
_  a
_  a
_  a
_  a
_  a
_  a
_  a
_  a
_  []     = []
df1 a
a1 a
a2 a
b0 a
b1 a
b2 a
x1 a
x2 a
y1 a
y2 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a.
Num a =>
a -> a -> a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df1 a
a1 a
a2 a
b0 a
b1 a
b2 a
x a
x1 a
y a
y1 [a]
xs
    where v :: a
v = a
b0 forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ a
b1 forall a. Num a => a -> a -> a
* a
x1 forall a. Num a => a -> a -> a
+ a
b2 forall a. Num a => a -> a -> a
* a
x2
	  y :: a
y = a
v      forall a. Num a => a -> a -> a
- a
a1 forall a. Num a => a -> a -> a
* a
y1 forall a. Num a => a -> a -> a
- a
a2 forall a. Num a => a -> a -> a
* a
y2


-- | 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]@

{-# specialize biquad_df2 :: Float -> Float -> Float -> Float -> Float -> [Float] -> [Float] #-}
{-# specialize biquad_df2 :: Double -> Double -> Double -> Double -> Double -> [Double] -> [Double] #-}

biquad_df2 :: Num a => a -- ^ a_1
	   -> a -- ^ a_2
	   -> a -- ^ b_0
	   -> a -- ^ b_1
	   -> a -- ^ b_2
	   -> [a] -- ^ x[n]
	   -> [a] -- ^ y[n]

biquad_df2 :: forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
biquad_df2 a
a1 a
a2 a
b0 a
b1 a
b2 [a]
x = forall a. Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2 a
a1 a
a2 a
b0 a
b1 a
b2 a
0 a
0 [a]
x

df2 :: Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2 :: forall a. Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2 a
_  a
_  a
_  a
_  a
_  a
_  a
_  []     = []
df2 a
a1 a
a2 a
b0 a
b1 a
b2 a
w1 a
w2 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2 a
a1 a
a2 a
b0 a
b1 a
b2 a
w a
w1 [a]
xs
    where w :: a
w = a
x forall a. Num a => a -> a -> a
- a
a1 forall a. Num a => a -> a -> a
* a
w1 forall a. Num a => a -> a -> a
- a
a2 forall a. Num a => a -> a -> a
* a
w2
          y :: a
y = a
b0 forall a. Num a => a -> a -> a
* a
w forall a. Num a => a -> a -> a
+ a
b1 forall a. Num a => a -> a -> a
* a
w1 forall a. Num a => a -> a -> a
+ a
b2 forall a. Num a => a -> a -> a
* a
w2

-- | 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]@

{-# specialize biquad_df2t :: Float -> Float -> Float -> Float -> Float -> [Float] -> [Float] #-}
{-# specialize biquad_df2t :: Double -> Double -> Double -> Double -> Double -> [Double] -> [Double] #-}

biquad_df2t :: Num a => a -- ^ a_1
	    -> a -- ^ a_2
	    -> a -- ^ b_0
	    -> a -- ^ b_1
	    -> a -- ^ b_2
	    -> [a] -- ^ x[n]
	    -> [a] -- ^ y[n]

biquad_df2t :: forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
biquad_df2t a
a1 a
a2 a
b0 a
b1 a
b2 [a]
x = forall a. Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2t a
a1 a
a2 a
b0 a
b1 a
b2 a
0 a
0 [a]
x

df2t :: Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2t :: forall a. Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2t a
_  a
_  a
_  a
_  a
_  a
_   a
_   []     = []
df2t a
a1 a
a2 a
b0 a
b1 a
b2 a
v11 a
v21 (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => a -> a -> a -> a -> a -> a -> a -> [a] -> [a]
df2t a
a1 a
a2 a
b0 a
b1 a
b2 a
v1 a
v2 [a]
xs
    where v0 :: a
v0 = a
b0 forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ a
v11
          y :: a
y = a
v0
	  v1 :: a
v1 = -a
a1 forall a. Num a => a -> a -> a
* a
y forall a. Num a => a -> a -> a
+ a
b1 forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ a
v21
	  v2 :: a
v2 = -a
a2 forall a. Num a => a -> a -> a
* a
y forall a. Num a => a -> a -> a
+ a
b2 forall a. Num a => a -> a -> a
* a
x

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

{- specialize iir_df1 :: (Array Int Float, Array Int Float) -> [Float] -> [Float] -}
{- specialize iir_df1 :: (Array Int Double, Array Int Double) -> [Double] -> [Double] -}

iir_df1 :: (Num a, Eq a) => (Array Int a, Array Int a) -- ^ (b,a)
	-> [a] -- ^ x[n]
	-> [a] -- ^ y[n]

iir_df1 :: forall a. (Num a, Eq a) => (Array Int a, Array Int a) -> [a] -> [a]
iir_df1 (Array Int a
b,Array Int a
a) [a]
x = [a]
y
    where v :: [a]
v = forall a. (Num a, Eq a) => Array Int a -> [a] -> [a]
fir Array Int a
b [a]
x
	  y :: [a]
y = forall a. Num a => Array Int a -> Array Int a -> [a] -> [a]
iir'df1 Array Int a
a Array Int a
w [a]
v
	  w :: Array Int a
w = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
n) forall a b. (a -> b) -> a -> b
$ forall a. a -> [a]
repeat a
0
	  n :: Int
n = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int a
a

{- specialize iir'df1 :: Array Int Float -> Array Int Float -> [Float] -> [Float] -}
{- specialize iir'df1 :: Array Int Double -> Array Int Double -> [Double] -> [Double] -}

iir'df1 :: (Num a) => Array Int a -> Array Int a -> [a] -> [a]
iir'df1 :: forall a. Num a => Array Int a -> Array Int a -> [a] -> [a]
iir'df1 Array Int a
_ Array Int a
_ []  = []
iir'df1 Array Int a
a Array Int a
w (a
v:[a]
vs) = a
y forall a. a -> [a] -> [a]
: forall a. Num a => Array Int a -> Array Int a -> [a] -> [a]
iir'df1 Array Int a
a Array Int a
w' [a]
vs
    where y :: a
y  = a
v forall a. Num a => a -> a -> a
- forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array Int a
aforall i e. Ix i => Array i e -> i -> e
!Int
i forall a. Num a => a -> a -> a
* Array Int a
wforall i e. Ix i => Array i e -> i -> e
!Int
i | Int
i <- [Int
1..Int
n] ]
          w' :: Array Int a
w' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
n) forall a b. (a -> b) -> a -> b
$ a
y forall a. a -> [a] -> [a]
: forall i e. Array i e -> [e]
elems Array Int a
w
	  n :: Int
n  = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int 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]@

{- specialize iir_df2 :: (Array Int Float, Array Int Float) -> [Float] -> [Float] -}
{- specialize iir_df2 :: (Array Int Double, Array Int Double) -> [Double] -> [Double] -}

iir_df2 :: (Num a) => (Array Int a, Array Int a) -- ^ (b,a)
	-> [a] -- ^ x[n]
	-> [a] -- ^ y[n]

iir_df2 :: forall a. Num a => (Array Int a, Array Int a) -> [a] -> [a]
iir_df2 (Array Int a
b,Array Int a
a) [a]
x = [a]
y
    where y :: [a]
y = forall a.
Num a =>
(Array Int a, Array Int a) -> Array Int a -> [a] -> [a]
iir'df2 (Array Int a
b,Array Int a
a) Array Int a
w [a]
x
	  w :: Array Int a
w = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
mn) forall a b. (a -> b) -> a -> b
$ forall a. a -> [a]
repeat a
0
	  m :: Int
m = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int a
b
	  n :: Int
n = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int a
a
	  mn :: Int
mn = forall a. Ord a => a -> a -> a
max Int
m Int
n

{- specialize iir'df2 :: Array Int Float -> Array Int Float -> [Float] -> [Float] -}
{- specialize iir'df2 :: Array Int Double -> Array Int Double -> [Double] -> [Double] -}

iir'df2 :: (Num a) => (Array Int a,Array Int a) -> Array Int a -> [a] -> [a]
iir'df2 :: forall a.
Num a =>
(Array Int a, Array Int a) -> Array Int a -> [a] -> [a]
iir'df2 (Array Int a, Array Int a)
_     Array Int a
_ []     = []
iir'df2 (Array Int a
b,Array Int a
a) Array Int a
w (a
x:[a]
xs) = a
y forall a. a -> [a] -> [a]
: forall a.
Num a =>
(Array Int a, Array Int a) -> Array Int a -> [a] -> [a]
iir'df2 (Array Int a
b,Array Int a
a) Array Int a
w' [a]
xs
    where y :: a
y  = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array Int a
bforall i e. Ix i => Array i e -> i -> e
!Int
i forall a. Num a => a -> a -> a
* Array Int a
w'forall i e. Ix i => Array i e -> i -> e
!Int
i | Int
i <- [Int
0..Int
m] ]
          w0 :: a
w0 = a
x forall a. Num a => a -> a -> a
- forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [ Array Int a
aforall i e. Ix i => Array i e -> i -> e
!Int
i forall a. Num a => a -> a -> a
* Array Int a
w'forall i e. Ix i => Array i e -> i -> e
!Int
i | Int
i <- [Int
1..Int
m] ]
	  w' :: Array Int a
w' = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
mn) forall a b. (a -> b) -> a -> b
$ a
w0 forall a. a -> [a] -> [a]
: forall i e. Array i e -> [e]
elems Array Int a
w
	  m :: Int
m  = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int a
b
	  mn :: Int
mn = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall i e. Array i e -> (i, i)
bounds Array Int a
w

---------

-- test

xt :: [Double]
xt :: [Double]
xt = [ Double
1, Double
0, Double
0, Double
0, Double
0, Double
0, Double
0, Double
0 ] :: [Double]

yt :: [Double]
yt :: [Double]
yt = forall a. Num a => a -> [a] -> [a]
integrator Double
0.5 [Double]
xt

f1 :: Fractional a => [a] -> [a]
f1 :: forall a. Fractional a => [a] -> [a]
f1 [a]
x = forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
biquad_df1  (-a
0.4) a
0.3 a
0.5 a
0.4 (-a
0.3) [a]
x

f2 :: Fractional a => [a] -> [a]
f2 :: forall a. Fractional a => [a] -> [a]
f2 [a]
x = forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
biquad_df2  (-a
0.4) a
0.3 a
0.5 a
0.4 (-a
0.3) [a]
x

f3 :: Fractional a => [a] -> [a]
f3 :: forall a. Fractional a => [a] -> [a]
f3 [a]
x = forall a. Num a => a -> a -> a -> a -> a -> [a] -> [a]
biquad_df2t (-a
0.4) a
0.3 a
0.5 a
0.4 (-a
0.3) [a]
x

at :: Array Int Double
at :: Array Int Double
at = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
1,Int
2) [ -Double
0.4, Double
0.3 ]

bt :: Array Int Double
bt :: Array Int Double
bt = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0,Int
2) [ Double
0.5, Double
0.4, -Double
0.3 ]

f4 :: [Double] -> [Double]
f4 :: [Double] -> [Double]
f4 [Double]
x = forall a. (Num a, Eq a) => (Array Int a, Array Int a) -> [a] -> [a]
iir_df1 (Array Int Double
bt,Array Int Double
at) [Double]
x

f5 :: [Double] -> [Double]
f5 :: [Double] -> [Double]
f5 [Double]
x = forall a. Num a => (Array Int a, Array Int a) -> [a] -> [a]
iir_df2 (Array Int Double
bt,Array Int Double
at) [Double]
x