module Synthesizer.Plain.Filter.Recursive.Chebyshev 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 qualified Algebra.Additive as Additive
import Algebra.Module((*>))
import Number.Complex (real, imag, cis, )
import qualified Number.Complex as Complex
import qualified Data.StorableVector as SV
import Foreign.Storable (Storable)
import qualified Prelude as P
import NumericPrelude.Base
import NumericPrelude.Numeric
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)))
partialParameterA, partialParameterB :: (Trans.C a) =>
Int -> a -> a -> Complex.T a -> Filt2.Parameter a
partialParameterA order ratio freq =
let
bn = (log(1+ratio) log(1ratio^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)
partialParameterB order ratio freq =
let
bn = (log(1+sqrt(1ratio^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)
type ParameterA a = (a, Cascade.Parameter a)
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 ->
Filt2.adjustPassband kind
(flip (partialParameterA order ratio) c) freq) $
circleVec)
type ParameterB a = Cascade.Parameter a
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 ->
Filt2.adjustPassband kind
(flip (partialParameterB order ratio) c) freq) $
circleVec
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
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 -> Filt2.adjustPassband kind $
\freq -> partialParameterA order ratio freq c)
ratios freqs)
in foldl (.) (zipWith (*>) ratios)
(map makePartialFilter (makeCirclePoints order))
runBPole kind order ratios freqs =
let makePartialFilter c =
Filt2.run
(zipWith
(\ratio -> Filt2.adjustPassband kind $
\freq -> partialParameterB order ratio freq c)
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
makePartialFilter c =
Causal.first (Causal.map (\(Pole ratio freq) ->
Filt2.adjustPassband kind
(flip (partialParameterA order ratio) c) freq)) >>>
Filt2.causal
in (\(p, y) -> (p, poleResonance p *> y)) ^>>
(Causal.chainControlled $
map makePartialFilter $
makeCirclePoints order)
causalBPole kind order =
let
makePartialFilter c =
Causal.first (Causal.map (\(Pole ratio freq) ->
Filt2.adjustPassband kind
(flip (partialParameterB 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