{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-}moduleFilter.TwoWaywhereimportFilter.BasicimportSynthesizer.Plain.Interpolation(minLength,)importqualifiedSynthesizer.Plain.InterpolationasIpimportqualifiedSynthesizer.Plain.InterpolationasInterpolationimportAlgebra.Module(linearComb,(*>))importqualifiedAlgebra.ModuleasModuleimportqualifiedAlgebra.RealFieldasRealFieldimportqualifiedAlgebra.RingasRingimportqualifiedAlgebra.AdditiveasAdditiveimportNumber.Complex(cis)importSynthesizer.Utility(nest)importqualifiedData.ListasListimportPreludeBasehiding(take)importNumericPrelude{-| A TwoWay.Signal stores values of the past and the future -}dataSignalv=Signal{past,future::[v]}deriving(Show,Eq){-| Take n values starting from time zero.
If you want clips from elsewhere,
call 'take' after 'delay'. -}take::Int->Signalv->[v]taken(Signal_x)=List.takenxzipSignalWith::(a->b->c)->Signala->Signalb->SignalczipSignalWithf(SignalxPastxFuture)(SignalyPastyFuture)=(Signal(zipWithfxPastyPast)(zipWithfxFutureyFuture)){-| Take the value at time zero. -}origin::Ring.Ca=>Signala->aorigin(Signal_(x:_))=xorigin_=0{-| A signal that consists entirely of ones -}ones::Ring.Ca=>Signalaones=Signal(repeat1)(repeat1){-| shift signal in time,
keep all values but if required pad with zeros -}delay::(Additive.Cv)=>Int->Signalv->Signalvdelay=delayGendelayOncedelayPad::v->Int->Signalv->SignalvdelayPadz=delayGen(delayPadOncez){-| shift signal in time,
zero values at either ends will be flushed -}delayOpt::(Eqv,Additive.Cv)=>Int->Signalv->SignalvdelayOpt=delayGendelayOptOnce{-| Delay by one sample. -}delayOnce::(Additive.Cv)=>Signalv->Signalv--delayOnce (Signal [] []) = ([],[])delayOnce(Signal[]ys)=Signal[](zero:ys)delayOnce(Signal(x:xs)ys)=Signalxs(x:ys)delayPadOnce::v->Signalv->Signalv--delayPadOnce _ (Signal [] []) = ([],[])delayPadOncez(Signal[]ys)=Signal[](z:ys)delayPadOnce_(Signal(x:xs)ys)=Signalxs(x:ys)delayOptOnce::(Eqv,Additive.Cv)=>Signalv->Signalv--delayOptOnce (Signal [] []) = Signal [] []delayOptOnce(Signal[]ys)=Signal[](zero:ys)delayOptOnce(Signal(x:xs)[])=Signalxs(ifx==zerothen[]elsex:[])delayOptOnce(Signal(x:xs)ys)=Signalxs(x:ys){-| General routine that supports delaying and prefetching
using a general one-sample delaying routine. -}delayGen::(Signalv->Signalv)->Int->Signalv->Signalv{- Using this optimization applications of recursive filters
with zero initial conditions represented by an empty list will fail.
This is because in this case the value of the first item of the future list
depends on the first item of the input future list,
whereas normally the first future value depends on no input future value,
at all.
delayGen _ _ (Signal [] []) = Signal [] []
cf. the next example -}delayGendelOncet=ift<0thenreverseTwoWay.nest(negatet)delOnce.reverseTwoWayelsenesttdelOncereverseTwoWay::Signalv->SignalvreverseTwoWay(Signalxy)=Signalyxinstance(Additive.Cv)=>Additive.C(Signalv)wherezero=Signalzerozero(+)(Signaly0y1)(Signalx0x1)=Signal(y0+x0)(y1+x1)(-)(Signaly0y1)(Signalx0x1)=Signal(y0-x0)(y1-x1)negate(Signalx0x1)=Signal(negatex0)(negatex1)instance(Module.Cav)=>Module.Ca(Signalv)where(*>)s(Signalx0x1)=Signal(s*>x0)(s*>x1){-| for a Signal this means a reversion of the elements -}flipPair::(a,b)->(b,a)flipPair(x,y)=(y,x){- This example simulates what happens if you call
apply (Feedback (Serial []) (Prim (Delay 1)) []) ([],[1])
depending on the implementation of delayGen it may work or
loop infinitely when yFuture is computed.
It's even before the first element of yFuture is computed.
Note, that the equivalent
apply (Feedback (Serial []) (Prim (Delay 1)) [0]) ([],[1])
works! (i.e. set yPast = [0] -}testDelayGen::SignalDoubletestDelayGen=letyPast=[]x=Signal[][1]y=SignalyPastyFutureSignal_yFuture=delayOnce(x+y)-- Signal _ yFuture = delayOptOnce (add x y)-- Signal _ yFuture = delayGen delayOnce 1 (add x y)inSignalyPast(List.take10yFuture){-| Unmodulated non-recursive filter -}nonRecursiveFilter::Module.Cav=>[a]->Signalv->SignalvnonRecursiveFiltermx=linearCombm(iteratedelayOncex){-| Modulated non-recursive filter.
The number of values before time 0 (past) or
the filter mask lengths must be at most finite. -}nonRecursiveFilterMod::Module.Cav=>Signal[a]->Signalv->SignalvnonRecursiveFilterMod(Signalmpremsuf)x=let(pre,suf)=unzip(map(\(Signalab)->(a,b))(iteratedelayOncex))inSignal(zipWithlinearCombmprepre)(zipWithlinearCombmsufsuf){-| Interpolation allowing negative frequencies,
but requires storage of all past values. -}interpolatePaddedZero::(Orda,RealField.Ca)=>b->Interpolation.Tab->a->Signala->Signalb->SignalbinterpolatePaddedZerozipphasefs(SignalxPastxFuture)=let(phInt,phFrac)=splitFractionphasexPadded=Signal(xPast++repeatz)(xFuture++repeatz)ininterpolateCoreipphFracfs(delayPadz(Ip.offsetip-phInt)xPadded)interpolatePaddedCyclic::(Orda,RealField.Ca)=>Interpolation.Tab->a->Signala->Signalb->SignalbinterpolatePaddedCyclicipphasefs(SignalxPastxFuture)=let(phInt,phFrac)=splitFractionphasexCyclic=xFuture++reversexPastininterpolateCoreipphFracfs-- mod is for efficiency, only(delayPad(error"interpolate: infinite signal needs no zero padding")(mod(Ip.offsetip-phInt)(lengthxCyclic))(Signal(cycle(reversexCyclic))(cyclexCyclic)))-- note that the extrapolation may miss some of the first and some of the last pointsinterpolatePaddedExtrapolation::(Orda,RealField.Ca)=>Interpolation.Tab->a->Signala->Signalb->SignalbinterpolatePaddedExtrapolationipphasefsx=interpolateCoreip(phase-fromIntegral(Ip.offsetip))fsxinterpolateCore::(Orda,Ring.Ca)=>Interpolation.Tab->a->Signala->Signalb->SignalbinterpolateCoreipphase(SignalfreqPastfreqFuture)x=Signal(interpolateHalfWayip(1-phase)freqPast(delayPadOnce(error"interpolateCore: infinite signal needs no zero padding")(reverseTwoWayx)))(interpolateHalfWayipphasefreqFuturex)interpolateHalfWay::(Orda,Ring.Ca)=>Interpolation.Tab->a->[a]->Signalb->[b]interpolateHalfWayipphasefreqs(SignalxPastxFuture)=ifphase>=1&&minLength(1+Ip.numberip)xFuturetheninterpolateHalfWayip(phase-1)freqs(Signal(headxFuture:xPast)(tailxFuture))elseifphase<0&&minLength1xPasttheninterpolateHalfWayip(phase+1)freqs(Signal(tailxPast)(headxPast:xFuture))elseIp.funcipphasexFuture:interpolateHalfWayip(phase+headfreqs)(tailfreqs)(SignalxPastxFuture){-| Description of a basic filter that can be used in larger networks. -}dataTtav=Mask[a]{-^ A static filter described by its mask -}|ModMask(Signal[a]){-^ A modulated filter described by a list of masks -}|FracDelay(Interpolation.Ttv)t{-^ Delay the signal by a fractional amount of samples.
This is achieved by interpolation. -}|ModFracDelay(Interpolation.Ttv)(Signalt){-^ Delay with varying delay times. -}|DelayInt{-^ Delay the signal by given amount of samples. -}|Past[v]{-^ Replace the past by the given one.
This must be present in each recursive filter cycle
to let the magic work! -}instanceFilterSignalTwhereapply(Maskm)=nonRecursiveFiltermapply(ModMaskm)=nonRecursiveFilterModmapply(FracDelayipt)=interpolatePaddedZerozeroip(-t)onesapply(ModFracDelayipts)=interpolatePaddedZerozeroip(-origints)(ts-delay(-1)ts+ones)apply(Delayt)=delaytapply(Pastx)=Signalx.future{- This is in principle the same as for one way filters.
How can one merge them? -}transferFunction(Maskm)w=linearCombm(screw(negatew))transferFunction(FracDelay_t)w=cis(negatew*t)transferFunction(Delayt)w=cis(negatew*fromIntegralt)transferFunction(Past_)_=1transferFunction__=error"transfer function can't be computed for modulated filters"