```{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
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), Pole(Pole))
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 Synthesizer.Causal.Process as Causal
import Control.Arrow ((>>>), )

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 Data.StorableVector as SV
import Foreign.Storable (Storable)

import qualified Prelude as P
import NumericPrelude.Base
import NumericPrelude.Numeric

sineList, sineListSlow, sineListFast :: (Trans.C a) => a -> [a]
sineList = sineListFast

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

makeSines :: (Trans.C a) =>
Int -> [a]
makeSines order =
take (checkedHalf "makeSines" order) (sineList (pi / fromIntegral (2*order)))

partialRatio :: (Trans.C a) =>
Int -> a -> a
partialRatio order ratio =
(1/ratio^2-1) ** (-1 / fromIntegral (2*order))

partialParameterInstable, partialParameter :: (Trans.C a) =>
a -> a -> a -> Filt2.Parameter a

{- must handle infinite values when 'freq' approaches 0.5 -}
partialParameterInstable ratio freq sinw =
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
partialParameter ratio freq =
let phi      = pi*freq
rsin2phi = ratio * sin (2*phi)
cosphi2  = cos phi ^ 2
c0d      = (ratio * sin phi) ^ 2
d1d      = (cosphi2 - c0d) * 2
ratsin   = cosphi2 + c0d
in  \sinw ->
let rsinsin = rsin2phi * sinw
denom   = rsinsin + ratsin
d2d     = rsinsin - ratsin
c0      = c0d / denom
d1      = d1d / denom
d2      = d2d / denom
in  Filt2.Parameter c0 (2*c0) c0 d1 d2

-- * use second order filter parameters for control

type Parameter a = Cascade.Parameter a

{-# INLINE parameter #-}
parameter ::
(Trans.C a, Storable a) =>
Passband -> Int -> Pole a -> Parameter a
parameter kind order =
-- I hope that the 'let' is floated out of a 'map'
let sinesVec = SV.pack (makeSines order)
in  \ (Pole ratio freq) ->
SV.map (\sinw ->
(flip (partialParameter (partialRatio order ratio)) sinw) freq) \$
sinesVec

{-# INLINE modifier #-}
modifier ::
(Ring.C a, Module.C a v, Storable a, Storable v) =>
Int ->
Modifier.Simple (Cascade.State v) (Parameter a) v v
modifier =

{-# INLINE causal #-}
causal :: (Ring.C a, Module.C a v, Storable a, Storable v) =>
Int ->
Causal.T (Parameter a, v) v
causal order =
Cascade.causal (checkedHalf "causal" order)

{-# INLINE checkedHalf #-}
checkedHalf :: String -> Int -> Int
checkedHalf funcName order =
let (order2,r) = divMod order 2
in  if r==0
then order2
else error \$ "Butterworth."++funcName++": order must be even"

{-
lowpassCausal, highpassCausal :: (Trans.C a, Module.C a v) =>
Int -> Causal.T (Parameter a, v) v
lowpassCausal  = causal Lowpass
highpassCausal = causal Highpass

lowpass, highpass :: (Trans.C a, Module.C a v) =>
Int -> Sig.T (Parameter a) -> Sig.T v -> Sig.T v
lowpass  = run Lowpass
highpass = run Highpass
-}

-- * directly use frequency as control parameter

{- |
When called as @runPole kind order ratio freqs@,
the filter amplifies frequency 0 with factor 1
and frequency @freq@ with factor @ratio@.

It uses the frequency and ratio information directly
and thus cannot benefit from efficient parameter interpolation
(asynchronous run of a ControlledProcess.
-}
runPole :: (Trans.C a, Module.C a v) =>
Passband -> Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v
runPole kind order ratios freqs =
let makePartialFilter s =
Filt2.run (zipWith (\ratio ->
Filt2.adjustPassband kind \$ \freq ->
partialParameter (partialRatio order ratio) freq s) ratios freqs)
in  foldl (.) id (map makePartialFilter (makeSines order))

causalPole :: (Trans.C a, Module.C a v) =>
Passband -> Int -> Causal.T (Pole a, v) v
causalPole kind order =
let {-# INLINE makePartialFilter #-}
makePartialFilter s =
Causal.first (Causal.map (\(Pole ratio freq) ->
(flip (partialParameter (partialRatio order ratio)) s) freq)) >>>
Filt2.causal
in  Causal.chainControlled \$ map makePartialFilter \$ makeSines order

lowpassCausalPole, highpassCausalPole :: (Trans.C a, Module.C a v) =>
Int -> Causal.T (Pole a, v) v
lowpassCausalPole  = causalPole Lowpass
highpassCausalPole = causalPole Highpass

lowpassPole, highpassPole :: (Trans.C a, Module.C a v) =>
Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v
lowpassPole  = runPole Lowpass
highpassPole = runPole Highpass
```