-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.Analog.Prototype
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Module for generating analog filter prototypes
--
-----------------------------------------------------------------------------

-- Notes (mainly for self):

-- The gain of an analog filter is

--    gain = abs $ realPart $ product zeros / product poles
--         = abs $ b_m / a_n

-- For a Butterworth filter, the product of the poles is one, so we don't
-- have to worry about any gain.

-- For a Chebyshev 1 filter, the product of the poles is a_n, which is
-- the head of the polynomial.  We make this b_0 to set the gain in the
-- passband.

-- For a Chebyshev 2 filter, we use the full gain formula because we want
-- to set the gain to unity at DC.

-- TODO: Do we want to include Bessel filters?

module DSP.Filter.Analog.Prototype where

import Data.Complex (Complex((:+)), realPart)

import Polynomial.Basic (roots2poly)

-- | Generates Butterworth filter prototype

butterworth :: Int -- ^ N
	    -> ([Double],[Double]) -- ^ (b,a)

butterworth :: Int -> ([Double], [Double])
butterworth Int
n = ([Double]
num, [Double]
den)
    where poles :: [Complex Double]
poles = [ (-forall {a} {a}. (Floating a, Integral a) => a -> a
u Int
k) forall a. a -> a -> Complex a
:+ (forall {a} {a}. (Floating a, Integral a) => a -> a
w Int
k) | Int
k <- [Int
0..(Int
nforall a. Num a => a -> a -> a
-Int
1)] ]
	  u :: a -> a
u a
k = forall a. Floating a => a -> a
sin (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  w :: a -> a
w a
k = forall a. Floating a => a -> a
cos (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  num :: [Double]
num = [ Double
1 ]
	  den :: [Double]
den = forall a b. (a -> b) -> [a] -> [b]
map forall a. Complex a -> a
realPart forall a b. (a -> b) -> a -> b
$ forall a. Num a => [a] -> [a]
roots2poly forall a b. (a -> b) -> a -> b
$ [Complex Double]
poles

-- | Generates Chebyshev filter prototype

chebyshev1 :: Double -- ^ epsilon
	   -> Int -- ^ N
	   -> ([Double],[Double]) -- ^ (b,a)

chebyshev1 :: Double -> Int -> ([Double], [Double])
chebyshev1 Double
eps Int
n = ([Double]
num, [Double]
den)
    where poles :: [Complex Double]
poles = [ (-forall {a}. Integral a => a -> Double
u Int
k) forall a. a -> a -> Complex a
:+ (forall {a}. Integral a => a -> Double
w Int
k) | Int
k <- [Int
0..(Int
nforall a. Num a => a -> a -> a
-Int
1)] ]
	  u :: a -> Double
u a
k = forall a. Floating a => a -> a
sinh Double
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  w :: a -> Double
w a
k = forall a. Floating a => a -> a
cosh Double
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  num :: [Double]
num = [ Double
gain ]
	  den :: [Double]
den = forall a b. (a -> b) -> [a] -> [b]
map forall a. Complex a -> a
realPart forall a b. (a -> b) -> a -> b
$ forall a. Num a => [a] -> [a]
roots2poly forall a b. (a -> b) -> a -> b
$ [Complex Double]
poles
	  v0 :: Double
v0 = forall a. Floating a => a -> a
asinh (Double
1forall a. Fractional a => a -> a -> a
/Double
eps) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
	  gain :: Double
gain =
             if forall a. Integral a => a -> Bool
even Int
n
               then forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head [Double]
den forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
1 forall a. Num a => a -> a -> a
+ Double
epsforall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int))
	       else forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head [Double]
den

-- | Generates Inverse Chebyshev filter prototype

chebyshev2 :: Double -- ^ epsilon
	   -> Int -- ^ N
	   -> ([Double],[Double]) -- ^ (b,a)

chebyshev2 :: Double -> Int -> ([Double], [Double])
chebyshev2 Double
eps Int
n = ([Double]
num, [Double]
den)
    where zeros :: [Complex Double]
zeros = [ Double
0 forall a. a -> a -> Complex a
:+ Double
1 forall a. Fractional a => a -> a -> a
/ forall {a} {a}. (Floating a, Integral a) => a -> a
wz Int
k | Int
k <- [Int
0..(Int
nforall a. Num a => a -> a -> a
-Int
1)], Int
2forall a. Num a => a -> a -> a
*Int
kforall a. Num a => a -> a -> a
+Int
1 forall a. Eq a => a -> a -> Bool
/= Int
n ]
	  poles :: [Complex Double]
poles = [ Complex Double
1 forall a. Fractional a => a -> a -> a
/ ((-forall {a}. Integral a => a -> Double
u Int
k) forall a. a -> a -> Complex a
:+ (forall {a}. Integral a => a -> Double
w Int
k)) | Int
k <- [Int
0..(Int
nforall a. Num a => a -> a -> a
-Int
1)] ]
	  wz :: a -> a
wz a
k = forall a. Floating a => a -> a
cos (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  u :: a -> Double
u a
k = forall a. Floating a => a -> a
sinh Double
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  w :: a -> Double
w a
k = forall a. Floating a => a -> a
cosh Double
v0 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (forall a b. (Integral a, Num b) => a -> b
fromIntegral (a
2forall a. Num a => a -> a -> a
*a
kforall a. Num a => a -> a -> a
+a
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2forall a. Num a => a -> a -> a
*Int
n))
	  num :: [Double]
num = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
*Double
gain) forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a. Complex a -> a
realPart forall a b. (a -> b) -> a -> b
$ forall a. Num a => [a] -> [a]
roots2poly forall a b. (a -> b) -> a -> b
$ [Complex Double]
zeros
	  den :: [Double]
den =               forall a b. (a -> b) -> [a] -> [b]
map forall a. Complex a -> a
realPart forall a b. (a -> b) -> a -> b
$ forall a. Num a => [a] -> [a]
roots2poly forall a b. (a -> b) -> a -> b
$ [Complex Double]
poles
	  v0 :: Double
v0 = forall a. Floating a => a -> a
asinh (Double
1forall a. Fractional a => a -> a -> a
/Double
eps) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
	  gain :: Double
gain = forall a. Num a => a -> a
abs forall a b. (a -> b) -> a -> b
$ forall a. Complex a -> a
realPart forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Complex Double]
poles forall a. Fractional a => a -> a -> a
/ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Complex Double]
zeros