module Synthesizer.Plain.Filter.Recursive.Butterworth where
import Synthesizer.Plain.Filter.Recursive (Passband(Lowpass,Highpass), Pole(Pole))
import qualified Synthesizer.Plain.Filter.Recursive.SecondOrderCascade as Cascade
import qualified Synthesizer.Plain.Filter.Recursive.SecondOrder as Filt2
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 Algebra.Additive as Additive
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
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^21) ** (1 / fromIntegral (2*order))
partialParameterInstable, partialParameter :: (Trans.C a) =>
a -> a -> a -> Filt2.Parameter a
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*(1wc2)/denom) ((wc2+sinw21)/denom)
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
type Parameter a = Cascade.Parameter a
parameter ::
(Trans.C a, Storable a) =>
Passband -> Int -> Pole a -> Parameter a
parameter kind order =
let sinesVec = SV.pack (makeSines order)
in \ (Pole ratio freq) ->
Cascade.Parameter $
SV.map (\sinw ->
Filt2.adjustPassband kind
(flip (partialParameter (partialRatio order ratio)) sinw) freq) $
sinesVec
modifier ::
(Ring.C a, Module.C a v, Storable a, Storable v) =>
Int ->
Modifier.Simple (Cascade.State v) (Parameter a) v v
modifier =
Cascade.modifier
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)
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"
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
makePartialFilter s =
Causal.first (Causal.map (\(Pole ratio freq) ->
Filt2.adjustPassband kind
(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