-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.IIR.Design
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Lowpass, Highpass, Bandpass IIR design functions
--
-- Method:
--
-- (1) Design analog prototype
--
-- 2.  Perform analog-to-analog frequency transformation
--
-- 3.  Perform bilinear transform
--
-----------------------------------------------------------------------------

module DSP.Filter.IIR.Design (
   poly2iir,
   butterworthLowpass,  butterworthHighpass, butterworthBandpass,
   chebyshev1Lowpass, chebyshev2Lowpass,
   mkButterworth, mkChebyshev1, mkChebyshev2,
   ) where

import qualified DSP.Filter.Analog.Prototype as Analog
import DSP.Filter.Analog.Transform (a_lp2lp, a_lp2hp, a_lp2bp)
import DSP.Filter.IIR.Bilinear (bilinear, prewarp)

import DSP.Basic ((^!))

import Data.Array (Array, listArray)


poly2iir :: ([a], [b]) -> (Array Int a, Array Int b)
poly2iir :: forall a b. ([a], [b]) -> (Array Int a, Array Int b)
poly2iir ([a]
a,[b]
b) =
   let toArray :: [e] -> Array Int e
toArray [e]
x = forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int
0, forall (t :: * -> *) a. Foldable t => t a -> Int
length [e]
x forall a. Num a => a -> a -> a
- Int
1) forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse [e]
x
   in  (forall {e}. [e] -> Array Int e
toArray [a]
a, forall {e}. [e] -> Array Int e
toArray [b]
b)


-- | Generates lowpass Butterworth IIR filters

butterworthLowpass, mkButterworth ::
      (Double, Double) -- ^ (wp,dp)
   -> (Double, Double) -- ^ (ws,ds)
   -> (Array Int Double, Array Int Double) -- ^ (b,a)
butterworthLowpass :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
butterworthLowpass (Double, Double)
p (Double, Double)
s =
   let (Int
n, (Double, Double)
s') = (Double, Double) -> (Double, Double) -> (Int, (Double, Double))
butterworthParams (Double, Double)
p (Double, Double)
s
   in  (([Double], [Double]) -> ([Double], [Double]))
-> Int -> (Array Int Double, Array Int Double)
butterworth (Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2lp forall a b. (a -> b) -> a -> b
$ forall a. Floating a => Int -> (a, a) -> a
wc Int
n (Double, Double)
s') Int
n

butterworthHighpass, butterworthBandpass ::
   (Double, Double) ->
   (Double, Double) ->
   (Array Int Double, Array Int Double)
butterworthHighpass :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
butterworthHighpass (Double, Double)
p (Double, Double)
s =
   let (Int
n, (Double, Double)
s') = (Double, Double) -> (Double, Double) -> (Int, (Double, Double))
butterworthParams (Double, Double)
p (Double, Double)
s
   in  (([Double], [Double]) -> ([Double], [Double]))
-> Int -> (Array Int Double, Array Int Double)
butterworth (Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2hp forall a b. (a -> b) -> a -> b
$ forall a. Floating a => Int -> (a, a) -> a
wc Int
n (Double, Double)
s') Int
n

butterworthBandpass :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
butterworthBandpass p :: (Double, Double)
p@(Double
wp, Double
_dp) s :: (Double, Double)
s@(Double
ws, Double
_ds) =
   let (Int
n, (Double, Double)
_s') = (Double, Double) -> (Double, Double) -> (Int, (Double, Double))
butterworthParams (Double, Double)
p (Double, Double)
s
   in  (([Double], [Double]) -> ([Double], [Double]))
-> Int -> (Array Int Double, Array Int Double)
butterworth (Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2bp Double
wp Double
ws) Int
n

butterworth ::
   (([Double], [Double]) -> ([Double], [Double])) ->
   Int -> (Array Int Double, Array Int Double)
butterworth :: (([Double], [Double]) -> ([Double], [Double]))
-> Int -> (Array Int Double, Array Int Double)
butterworth ([Double], [Double]) -> ([Double], [Double])
analogToAnalog Int
n =
   forall a b. ([a], [b]) -> (Array Int a, Array Int b)
poly2iir forall a b. (a -> b) -> a -> b
$ Double -> ([Double], [Double]) -> ([Double], [Double])
bilinear Double
1 forall a b. (a -> b) -> a -> b
$ ([Double], [Double]) -> ([Double], [Double])
analogToAnalog forall a b. (a -> b) -> a -> b
$ Int -> ([Double], [Double])
Analog.butterworth Int
n

butterworthParams ::
   (Double, Double) ->
   (Double, Double) ->
   (Int, (Double, Double))
butterworthParams :: (Double, Double) -> (Double, Double) -> (Int, (Double, Double))
butterworthParams (Double
wp, Double
dp) (Double
ws, Double
ds) =
   let n :: Int
n = forall a b. (RealFrac a, Integral b) => a -> b
ceiling forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
log (((Double
1forall a. Fractional a => a -> a -> a
/Double
ds)forall a. Num a => a -> Int -> a
^!Int
2forall a. Num a => a -> a -> a
-Double
1) forall a. Fractional a => a -> a -> a
/ ((Double
1forall a. Fractional a => a -> a -> a
/(Double
1forall a. Num a => a -> a -> a
-Double
dp))forall a. Num a => a -> Int -> a
^!Int
2forall a. Num a => a -> a -> a
-Double
1)) forall a. Fractional a => a -> a -> a
/ Double
2 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
log (Double
ws' forall a. Fractional a => a -> a -> a
/ Double
wp')
       wp' :: Double
wp' = Double -> Double -> Double
prewarp Double
wp Double
1
       ws' :: Double
ws' = Double -> Double -> Double
prewarp Double
ws Double
1
   in  (Int
n, (Double
ws', Double
ds))

wc :: Floating a => Int -> (a, a) -> a
wc :: forall a. Floating a => Int -> (a, a) -> a
wc Int
n (a
ws', a
ds) = a
ws' forall a. Fractional a => a -> a -> a
/ ((a
1forall a. Fractional a => a -> a -> a
/a
ds)forall a. Num a => a -> Int -> a
^!Int
2 forall a. Num a => a -> a -> a
- a
1) forall a. Floating a => a -> a -> a
** (a
1forall a. Fractional a => a -> a -> a
/a
2forall a. Fractional a => a -> a -> a
/forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)


{-# DEPRECATED mkButterworth "Use butterworthLowpass instead" #-}
mkButterworth :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
mkButterworth = (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
butterworthLowpass


-- | Generates lowpass Chebyshev IIR filters

chebyshev1Lowpass, mkChebyshev1 ::
      (Double, Double) -- ^ (wp,dp)
   -> (Double, Double) -- ^ (ws,ds)
   -> (Array Int Double, Array Int Double) -- ^ (b,a)

chebyshev1Lowpass :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
chebyshev1Lowpass (Double
wp,Double
dp) (Double
ws,Double
ds) = forall a b. ([a], [b]) -> (Array Int a, Array Int b)
poly2iir    forall a b. (a -> b) -> a -> b
$
			       Double -> ([Double], [Double]) -> ([Double], [Double])
bilinear Double
1  forall a b. (a -> b) -> a -> b
$
			       Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2lp Double
wp' forall a b. (a -> b) -> a -> b
$
			       Double -> Int -> ([Double], [Double])
Analog.chebyshev1 Double
eps Int
n
    where wp' :: Double
wp' = Double -> Double -> Double
prewarp Double
wp Double
1
          ws' :: Double
ws' = Double -> Double -> Double
prewarp Double
ws Double
1
	  eps :: Double
eps = forall a. Floating a => a -> a
sqrt ((Double
2 forall a. Num a => a -> a -> a
- Double
dp)forall a. Num a => a -> a -> a
*Double
dp) forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
- Double
dp)
	  a :: Double
a   = Double
1 forall a. Fractional a => a -> a -> a
/ Double
ds
	  k1 :: Double
k1  = Double
eps forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
aforall a. Num a => a -> Int -> a
^!Int
2 forall a. Num a => a -> a -> a
- Double
1)
	  k :: Double
k   = Double
wp' forall a. Fractional a => a -> a -> a
/ Double
ws'
	  n :: Int
n   = forall a b. (RealFrac a, Integral b) => a -> b
ceiling forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
acosh (Double
1forall a. Fractional a => a -> a -> a
/Double
k1) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
log ((Double
1 forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
sqrt (Double
1 forall a. Num a => a -> a -> a
- Double
kforall a. Num a => a -> Int -> a
^!Int
2)) forall a. Fractional a => a -> a -> a
/ Double
k)

{-# DEPRECATED mkChebyshev1 "Use chebyshev1Lowpass instead" #-}
mkChebyshev1 :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
mkChebyshev1 = (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
chebyshev1Lowpass


-- | Generates lowpass Inverse Chebyshev IIR filters

chebyshev2Lowpass, mkChebyshev2 ::
      (Double, Double) -- ^ (wp,dp)
   -> (Double, Double) -- ^ (ws,ds)
   -> (Array Int Double, Array Int Double) -- ^ (b,a)

chebyshev2Lowpass :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
chebyshev2Lowpass (Double
wp,Double
dp) (Double
ws,Double
ds) = forall a b. ([a], [b]) -> (Array Int a, Array Int b)
poly2iir    forall a b. (a -> b) -> a -> b
$
			       Double -> ([Double], [Double]) -> ([Double], [Double])
bilinear Double
1  forall a b. (a -> b) -> a -> b
$
			       Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2lp Double
ws' forall a b. (a -> b) -> a -> b
$
			       Double -> Int -> ([Double], [Double])
Analog.chebyshev2 Double
eps Int
n
    where wp' :: Double
wp' = Double -> Double -> Double
prewarp Double
wp Double
1
          ws' :: Double
ws' = Double -> Double -> Double
prewarp Double
ws Double
1
	  eps :: Double
eps = Double
ds forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
1 forall a. Num a => a -> a -> a
- Double
dsforall a. Num a => a -> Int -> a
^!Int
2)
	  g :: Double
g = Double
1 forall a. Num a => a -> a -> a
- Double
dp
	  n :: Int
n   = forall a b. (RealFrac a, Integral b) => a -> b
ceiling forall a b. (a -> b) -> a -> b
$ forall a. Floating a => a -> a
acosh (Double
g forall a. Fractional a => a -> a -> a
/ Double
eps forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
1 forall a. Num a => a -> a -> a
- Double
gforall a. Num a => a -> Int -> a
^!Int
2)) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
acosh (Double
ws' forall a. Fractional a => a -> a -> a
/ Double
wp')

{-# DEPRECATED mkChebyshev2 "Use chebyshev2Lowpass instead" #-}
mkChebyshev2 :: (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
mkChebyshev2 = (Double, Double)
-> (Double, Double) -> (Array Int Double, Array Int Double)
chebyshev2Lowpass