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

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 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)
```