{-# OPTIONS -fno-implicit-prelude #-}
{- |
Copyright   :  (c) Henning Thielemann 2008
License     :  GPL

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes
-}
module Synthesizer.Dimensional.RateAmplitude.Filter (
   {- * Non-recursive -}

   {- ** Amplification -}
   amplify,
   amplifyDimension,
   negate,
   envelope,
   envelopeVector,
   envelopeVectorDimension,
   {- ** Filter operators from calculus -}
   differentiate,

   {- ** Smooth -}
   meanStatic,
   mean,

   {- ** Delay -}
   delay,
   phaseModulation,
   frequencyModulation,
   frequencyModulationDecoupled,
   phaser,
   phaserStereo,


   {- * Recursive -}

   {- ** Without resonance -}
   firstOrderLowpass,
   firstOrderHighpass,
   butterworthLowpass,
   butterworthHighpass,
   chebyshevALowpass,
   chebyshevAHighpass,
   chebyshevBLowpass,
   chebyshevBHighpass,
   {- ** With resonance -}
   universal,
   FiltR.highpassFromUniversal,
   FiltR.bandpassFromUniversal,
   FiltR.lowpassFromUniversal,
   moogLowpass,

   {- ** Allpass -}
   allpassCascade,

   {- ** Reverb -}
   comb,
   combProc,

   {- ** Filter operators from calculus -}
   integrate,
) where

import qualified Synthesizer.Dimensional.Rate.Filter as FiltR
import qualified Synthesizer.Dimensional.Amplitude.Filter       as FiltV
-- import qualified Synthesizer.Dimensional.Amplitude.Displacement as MiscV
-- import qualified Synthesizer.Dimensional.Amplitude.Cut          as CutV
import qualified Synthesizer.Dimensional.ControlledProcess as CProc
import qualified Synthesizer.Dimensional.Process as Proc
-- import qualified Synthesizer.Dimensional.Rate as Rate

-- import Synthesizer.Dimensional.Process ((.:), (.^), )

import qualified Synthesizer.Dimensional.Abstraction.Flat as Flat
import qualified Synthesizer.Dimensional.Abstraction.Homogeneous as Hom

import qualified Synthesizer.Dimensional.RatePhantom as RP

import qualified Synthesizer.Dimensional.Straight.Signal      as SigS
import qualified Synthesizer.Dimensional.RateAmplitude.Signal as SigA
import qualified Synthesizer.Dimensional.RateWrapper      as SigP
-- import qualified Synthesizer.Dimensional.Amplitude.Signal as SigPA
import qualified Synthesizer.State.Signal as Sig
import Synthesizer.Plain.Signal (Modifier)

import Synthesizer.Dimensional.RateAmplitude.Signal
   (toTimeScalar, toFrequencyScalar, DimensionGradient, )

import qualified Synthesizer.Generic.SampledValue as Sample

import qualified Synthesizer.Frame.Stereo as Stereo

-- import qualified Synthesizer.State.Displacement as Disp
import qualified Synthesizer.State.Interpolation as Interpolation
import qualified Synthesizer.State.Filter.Delay as Delay
import qualified Synthesizer.Plain.Filter.Recursive.FirstOrder  as Filt1
import qualified Synthesizer.Plain.Filter.Recursive.Allpass     as Allpass
import qualified Synthesizer.Plain.Filter.Recursive.Universal   as UniFilter
import qualified Synthesizer.Plain.Filter.Recursive.Moog        as Moog
import qualified Synthesizer.Plain.Filter.Recursive.Butterworth as Butter
import qualified Synthesizer.Plain.Filter.Recursive.Chebyshev   as Cheby
import qualified Synthesizer.Plain.Filter.Recursive    as FiltR
import qualified Synthesizer.State.Filter.Recursive.Integration as Integrate
import qualified Synthesizer.State.Filter.Recursive.MovingAverage as MA
import qualified Synthesizer.State.Filter.NonRecursive as FiltNR

import qualified Synthesizer.Storable.Signal as SigSt
import qualified Synthesizer.Storable.Filter.Recursive.Comb as Comb

import qualified Number.DimensionTerm        as DN
import qualified Algebra.DimensionTerm       as Dim

import Number.DimensionTerm ((&*&), (&/&))

import qualified Number.NonNegative     as NonNeg

import qualified Algebra.Transcendental as Trans
import qualified Algebra.RealField      as RealField
import qualified Algebra.Field          as Field
import qualified Algebra.Real           as Real
import qualified Algebra.Ring           as Ring
import qualified Algebra.Additive       as Additive
import qualified Algebra.VectorSpace    as VectorSpace
import qualified Algebra.Module         as Module

-- import Control.Monad(liftM2)

import NumericPrelude hiding (negate)
import PreludeBase as P
import Prelude ()


{- | The amplification factor must be positive. -}
{-# INLINE amplify #-}
amplify :: (Ring.C y, Dim.C u, Dim.C v) =>
      y
   -> Proc.T s u t (
        SigA.R s v y yv
     -> SigA.R s v y yv)
amplify volume = Proc.pure $ FiltV.amplify volume

{-# INLINE amplifyDimension #-}
amplifyDimension :: (Ring.C y, Dim.C u, Dim.C v0, Dim.C v1) =>
      DN.T v0 y
   -> Proc.T s u t (
        SigA.R s v1 y yv
     -> SigA.R s (Dim.Mul v0 v1) y yv)
amplifyDimension volume = Proc.pure $ FiltV.amplifyDimension volume


{-# INLINE negate #-}
negate :: (Additive.C yv, Dim.C u, Dim.C v) =>
      Proc.T s u t (
        SigA.R s v y yv
     -> SigA.R s v y yv)
negate = Proc.pure FiltV.negate


{-# INLINE envelope #-}
envelope :: (Flat.C flat y0, Ring.C y0, Dim.C u, Dim.C v) =>
      Proc.T s u t (
        RP.T s flat y0   {- v the envelope -}
     -> SigA.R s v y y0  {- v the signal to be enveloped -}
     -> SigA.R s v y y0)
envelope = Proc.pure FiltV.envelope

{-# INLINE envelopeVector #-}
envelopeVector :: (Flat.C flat y0, Module.C y0 yv, Ring.C y, Dim.C u, Dim.C v) =>
      Proc.T s u t (
        RP.T s flat y0   {- v the envelope -}
     -> SigA.R s v y yv  {- v the signal to be enveloped -}
     -> SigA.R s v y yv)
envelopeVector = Proc.pure FiltV.envelopeVector

{-# INLINE envelopeVectorDimension #-}
envelopeVectorDimension ::
   (Module.C y0 yv, Ring.C y, Dim.C u, Dim.C v0, Dim.C v1) =>
      Proc.T s u t (
        SigA.R s v0 y y0  {-  the envelope -}
     -> SigA.R s v1 y yv  {-  the signal to be enveloped -}
     -> SigA.R s (Dim.Mul v0 v1) y yv)
envelopeVectorDimension = Proc.pure FiltV.envelopeVectorDimension


{-# INLINE differentiate #-}
differentiate :: (Additive.C yv, Ring.C q, Dim.C u, Dim.C v) =>
      Proc.T s u q (
        SigA.R s v q yv
     -> SigA.R s (DimensionGradient u v) q yv)
differentiate =
   do rate <- Proc.getSampleRate
      return $ \ x ->
         SigA.fromSamples
            (rate &*& SigA.amplitude x)
            (FiltNR.differentiate (SigA.samples x))


{- | needs a good handling of boundaries, yet -}
{-# INLINE meanStatic #-}
meanStatic ::
   (RealField.C q, Module.C q yv, Dim.C u, Dim.C v) =>
      DN.T (Dim.Recip u) q   {- ^ cut-off freqeuncy -}
   -> Proc.T s u q (
        SigA.R s v q yv
     -> SigA.R s v q yv)
meanStatic time =
   FiltR.meanStatic time

meanStaticSeparateTY :: (Additive.C yv, Field.C y, RealField.C t,
         Module.C y yv, Dim.C u, Dim.C v) =>
      DN.T (Dim.Recip u) t   {- ^ cut-off freqeuncy -}
   -> Proc.T s u t (
        SigA.R s v y yv
     -> SigA.R s v y yv)
meanStaticSeparateTY time =
   -- FiltR.meanStatic time, means that 't' = 'y'
   do f <- toFrequencyScalar time
      return $ \ x ->
         let tInt  = round ((recip f - 1)/2)
             width = tInt*2+1
         in  SigA.processSamples
                ((SigA.asTypeOfAmplitude (recip (fromIntegral width)) x *> ) .
                 Delay.staticNeg tInt .
                 MA.sumsStaticInt width) x


{- | needs a better handling of boundaries, yet -}
{-# INLINE mean #-}
mean :: (Additive.C yv, RealField.C q,
         Module.C q yv, Dim.C u, Dim.C v, Sample.C q, Sample.C yv) =>
      DN.T (Dim.Recip u) q    {- ^ minimum cut-off freqeuncy -}
   -> Proc.T s u q (
        SigA.R s (Dim.Recip u) q q
                              {- v cut-off freqeuncies -}
     -> SigA.R s v q yv
     -> SigA.R s v q yv)
mean minFreq =
   FiltR.mean minFreq


{-# INLINE delay #-}
delay :: (Additive.C yv, Field.C y, RealField.C t, Dim.C u, Dim.C v) =>
      DN.T u t
   -> Proc.T s u t (
        SigA.R s v y yv
     -> SigA.R s v y yv)
delay time =
   do t <- toTimeScalar time
      return $ SigA.processSamples (Delay.static (round t))


{-# INLINE phaseModulation #-}
phaseModulation ::
   (Additive.C yv, RealField.C q, Dim.C u, Dim.C v,
    Sample.C q, Sample.C yv) =>
      Interpolation.T q yv
   -> DN.T u q
          {- ^ minDelay, minimal delay, may be negative -}
   -> DN.T u q
          {- ^ maxDelay, maximal delay, it must be @minDelay <= maxDelay@
               and the modulation must always be
               in the range [minDelay,maxDelay]. -}
   -> Proc.T s u q (
        SigA.R s u q q
          {- v delay control, positive numbers meanStatic delay,
               negative numbers meanStatic prefetch -}
     -> SigA.R s v q yv
     -> SigA.R s v q yv)
phaseModulation ip minDelay maxDelay =
   FiltR.phaseModulation ip minDelay maxDelay

{-# INLINE frequencyModulation #-}
frequencyModulation ::
   (Flat.C flat q, Additive.C yv, RealField.C q, Dim.C u, Dim.C v) =>
      Interpolation.T q yv
   -> Proc.T s u q (
        RP.T s flat q    {- v frequency factors -}
     -> SigA.R s v q yv
     -> SigA.R s v q yv)
frequencyModulation ip =
   Proc.pure $
      \ factors ->
          SigA.processSamples
             (FiltR.interpolateMultiRelativeZeroPad ip (Flat.toSamples factors))

{- |
Frequency modulation where the input signal can have a sample rate
different from the output.
(The sample rate values can differ, the unit must be the same.
We could lift that restriction,
but then the unit handling becomes more complicated,
and I didn't have a use for it so far.)

The function can be used for resampling.
-}
{-# INLINE frequencyModulationDecoupled #-}
frequencyModulationDecoupled ::
   (Flat.C flat q, Additive.C yv, RealField.C q, Dim.C u, Dim.C v) =>
      Interpolation.T q yv
   -> Proc.T s u q (
        RP.T s flat q    {- v frequency factors -}
     -> SigP.T u q (SigA.T v q (SigS.T Sig.T)) yv
     -> SigA.R s v q yv)
frequencyModulationDecoupled ip =
   fmap
      (\toFreq factors y ->
         flip SigA.processSamples (RP.fromSignal (SigP.signal y)) $
            FiltR.interpolateMultiRelativeZeroPad ip
               (SigA.scalarSamples toFreq
                  (SigA.fromSamples (SigP.sampleRate y) (Flat.toSamples factors))))
      (Proc.withParam Proc.toFrequencyScalar)


{- | symmetric phaser -}
{-# INLINE phaser #-}
phaser ::
   (Additive.C yv, RealField.C q,
    Module.C q yv, Dim.C u, Dim.C v,
    Sample.C q, Sample.C yv) =>
      Interpolation.T q yv
   -> DN.T u q  {- ^ maxDelay, must be positive -}
   -> Proc.T s u q (
        SigA.R s u q q
                {- v delay control -}
     -> SigA.R s v q yv
     -> SigA.R s v q yv)
phaser = FiltR.phaser
{-
phaser ip maxDelay =
   do p <- phaserCore ip maxDelay
      return $ \ delays x ->
         FiltV.amplify 0.5 .
         uncurry MiscV.mix . p delays $ x
-}

{-# INLINE phaserStereo #-}
phaserStereo ::
   (Additive.C yv, RealField.C q,
    Module.C q yv, Dim.C u, Dim.C v,
    Sample.C q, Sample.C yv) =>
      Interpolation.T q yv
   -> DN.T u q   {- ^ maxDelay, must be positive -}
   -> Proc.T s u q (
        SigA.R s u q q
                 {- v delay control -}
     -> SigA.R s v q yv
     -> SigA.R s v q (Stereo.T yv))
phaserStereo = FiltR.phaserStereo
{-
phaserStereo ip maxDelay =
   do p <- phaserCore ip maxDelay
      return $ \ delays -> uncurry CutV.zip . p delays
-}

{-
{-# INLINE phaserCore #-}
phaserCore ::
   (Additive.C yv, RealField.C q,
    Module.C q yv, Dim.C u, Dim.C v,
    Sample.C q, Sample.C yv) =>
      Interpolation.T q yv
   -> DN.T u q   {- ^ maxDelay, must be positive -}
   -> Proc.T s u q (
        SigA.R s u q q
                 {- v delay control -}
     -> SigA.R s v q yv
     -> (SigA.R s v q yv, SigA.R s v q yv))
phaserCore ip maxDelay =
   do let minDelay  = Additive.negate maxDelay
      pm <- phaseModulation ip minDelay maxDelay
      return $ \ delays x ->
         let negDelays = FiltV.negate delays
         in  (pm delays x,
              pm negDelays x)
-}



{-# INLINE firstOrderLowpass #-}
{-# INLINE firstOrderHighpass #-}
firstOrderLowpass, firstOrderHighpass ::
   (Trans.C q, Module.C q yv, Dim.C u, Dim.C v) =>
      CProc.T s u q
          (SigA.R s (Dim.Recip u) q q
                    {- v Control signal for the cut-off frequency. -} )
          (Filt1.Parameter q) (
        SigA.R s v q yv
                    {- v Input signal -}
     -> SigA.R s v q yv)
firstOrderLowpass  = firstOrderGen Filt1.lowpassModifier
firstOrderHighpass = firstOrderGen Filt1.highpassModifier

{-# INLINE firstOrderGen #-}
firstOrderGen ::
   (Trans.C q, Module.C q yv, Dim.C u, Dim.C v) =>
      (Modifier yv (Filt1.Parameter q) yv yv)
--      (Sig.T (Filt1.Parameter q) -> Sig.T yv -> Sig.T yv)
   -> CProc.T s u q (SigA.R s (Dim.Recip u) q q) (Filt1.Parameter q) (
        SigA.R s v q yv
     -> SigA.R s v q yv)
firstOrderGen modif =
   frequencyControl Filt1.parameter
      (Sig.modifyModulated modif)



{-# INLINE butterworthLowpass #-}
{-# INLINE butterworthHighpass #-}
{-# INLINE chebyshevALowpass #-}
{-# INLINE chebyshevAHighpass #-}
{-# INLINE chebyshevBLowpass #-}
{-# INLINE chebyshevBHighpass #-}

butterworthLowpass, butterworthHighpass,
   chebyshevALowpass, chebyshevAHighpass,
   chebyshevBLowpass, chebyshevBHighpass ::
      (Trans.C q, VectorSpace.C q yv, Dim.C u, Dim.C v) =>
      NonNeg.Int   {- ^ Order of the filter, must be even,
                        the higher the order, the sharper is the separation of frequencies. -}
   -> q            {- ^ The attenuation at the cut-off frequency.
                        Should be between 0 and 1. -}
   -> CProc.T s u q
         (SigA.R s (Dim.Recip u) q q
                      {- v Control signal for the cut-off frequency. -}  )
         q (
        SigA.R s v q yv {- v Input signal -}
     -> SigA.R s v q yv)

butterworthLowpass  = higherOrderNoResoGen Butter.lowpass
butterworthHighpass = higherOrderNoResoGen Butter.highpass
chebyshevALowpass   = higherOrderNoResoGen Cheby.lowpassA
chebyshevAHighpass  = higherOrderNoResoGen Cheby.highpassA
chebyshevBLowpass   = higherOrderNoResoGen Cheby.lowpassB
chebyshevBHighpass  = higherOrderNoResoGen Cheby.highpassB

{- FIXME:
currently only frequencies can be interpolated not the filter parameters,
this is not very efficient
-}
{- TODO:
initial value
-}
{-# INLINE higherOrderNoResoGen #-}
higherOrderNoResoGen ::
   (Field.C q, Dim.C u, Dim.C v) =>
      (Int -> q -> [q] -> [yv] -> [yv])
   -> NonNeg.Int
   -> q
   -> CProc.T s u q (SigA.R s (Dim.Recip u) q q) q (
        SigA.R s v q yv
     -> SigA.R s v q yv)
higherOrderNoResoGen filt order ratio =
   frequencyControl id
      (\ cs xs ->
           Sig.fromList (filt (NonNeg.toNumber order) ratio
               (Sig.toList cs) (Sig.toList xs)))



{-# INLINE universal #-}
universal ::
   (Flat.C flat q, Trans.C q, Module.C q yv, Dim.C u, Dim.C v) =>
      CProc.T s u q
         (RP.T s flat q
                   {- v signal for resonance,
                        i.e. factor of amplification at the resonance frequency
                        relatively to the transition band. -},
          SigA.R s (Dim.Recip u) q q
                   {- v signal for cut off and band center frequency -} )
         (UniFilter.Parameter q) (
        SigA.R s v q yv
                    {- v input signal -}
     -> SigA.R s v q (UniFilter.Result yv))
                    {- ^ highpass, bandpass, lowpass filter -}
universal =
   frequencyResonanceControl
      UniFilter.parameter
      (Sig.modifyModulated UniFilter.modifier)

{-# INLINE moogLowpass #-}
moogLowpass :: (Flat.C flat q, Trans.C q, Module.C q yv, Dim.C u, Dim.C v) =>
      NonNeg.Int
   -> CProc.T s u q
         (RP.T s flat q
                   {- v signal for resonance,
                        i.e. factor of amplification at the resonance frequency
                        relatively to the transition band. -},
          SigA.R s (Dim.Recip u) q q
                   {- v signal for cut off frequency -} )
         (Moog.Parameter q) (
        SigA.R s v q yv
     -> SigA.R s v q yv)
moogLowpass order =
   let orderInt = NonNeg.toNumber order
   in  frequencyResonanceControl
          (Moog.parameter orderInt)
          (Sig.modifyModulated (Moog.lowpassModifier orderInt))

{-# INLINE allpassCascade #-}
allpassCascade :: (Trans.C q, Module.C q yv, Dim.C u, Dim.C v) =>
      NonNeg.Int  {- ^ order, number of filters in the cascade -}
   -> q           {- ^ the phase shift to be achieved for the given frequency -}
   -> CProc.T s u q
         (SigA.R s (Dim.Recip u) q q {- v lowest comb frequency -})
         (Allpass.Parameter q) (
        SigA.R s v q yv
     -> SigA.R s v q yv)
allpassCascade order phase =
   let orderInt = NonNeg.toNumber order
   in  frequencyControl
          (Allpass.parameter orderInt phase)
          (Sig.modifyModulated (Allpass.cascadeModifier orderInt))


{-# INLINE frequencyControl #-}
frequencyControl ::
   (Field.C y, Dim.C u, Dim.C v) =>
   (y -> ic) ->
   (Sig.T ic -> Sig.T yv0 -> Sig.T yv1) ->
   CProc.T s u y
      (SigA.R s (Dim.Recip u) y y) ic
      (SigA.R s v y1 yv0 -> SigA.R s v y1 yv1)

frequencyControl mkParam filt = CProc.Cons $
   do toFreq <- Proc.withParam toFrequencyScalar
      return
         (\ freqs -> Sig.map mkParam (SigA.scalarSamples toFreq freqs),
          \ params -> SigA.processSamples (filt params))


{-# INLINE frequencyResonanceControl #-}
frequencyResonanceControl ::
   (Flat.C flat y, Field.C y, Dim.C u, Dim.C v) =>
   (FiltR.Pole y -> ic) ->
   (Sig.T ic -> Sig.T yv0 -> Sig.T yv1) ->
   CProc.T s u y
      (RP.T s flat y, SigA.R s (Dim.Recip u) y y) ic
      (SigA.R s v y1 yv0 -> SigA.R s v y1 yv1)

frequencyResonanceControl mkParam filt = CProc.Cons $
   do toFreq <- Proc.withParam toFrequencyScalar
      return
         (\ (resos, freqs) ->
               Sig.map mkParam $
               Sig.zipWith FiltR.Pole
                  (Flat.toSamples resos)
                  (SigA.scalarSamples toFreq freqs),
          \ params -> SigA.processSamples (filt params))


{- | Infinitely many equi-delayed exponentially decaying echos. -}
{-# INLINE comb #-}
comb :: (RealField.C t, Module.C y yv, Dim.C u, Dim.C v, Sample.C yv) =>
   DN.T u t -> y -> Proc.T s u t (SigA.R s v y yv -> SigA.R s v y yv)
comb = FiltR.comb


{- | Infinitely many equi-delayed echos processed by an arbitrary time-preserving signal processor. -}
{-# INLINE combProc #-}
combProc ::
   (RealField.C t, Real.C y, Field.C y, Module.C y yv, Sample.C yv,
    Dim.C u, Dim.C v) =>
   DN.T u t ->
   Proc.T s u t (SigA.R s v y yv -> SigA.R s v y yv) ->
   Proc.T s u t (SigA.R s v y yv -> SigA.R s v y yv)
combProc time proc =
   do f <- proc
      t <- fmap round $ toTimeScalar time
      let chunkSize = SigSt.chunkSize t
      return $ \x ->
         SigA.processSamples
            (Sig.fromStorableSignal .
             Comb.runProc t
                (Sig.toStorableSignal chunkSize .
                 SigA.vectorSamples (SigA.toAmplitudeScalar x) .
                 f .
                 SigA.fromSamples (SigA.amplitude x) .
                 Sig.fromStorableSignal) .
             Sig.toStorableSignal chunkSize) x

{-
combProc time proc sr x =
   Rate.loop (\sr' y -> MiscV.mixVolume (SigA.amplitude x) x (delay time sr' (proc sr' y))) sr
-}


{-# INLINE integrate #-}
integrate :: (Additive.C yv, Field.C q, Dim.C u, Dim.C v) =>
      Proc.T s u q (
        SigA.R s v q yv
     -> SigA.R s (Dim.Mul u v) q yv)
integrate =
   do rate <- Proc.getSampleRate
      return $ \ x ->
         SigA.replaceAmplitude
            (DN.rewriteDimension (Dim.commute . Dim.applyRightMul Dim.invertRecip) $
             SigA.amplitude x &/& rate)
            (Hom.processSamples Integrate.run x)