-----------------------------------------------------------------------------
-- |
-- Module      :  DSP.Filter.Analog.Transform
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Analog prototype filter transforms
---
-- Reference: R&G, pg 258; P&M, pg 698
--
-----------------------------------------------------------------------------

module DSP.Filter.Analog.Transform (
   a_lp2lp, a_lp2hp, a_lp2bp, a_lp2bs,
   substitute, propSubstituteRecip, propSubstituteAlt,
  ) where

import Polynomial.Basic

-- Normalizes a filter

normalize :: ([Double],[Double]) -> ([Double],[Double])
normalize :: ([Double], [Double]) -> ([Double], [Double])
normalize ([Double]
num,[Double]
den) = ([Double]
num',[Double]
den')
    where a0 :: Double
a0 = forall a. [a] -> a
last [Double]
den
          num' :: [Double]
num' = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/ Double
a0) [Double]
num
          den' :: [Double]
den' = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Fractional a => a -> a -> a
/ Double
a0) [Double]
den

-- | Lowpass to lowpass: @s --> s\/wc@

a_lp2lp :: Double -- ^ wc
        -> ([Double],[Double]) -- ^ (b,a)
        -> ([Double],[Double]) -- ^ (b',a')

a_lp2lp :: Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2lp Double
wu ([Double]
num,[Double]
den) = ([Double], [Double]) -> ([Double], [Double])
normalize ([Double]
num',[Double]
den')
    where num' :: [Double]
num' = forall a. Num a => [a] -> [a] -> [a]
polysubst [ Double
0, Double
1forall a. Fractional a => a -> a -> a
/Double
wu ] [Double]
num
          den' :: [Double]
den' = forall a. Num a => [a] -> [a] -> [a]
polysubst [ Double
0, Double
1forall a. Fractional a => a -> a -> a
/Double
wu ] [Double]
den

-- | Lowpass to highpass: @s --> wc\/s@

a_lp2hp :: Double -- ^ wc
        -> ([Double],[Double]) -- ^ (b,a)
        -> ([Double],[Double]) -- ^ (b',a')

a_lp2hp :: Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2hp Double
wu ([Double]
num,[Double]
den) = ([Double], [Double]) -> ([Double], [Double])
normalize ([Double]
num',[Double]
den')
    where nn :: Int
nn   = forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
num
          nd :: Int
nd   = forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
den
          n :: Int
n    = forall a. Ord a => a -> a -> a
max Int
nn Int
nd
          num' :: [Double]
num' = forall a. Num a => [a] -> [a] -> [a]
polysubst [ Double
0, Double
1forall a. Fractional a => a -> a -> a
/Double
wu ] forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ [Double]
num forall a. [a] -> [a] -> [a]
++ forall a. Int -> a -> [a]
replicate (Int
nforall a. Num a => a -> a -> a
-Int
nn) Double
0
          den' :: [Double]
den' = forall a. Num a => [a] -> [a] -> [a]
polysubst [ Double
0, Double
1forall a. Fractional a => a -> a -> a
/Double
wu ] forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ [Double]
den forall a. [a] -> [a] -> [a]
++ forall a. Int -> a -> [a]
replicate (Int
nforall a. Num a => a -> a -> a
-Int
nd) Double
0


-- | Lowpass to bandpass: @s --> (s^2 + wl*wu) \/ (s(wu-wl))@

a_lp2bp :: Double -- ^ wl
        -> Double -- ^ wu
        -> ([Double],[Double]) -- ^ (b,a)
        -> ([Double],[Double]) -- ^ (b',a')

a_lp2bp :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2bp Double
wl Double
wu = ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([ Double
wlforall a. Num a => a -> a -> a
*Double
wu, Double
0, Double
1 ], [ Double
0, Double
wuforall a. Num a => a -> a -> a
-Double
wl ])


-- | Lowpass to bandstop: @s --> (s(wu-wl)) \/ (s^2 + wl*wu)@

a_lp2bs :: Double -- ^ wl
        -> Double -- ^ wu
        -> ([Double],[Double]) -- ^ (b,a)
        -> ([Double],[Double]) -- ^ (b',a')

a_lp2bs :: Double -> Double -> ([Double], [Double]) -> ([Double], [Double])
a_lp2bs Double
wl Double
wu = ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([ Double
0, Double
wuforall a. Num a => a -> a -> a
-Double
wl ], [ Double
wuforall a. Num a => a -> a -> a
*Double
wl, Double
0, Double
1 ])



substitute ::
   ([Double],[Double]) -> ([Double],[Double]) -> ([Double],[Double])
substitute :: ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den) = ([Double], [Double]) -> ([Double], [Double])
normalize ([Double]
num',[Double]
den')
    where num' :: [Double]
num' = forall a. Num a => [a] -> [[a]] -> [a]
polyPolySubst [Double]
nsub forall a b. (a -> b) -> a -> b
$ [Double] -> [[Double]]
weightedPowers forall a b. (a -> b) -> a -> b
$ [Double]
num
          den' :: [Double]
den' = forall a. Num a => [a] -> [[a]] -> [a]
polyPolySubst [Double]
nsub forall a b. (a -> b) -> a -> b
$ [Double] -> [[Double]]
weightedPowers forall a b. (a -> b) -> a -> b
$ [Double]
den
          weightedPowers :: [Double] -> [[Double]]
weightedPowers = forall a b c. (a -> b -> c) -> b -> a -> c
flip (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> [a] -> [a]
polyscale) [[Double]]
dsubPowers
          dsubPowers :: [[Double]]
dsubPowers = forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
take Int
m forall a b. (a -> b) -> a -> b
$ forall a. (a -> a) -> a -> [a]
iterate (forall a. Num a => [a] -> [a] -> [a]
polymult [Double]
dsub) [Double
1]
          m :: Int
m = forall a. Ord a => a -> a -> a
max (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
num) (forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
den)

substituteAlt ::
   ([Double],[Double]) -> ([Double],[Double]) -> ([Double],[Double])
substituteAlt :: ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substituteAlt ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den) = ([Double], [Double]) -> ([Double], [Double])
normalize ([Double]
num',[Double]
den')
    where m :: Int
m    = 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 {t :: * -> *} {a}. (Foldable t, Num a) => t [a] -> [a]
step3 forall a b. (a -> b) -> a -> b
$ forall {t}. Integral t => t -> [[Double]] -> [[Double]]
step2 (Int
0::Int) forall a b. (a -> b) -> a -> b
$ forall {t}. Integral t => t -> [Double] -> [[Double]]
step1 Int
m forall a b. (a -> b) -> a -> b
$ [Double]
num
          den' :: [Double]
den' = forall {t :: * -> *} {a}. (Foldable t, Num a) => t [a] -> [a]
step3 forall a b. (a -> b) -> a -> b
$ forall {t}. Integral t => t -> [[Double]] -> [[Double]]
step2 (Int
0::Int) forall a b. (a -> b) -> a -> b
$ forall {t}. Integral t => t -> [Double] -> [[Double]]
step1 Int
m forall a b. (a -> b) -> a -> b
$ [Double]
den
          step1 :: t -> [Double] -> [[Double]]
step1 t
_ []     = []
          step1 t
n (Double
x:[Double]
xs) = forall a b. (a -> b) -> [a] -> [b]
map (Double
xforall a. Num a => a -> a -> a
*) (forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [Double]
dsub t
n) forall a. a -> [a] -> [a]
: t -> [Double] -> [[Double]]
step1 (t
nforall a. Num a => a -> a -> a
-t
1) [Double]
xs
          step2 :: t -> [[Double]] -> [[Double]]
step2 t
_ []     = []
          step2 t
n ([Double]
x:[[Double]]
xs) = forall a. Num a => [a] -> [a] -> [a]
polymult (forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [Double]
nsub t
n) [Double]
x forall a. a -> [a] -> [a]
: t -> [[Double]] -> [[Double]]
step2 (t
nforall a. Num a => a -> a -> a
+t
1) [[Double]]
xs
          step3 :: t [a] -> [a]
step3 t [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] t [a]
x


propSubstituteRecip :: ([Double],[Double]) -> ([Double],[Double]) -> Bool
propSubstituteRecip :: ([Double], [Double]) -> ([Double], [Double]) -> Bool
propSubstituteRecip ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den) =
   let ([Double]
x,[Double]
y) =  ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
nsub,[Double]
dsub) ([Double]
num,[Double]
den)
   in  ([Double]
y,[Double]
x) forall a. Eq a => a -> a -> Bool
== ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double]
dsub,[Double]
nsub) ([Double]
den,[Double]
num)


propSubstituteAlt :: ([Double],[Double]) -> ([Double],[Double]) -> Bool
propSubstituteAlt :: ([Double], [Double]) -> ([Double], [Double]) -> Bool
propSubstituteAlt ([Double], [Double])
q ([Double], [Double])
p   =   ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substitute ([Double], [Double])
q ([Double], [Double])
p forall a. Eq a => a -> a -> Bool
== ([Double], [Double])
-> ([Double], [Double]) -> ([Double], [Double])
substituteAlt ([Double], [Double])
q ([Double], [Double])
p