```{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-}
{- |
Copyright   :  (c) Henning Thielemann 2008

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes

Chebyshev lowpass and highpass
-}
module Synthesizer.Plain.Filter.Recursive.Chebyshev where

import Synthesizer.Plain.Filter.Recursive (Passband(Lowpass,Highpass))
import qualified Synthesizer.Plain.Filter.Recursive.SecondOrder as Filt2
-- import qualified Synthesizer.Plain.Filter.Recursive.FirstOrder  as Filt1
import qualified Synthesizer.Plain.Signal   as Sig
-- import qualified Synthesizer.Plain.Modifier as Modifier

import qualified Algebra.VectorSpace           as VectorSpace
import qualified Algebra.Module                as Module
import qualified Algebra.Transcendental        as Trans
import qualified Algebra.Field                 as Field
import qualified Algebra.Ring                  as Ring

import Algebra.Module((*>))

import Number.Complex ((+:),real,imag)
import qualified Number.Complex as Complex

import qualified Prelude as P
import PreludeBase
import NumericPrelude

-- compute the partial filter of the second order from the pole information
{-
It's worth to think it over whether this routine could be used for the Butterworth filter.
But whereas this function is specialized to the zeros of the denominator polynomial
for the Butterworth filter the quadratic factors of the polynomial can be determined
more efficiently than the zeros.
-}
parameterA, parameterB :: (Trans.C a) =>
a -> Complex.T a -> a -> Filt2.Parameter a
parameterA vol z freq =
let re      = real z
im      = imag z
phi     = pi*freq
sinphi  = sin phi
cosphi  = cos phi
cpims   = cosphi + im*sinphi
cmims   = cosphi - im*sinphi
resin2  = (re*sinphi)^2
denom   = - cmims^2 - resin2
c0      = vol * sinphi^2 / denom
in  Filt2.Parameter
c0 (2*c0) c0
(-2*(cpims*cmims - resin2)/denom) ((cpims^2 + resin2)/denom)

runA :: (Trans.C a, Module.C a v) =>
Passband -> Int -> a -> Sig.T a -> Sig.T v -> Sig.T v
runA kind order ratio freqs =
let orderFloat = fromIntegral (2*order)
bn = asinh (ratio/sqrt(1-ratio^2)) / orderFloat
pikn  = take order (map (pi/(2*orderFloat)*) (iterate (2+) 1))
zeros = map (\angle -> cos angle * cosh bn +: (- sin angle * sinh bn)) pikn
{- if ratio == (sqrt 2) then the product of the normalization factors is
2^(1-2*order) -}
makePartialFilter z =
Filt2.run
parameterA (sqrt
((1-real z^2-imag z^2)^2 + 4*imag z^2)) z)
freqs)
in  foldl (.) (ratio *>) (map makePartialFilter zeros)

parameterB a0 z freq =
let re      = real z
im      = imag z
phi     = pi*freq
sinphi  = sin phi
cosphi  = cos phi
spimc   = sinphi + im*cosphi
smimc   = sinphi - im*cosphi
recos2  = (re*cosphi)^2
denom   = smimc^2 + recos2
c0      = (sinphi^2 + a0^2*cosphi^2) / denom
c1      = (sinphi^2 - a0^2*cosphi^2) / denom
in  Filt2.Parameter
c0 (2*c1) c0
(-2*(spimc*smimc - recos2)/denom) (-(spimc^2 + recos2)/denom)

runB :: (Trans.C a, Module.C a v) =>
Passband -> Int -> a -> Sig.T a -> Sig.T v -> Sig.T v
runB kind order ratio freqs =
let orderFloat = fromIntegral (2*order)
bn = (asinh (sqrt(1-ratio^2)/ratio)) / orderFloat
pikn  = take order (map (pi/(2*orderFloat)*) (iterate (2+) 1))
zeros = map (\angle -> (cos angle * cosh bn +: (- sin angle * sinh bn))) pikn
a0s   = map cos pikn
makePartialFilter a0 z =
Filt2.run
parameterB a0 z) freqs)
in  foldl (.) id (zipWith makePartialFilter a0s zeros)

lowpassA, highpassA, lowpassB, highpassB ::
(Trans.C a, VectorSpace.C a v) =>
Int -> a -> Sig.T a -> Sig.T v -> Sig.T v
lowpassA = runA Lowpass
highpassA order ratio freqs =
runA Highpass order ratio (map (0.5-) freqs)

lowpassB = runB Lowpass
highpassB order ratio freqs =
runB Highpass order ratio (map (0.5-) freqs)
```