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

Butterworth lowpass and highpass
-}
module Synthesizer.Plain.Filter.Recursive.Butterworth 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.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 Control.Monad.State (State(..), evalState)

import qualified Prelude as P
import PreludeBase
import NumericPrelude



sineListSlow, sineListFast :: (Trans.C a) => a -> [a]
sineListSlow x = map sin (map (x*) (iterate (2+) 1))

sineListFast x =
   let sinx  = sin x
       cos2x = 2 - 4*sinx*sinx
       --cos2x = 2 * cos (2*x)
       sines = (-sinx) : sinx :
                  zipWith (\y1 y0 -> cos2x * y0 - y1) sines (tail sines)
   in  tail sines

parameterInstable, parameter :: (Trans.C a) =>
   a -> a -> a -> Filt2.Parameter a

{- must handle infinite values when 'freq' approaches 0.5 -}
parameterInstable ratio sinw freq =
   let wc    = ratio * tan (pi*freq)
       sinw2 = 2 * wc * sinw
       wc2   = wc * wc
       denom = wc2 + sinw2 + 1
       c0    = wc2 / denom
   in  Filt2.Parameter c0 (2*c0) c0
          (2*(1-wc2)/denom) ((-wc2+sinw2-1)/denom)

-- using ratio disallows simplification by trigonometric Pythagoras' theorem
parameter ratio sinw freq =
   let phi     = pi*freq
       sinsin  = sin (2*phi) * sinw
       sinphi2 = (sin phi)^2
       ratsin  = 1 + (ratio^2-1) * sinphi2
       denom   = ratsin + ratio*sinsin
       c0      = ratio^2 * sinphi2 / denom
   in  Filt2.Parameter c0 (2*c0) c0
          (2*(1-(ratio^2+1)*sinphi2)/denom) ((ratio*sinsin-ratsin)/denom)

run :: (Trans.C a, Module.C a v) =>
   Passband -> Int -> a -> Sig.T a -> Sig.T v -> Sig.T v
run kind order ratio freqs =
   let orderFloat = fromIntegral (2*order)
       sines = take (div order 2) (sineListFast (pi/orderFloat))
       -- the filter amplifies frequency 0 with factor 1
       -- and frequency 'freq' with factor 'ratio'
       wcoef = (1/ratio^2-1)**(-1/orderFloat)
       makePartialFilter s =
          Filt2.run (map (Filt2.adjustPassband kind . parameter wcoef s) freqs)
   in  foldl (.) id (map makePartialFilter sines)

lowpass, highpass :: (Trans.C a, Module.C a v) =>
   Int -> a -> Sig.T a -> Sig.T v -> Sig.T v
lowpass = run Lowpass
highpass order ratio freqs =
   run Highpass order ratio (map (0.5-) freqs)