{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Synthesizer.Plain.Filter.Recursive.Chebyshev (
ParameterA, parameterA, partialParameterA,
ParameterB, parameterB, partialParameterB,
canonicalizeParameterA,
causalA, runAPole, causalAPole,
causalB, runBPole, causalBPole,
lowpassACausalPole, highpassACausalPole,
lowpassBCausalPole, highpassBCausalPole,
lowpassAPole, highpassAPole,
lowpassBPole, highpassBPole,
makeCirclePoints,
) where
import Synthesizer.Plain.Filter.Recursive (Passband(Lowpass,Highpass), Pole(Pole, poleResonance))
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.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 Number.Complex (real, imag, cis, )
import qualified Number.Complex as Complex
import qualified Data.StorableVector as SV
import Foreign.Storable (Storable)
import NumericPrelude.Numeric
import NumericPrelude.Base
circleList, circleListSlow, _circleListFast :: (Trans.C a) => a -> [Complex.T a]
circleList = circleListSlow
circleListSlow x =
map cis $ map (x*) $ iterate (2+) 1
_circleListFast x =
let z1 = cis x
z2 = z1^2
in iterate (z2*) z1
makeCirclePoints :: (Trans.C a) => Int -> [Complex.T a]
makeCirclePoints order =
take order (circleList (pi / (4 * fromIntegral order)))
partialLowpassParameterA, partialLowpassParameterB :: (Trans.C a) =>
Int -> a -> a -> Complex.T a -> Filt2.Parameter a
partialLowpassParameterA order ratio freq =
let
bn = (log(1+ratio) - log(1-ratio^2)/2) / fromIntegral (2*order)
coshbn = cosh bn
sinhbn = sinh bn
coshbn2 = coshbn^2
sinhbn2 = sinhbn^2
phi = pi*freq
sinphi = sin phi
cosphi = cos phi
sinphi2 = sinphi^2
sinhbnsinphi = sinhbn*sinphi
coshbn2sinphi2 = coshbn2*sinphi2
in \c ->
let re = real c
im = imag c
imss = im * sinhbnsinphi
cpims = cosphi - imss
cmims = cosphi + imss
resin2 = re^2 * coshbn2sinphi2
denom = - cmims^2 - resin2
c0 = (im^2 + sinhbn2)*sinphi2 / denom
in Filt2.Parameter
c0 (2*c0) c0
(-2*(cpims*cmims - resin2)/denom) ((cpims^2 + resin2)/denom)
partialLowpassParameterB order ratio freq =
let
bn = (log(1+sqrt(1-ratio^2)) - log ratio) / fromIntegral (2*order)
coshbn = cosh bn
sinhbn = sinh bn
coshbn2 = coshbn^2
phi = pi*freq
sinphi = sin phi
cosphi = cos phi
sinphi2 = sinphi^2
cosphi2 = cosphi^2
sinhbncosphi = sinhbn*cosphi
in \c ->
let a02cosphi2 = real c ^ 2 * cosphi2
imsc = imag c * sinhbncosphi
spimc = sinphi - imsc
smimc = sinphi + imsc
recos2 = a02cosphi2 * coshbn2
denom = smimc^2 + recos2
c0 = (sinphi2 + a02cosphi2) / denom
c1 = (sinphi2 - a02cosphi2) / denom
in Filt2.Parameter
c0 (2*c1) c0
(-2*(spimc*smimc - recos2)/denom) (-(spimc^2 + recos2)/denom)
{-# INLINE partialParameter #-}
partialParameter ::
(Field.C a) =>
(a -> a -> Complex.T a -> Filt2.Parameter a) ->
Passband -> a -> Complex.T a -> a -> Filt2.Parameter a
partialParameter lowpassParameter kind ratio c freq =
Filt2.adjustPassband kind
(flip (lowpassParameter ratio) c)
freq
{-# INLINE partialParameterA #-}
{-# INLINE partialParameterB #-}
partialParameterA, partialParameterB ::
(Trans.C a) =>
Passband -> Int -> a -> Complex.T a -> a -> Filt2.Parameter a
partialParameterA kind order =
partialParameter (partialLowpassParameterA order) kind
partialParameterB kind order =
partialParameter (partialLowpassParameterB order) kind
type ParameterA a = (a, Cascade.Parameter a)
{-# INLINE parameterA #-}
parameterA ::
(Trans.C a, Storable a) =>
Passband -> Int -> Pole a -> ParameterA a
parameterA kind order =
let circleVec = SV.pack (makeCirclePoints order)
in \ (Pole ratio freq) ->
(ratio,
Cascade.Parameter $
SV.map (\c -> partialParameterA kind order ratio c freq) $
circleVec)
{-# INLINE canonicalizeParameterA #-}
canonicalizeParameterA ::
(Ring.C a, Storable a) =>
ParameterA a -> Cascade.Parameter a
canonicalizeParameterA (amp, Cascade.Parameter p) =
Cascade.Parameter
(SV.switchL SV.empty (\h -> SV.cons (Filt2.amplify amp h)) p)
type ParameterB a = Cascade.Parameter a
{-# INLINE parameterB #-}
parameterB ::
(Trans.C a, Storable a) =>
Passband -> Int -> Pole a -> ParameterB a
parameterB kind order =
let circleVec = SV.pack (makeCirclePoints order)
in \ (Pole ratio freq) ->
Cascade.Parameter $
SV.map (\c -> partialParameterB kind order ratio c freq) $
circleVec
{-# INLINE causalA #-}
causalA :: (Ring.C a, Module.C a v, Storable a, Storable v) =>
Int ->
Causal.T (ParameterA a, v) v
causalA order =
Causal.map (snd.fst) &&& Causal.map (\((ratio,_), y) -> ratio *> y)
>>> Cascade.causal order
{-# INLINE causalB #-}
causalB :: (Ring.C a, Module.C a v, Storable a, Storable v) =>
Int ->
Causal.T (ParameterB a, v) v
causalB =
Cascade.causal
runAPole, runBPole :: (Trans.C a, Module.C a v) =>
Passband -> Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v
runAPole kind order ratios freqs =
let makePartialFilter c =
Filt2.run
(zipWith
(\ratio freq -> partialParameterA kind order ratio c freq)
ratios freqs)
in foldl (.) (zipWith (*>) ratios)
(map makePartialFilter (makeCirclePoints order))
runBPole kind order ratios freqs =
let makePartialFilter c =
Filt2.run
(zipWith
(\ratio freq -> partialParameterB kind order ratio c freq)
ratios freqs)
in foldl (.) id (map makePartialFilter (makeCirclePoints order))
causalAPole, causalBPole :: (Trans.C a, Module.C a v) =>
Passband -> Int -> Causal.T (Pole a, v) v
causalAPole kind order =
let {-# INLINE makePartialFilter #-}
makePartialFilter c =
Causal.first (Causal.map (\(Pole ratio freq) ->
partialParameterA kind order ratio c freq)) >>>
Filt2.causal
in (\(p, y) -> (p, poleResonance p *> y)) ^>>
(Causal.chainControlled $
map makePartialFilter $
makeCirclePoints order)
causalBPole kind order =
let {-# INLINE makePartialFilter #-}
makePartialFilter c =
Causal.first (Causal.map (\(Pole ratio freq) ->
partialParameterB kind order ratio c freq)) >>>
Filt2.causal
in Causal.chainControlled $
map makePartialFilter $
makeCirclePoints order
lowpassACausalPole, highpassACausalPole,
lowpassBCausalPole, highpassBCausalPole ::
(Trans.C a, Module.C a v) =>
Int -> Causal.T (Pole a, v) v
lowpassACausalPole = causalAPole Lowpass
highpassACausalPole = causalAPole Highpass
lowpassBCausalPole = causalBPole Lowpass
highpassBCausalPole = causalBPole Highpass
lowpassAPole, highpassAPole, lowpassBPole, highpassBPole ::
(Trans.C a, Module.C a v) =>
Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v
lowpassAPole = runAPole Lowpass
highpassAPole = runAPole Highpass
lowpassBPole = runBPole Lowpass
highpassBPole = runBPole Highpass