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

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 qualified Algebra.Additive              as Additive

import Algebra.Module((*>))

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

-- import Control.Monad.State (State(..), evalState)

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
             (map (Filt2.adjustPassband kind .
                   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
             (map (Filt2.adjustPassband kind .
                   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)