{-# 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

All recursive filters with real coefficients
can be decomposed into first order and second order filters with real coefficients.
This follows from the Fundamental theorem of algebra.
-}
module Synthesizer.Plain.Filter.Recursive.SecondOrder where

import Synthesizer.Plain.Filter.Recursive (Passband(Lowpass,Highpass))
import qualified Synthesizer.Plain.Signal   as Sig
-- import qualified Synthesizer.Plain.Modifier as Modifier
-- import qualified Synthesizer.Plain.Control as Ctrl

-- 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 Data.List (zipWith6)

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

import qualified Prelude as P
import PreludeBase
import NumericPrelude


{- | Parameters for a general recursive filter of 2nd order. -}
data Parameter a =
   Parameter {c0, c1, c2, d1, d2 :: !a}
       deriving Show

{- | Given the filter parameters of a lowpass filter,
     turn them into highpass parameters, if requested filter type is Highpass -}
{-# INLINE adjustPassband #-}
adjustPassband :: (Ring.C a) =>
   Passband -> Parameter a -> Parameter a
adjustPassband kind p =
   case kind of
      Lowpass  -> p
      Highpass -> Parameter (c0 p) (- c1 p) (c2 p) (- d1 p) (d2 p)

{-# INLINE step #-}
step :: (Ring.C a, Module.C a v) =>
   Parameter a -> v -> State ((v,v),(v,v)) v
step c u0 = State $ \((u1,u2),(y1,y2)) ->
   let y0 =
          c0 c *> u0 +
          c1 c *> u1 + d1 c *> y1 +
          c2 c *> u2 + d2 c *> y2
   in  (y0, ((u0,u1),(y0,y1)))


{-# INLINE runInit #-}
runInit :: (Ring.C a, Module.C a v) =>
   ((v,v),(v,v)) -> Sig.T (Parameter a) -> Sig.T v -> Sig.T v
runInit ((u0init,u1init),(y0init,y1init)) control input =
   let u0s = input
       u1s = u0init:u0s
       u2s = u1init:u1s
       y1s = y0init:y0s
       y2s = y1init:y1s
       y0s = zipWith6
          (\c u0 u1 u2 y1 y2 ->
              c0 c *> u0 +
              c1 c *> u1 + d1 c *> y1 +
              c2 c *> u2 + d2 c *> y2)
          control u0s u1s u2s y1s y2s
   in  y0s

{-# INLINE run #-}
run :: (Ring.C a, Module.C a v) =>
   Sig.T (Parameter a) -> Sig.T v -> Sig.T v
run = runInit zero