{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Copyright   :  (c) Henning Thielemann 2006
License     :  GPL

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

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

   {- ** Smooth -}
   mean,

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


   {- * Recursive -}

   {- ** Without resonance -}
   firstOrderLowpass,
   firstOrderHighpass,
   butterworthLowpass,
   butterworthHighpass,
   chebyshevALowpass,
   chebyshevAHighpass,
   chebyshevBLowpass,
   chebyshevBHighpass,
   {- ** With resonance -}
   universal,
   moogLowpass,
   {- ** Allpass -}
   allpassCascade,
   {- ** Reverb -}
   comb,
   {- ** Filter operators from calculus -}
   integrate
) where


import qualified UniqueLogicNP.Explicit.Process    as Process
import qualified UniqueLogicNP.Explicit.Expression as Expr
import qualified Synthesizer.Inference.Monad.Signal               as SigI
import qualified Synthesizer.Inference.Monad.Signal.Cut           as CutI
import qualified Synthesizer.Inference.Monad.Signal.Displacement as SynI

import UniqueLogicNP.Explicit.Expression ((=!=))
import Synthesizer.Inference.Monad.Signal
   (toTimeScalar, toFrequencyScalar, sampleRateExpr,
    amplitudeExpr)

import qualified Synthesizer.Physical.Signal as SigP
import qualified Synthesizer.Plain.Signal as Sig
import qualified Synthesizer.Plain.Displacement as Syn
import qualified Synthesizer.Plain.Interpolation as Interpolation
import qualified Synthesizer.Plain.Filter.Delay.ST 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.Comb        as Comb
import qualified Synthesizer.Plain.Filter.Recursive.Integration as Integrate
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.Plain.Filter.NonRecursive as FiltNR

import qualified Algebra.OccasionallyScalar as OccScalar
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.Module         as Module
import qualified Algebra.VectorSpace    as VectorSpace

import Control.Monad (liftM2, )

import Data.Ord.HT (limit, )

import NumericPrelude hiding (negate)
import PreludeBase as P


{- | The amplification factor must be positive. -}
amplify :: (Field.C q) =>
      q
   -> SigI.T a q v
   -> SigI.Process a q v
amplify volume x =
   do amplitude <- Process.fromExpr (Expr.constant volume * amplitudeExpr x)
      SigI.returnCons (SigP.sampleRate x) amplitude (SigP.samples x)

negate :: (Additive.C v, Eq q) =>
      SigI.T a q v
   -> SigI.Process a q v
negate =
   returnModified [] Additive.negate


envelope :: (Module.C y v, Field.C q, Eq q) =>
      SigI.T a q y  {- ^ the envelope -}
   -> SigI.T a q v  {- ^ the signal to be enveloped -}
   -> SigI.Process a q v
envelope y x =
   do sampleRate <- Process.fromExpr (sampleRateExpr x =!= sampleRateExpr y)
      amplitude  <- Process.fromExpr (amplitudeExpr  x  *  amplitudeExpr  y)
      SigI.returnCons sampleRate amplitude
         (FiltNR.envelopeVector (SigP.samples y) (SigP.samples x))


differentiate :: (Additive.C v, Field.C q, Eq q) =>
      SigI.T a q v
   -> SigI.Process a q v
differentiate x =
   do amp <- Process.fromExpr (amplitudeExpr x * sampleRateExpr x)
      SigI.returnCons
         (SigP.sampleRate x) amp
         (FiltNR.differentiate (SigP.samples x))


{- | needs a good handling of boundaries, yet -}
mean :: (Additive.C v, Field.C q, Eq q, RealField.C a,
         Module.C a v, OccScalar.C a q) =>
      q            {- ^ time length of the window -}
   -> SigI.T a q v
   -> SigI.Process a q v
mean time x =
   do t <- toTimeScalar x (Expr.constant time)
      let tInt  = round ((t-1)/2)
      let width = tInt*2+1
      returnModified []
         ((SigP.asTypeOfAmplitude (recip (fromIntegral width)) x *> ) .
          FiltNR.sums width . FiltNR.delay tInt) x


delay :: (Additive.C v, Field.C q, Eq q, RealField.C a, OccScalar.C a q) =>
      q
   -> SigI.T a q v
   -> SigI.Process a q v
delay time x =
   do t <- toTimeScalar x (Expr.constant time)
      returnModified [] (FiltNR.delay (round t)) x


phaseModulation ::
         (Additive.C v, Field.C q, Eq q, RealField.C a, OccScalar.C a q) =>
      Interpolation.T a v
   -> q   {- ^ minDelay, minimal delay, may be negative -}
   -> q   {- ^ maxDelay, maximal delay, it must be @minDelay <= maxDelay@
               and the modulation must always be
               in the range [minDelay,maxDelay]. -}
   -> SigI.T a q a
          {- ^ delay control, positive numbers mean delay,
               negative numbers mean prefetch -}
   -> SigI.T a q v
   -> SigI.Process a q v
phaseModulation ip minDelay maxDelay delays x =
   do t0 <- toTimeScalar x (Expr.constant minDelay)
      t1 <- toTimeScalar x (Expr.constant maxDelay)
      let tInt0 = floor   t0
      let tInt1 = ceiling t1
      let tInt0Neg = Additive.negate tInt0
      ds <- SigI.scalarSamples (toTimeScalar delays) delays
      returnModified [SigP.sampleRate delays]
         (FiltNR.delay tInt0 .
             Delay.modulated ip (tInt1-tInt0+1)
               (FiltNR.delay tInt0Neg
                  (Syn.raise (fromIntegral tInt0Neg)
                     (map (limit (t0,t1)) ds)))) x


{- | symmetric phaser -}
phaser :: (Additive.C v, Field.C q, Eq q, RealField.C a,
           Module.C a v, OccScalar.C a q) =>
      Interpolation.T a v
   -> q   {- ^ maxDelay, must be positive -}
   -> SigI.T a q a
          {- ^ delay control -}
   -> SigI.T a q v
   -> SigI.Process a q v
phaser ip maxDelay delays x =
   amplify (asTypeOf 0.5 maxDelay) =<<
      uncurry SynI.mix =<< phaserCore ip maxDelay delays x

phaserStereo :: (Additive.C v, Field.C q, Eq q, Real.C q, RealField.C a,
                 Module.C a v, OccScalar.C a q) =>
      Interpolation.T a v
   -> q   {- ^ maxDelay, must be positive -}
   -> SigI.T a q a
          {- ^ delay control -}
   -> SigI.T a q v
   -> SigI.Process a q (v,v)
phaserStereo ip maxDelay delays x =
   uncurry CutI.zip =<< phaserCore ip maxDelay delays x

phaserCore :: (Additive.C v, Field.C q, Eq q, RealField.C a,
               Module.C a v, OccScalar.C a q) =>
      Interpolation.T a v
   -> q   {- ^ maxDelay, must be positive -}
   -> SigI.T a q a
          {- ^ delay control -}
   -> SigI.T a q v
   -> Process.T q (SigI.T a q v, SigI.T a q v)
phaserCore ip maxDelay delays x =
   do let minDelay = Additive.negate maxDelay
      negDelays <- negate delays  -- FiltI.negate delays
      liftM2 (,)
         (phaseModulation ip minDelay maxDelay delays x)
         (phaseModulation ip minDelay maxDelay negDelays x)



firstOrderLowpass, firstOrderHighpass ::
   (Trans.C a, Trans.C q, Eq q, Module.C a v, OccScalar.C a q) =>
      SigI.T a q a {- ^ Control signal for the cut-off frequency. -}
   -> SigI.T a q v {- ^ Input signal -}
   -> SigI.Process a q v
firstOrderLowpass  = firstOrderGen Filt1.lowpass
firstOrderHighpass = firstOrderGen Filt1.highpass

firstOrderGen :: (Trans.C a, Trans.C q, Eq q, Module.C a v, OccScalar.C a q) =>
      (Sig.T (Filt1.Parameter a) -> Sig.T v -> Sig.T v)
   -> SigI.T a q a
   -> SigI.T a q v
   -> SigI.Process a q v
firstOrderGen filt freq x =
   do freqs <- SigI.scalarSamples (toFrequencyScalar x) freq
      returnModified [SigP.sampleRate freq]
         (filt (map Filt1.parameter freqs)) x


butterworthLowpass, butterworthHighpass,
   chebyshevALowpass, chebyshevAHighpass,
   chebyshevBLowpass, chebyshevBHighpass ::
      (Field.C q, Eq q, Trans.C a, VectorSpace.C a v, OccScalar.C a q) =>
      Int          {- ^ Order of the filter, must be even,
                        the higher the order, the sharper is the separation of frequencies. -}
   -> SigI.T a q a {- ^ The attenuation at the cut-off frequency.
                        Should be between 0 and 1. -}
   -> SigI.T a q a {- ^ Control signal for the cut-off frequency. -}
   -> SigI.T a q v {- ^ Input signal -}
   -> SigI.Process a q v

butterworthLowpass  = higherOrderNoResoGen Butter.lowpassPole
butterworthHighpass = higherOrderNoResoGen Butter.highpassPole
chebyshevALowpass   = higherOrderNoResoGen Cheby.lowpassAPole
chebyshevAHighpass  = higherOrderNoResoGen Cheby.highpassAPole
chebyshevBLowpass   = higherOrderNoResoGen Cheby.lowpassBPole
chebyshevBHighpass  = higherOrderNoResoGen Cheby.highpassBPole

higherOrderNoResoGen ::
   (Field.C q, Eq q, Ring.C a, OccScalar.C a q) =>
      (Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v)
   -> Int
   -> SigI.T a q a
   -> SigI.T a q a
   -> SigI.T a q v
   -> SigI.Process a q v
higherOrderNoResoGen filt order ratio freq x =
   do ratios <- SigI.scalarSamples (Process.exprToScalar) ratio
      freqs  <- SigI.scalarSamples (toFrequencyScalar x) freq
      returnModified [SigP.sampleRate freq]
         (filt order ratios freqs) x



universal :: (Trans.C a, Module.C a v, Field.C q, Eq q, OccScalar.C a q) =>
      SigI.T a q a {- ^ signal for resonance,
                        i.e. factor of amplification at the resonance frequency
                        relatively to the transition band. -}
   -> SigI.T a q a {- ^ signal for cut off and band center frequency -}
   -> SigI.T a q v {- ^ input signal -}
   -> SigI.Process a q (UniFilter.Result v)
                   {- ^ highpass, bandpass, lowpass filter -}
universal reso freq x =
   do resos <- SigI.scalarSamples (Process.exprToScalar) reso
      freqs <- SigI.scalarSamples (toFrequencyScalar x) freq
      let params =
             map UniFilter.parameter
                 (zipWith FiltR.Pole resos freqs)
      returnModified [SigP.sampleRate reso, SigP.sampleRate freq]
         (UniFilter.run params) x

moogLowpass :: (Trans.C a, Module.C a v, Field.C q, Eq q, OccScalar.C a q) =>
      Int
   -> SigI.T a q a {- ^ signal for resonance,
                        i.e. factor of amplification at the resonance frequency
                        relatively to the transition band. -}
   -> SigI.T a q a {- ^ signal for cut off and band center frequency -}
   -> SigI.T a q v
   -> SigI.Process a q v
moogLowpass order reso freq x =
   do resos <- SigI.scalarSamples (Process.exprToScalar) reso
      freqs <- SigI.scalarSamples (toFrequencyScalar x) freq
      let params =
             map (Moog.parameter order)
                 (zipWith FiltR.Pole resos freqs)
      returnModified [SigP.sampleRate reso, SigP.sampleRate freq]
         (Moog.lowpass order params) x

allpassCascade :: (Trans.C a, Module.C a v, Field.C q, Eq q, OccScalar.C a q) =>
      Int          {- ^ order, number of filters in the cascade -}
   -> a            {- ^ the phase shift to be achieved for the given frequency -}
   -> SigI.T a q a {- ^ lowest comb frequency -}
   -> SigI.T a q v
   -> SigI.Process a q v
allpassCascade order phase freq x =
   do freqs <- SigI.scalarSamples (toFrequencyScalar x) freq
      let params = map (Allpass.parameter order phase) freqs
      returnModified [SigP.sampleRate freq]
         (Allpass.cascade order params) x



{- | Infinitely many equi-delayed exponentially decaying echos. -}
comb :: (RealField.C a, Field.C q, Eq q, OccScalar.C a q, Module.C a v) =>
   q -> a -> SigI.T a q v -> SigI.Process a q v
comb time gain x =
   do t <- toTimeScalar x (Expr.constant time)
      returnModified [] (Comb.run (round t) gain) x


integrate :: (Additive.C v, Field.C q, Eq q) =>
      SigI.T a q v -> SigI.Process a q v
integrate x =
   do amp <- Process.fromExpr (amplitudeExpr x / sampleRateExpr x)
      SigI.returnCons
         (SigP.sampleRate x) amp
         (Integrate.run (SigP.samples x))


returnModified :: (Eq q) =>
   [Process.Atom q] -> (Sig.T v -> Sig.T w) -> SigI.T a q v -> SigI.Process a q w
returnModified sampleRates proc x =
   do let sampleRate = SigP.sampleRate x
      mapM_ (Process.equalValue sampleRate) sampleRates
      SigI.returnCons
         sampleRate (SigP.amplitude x)
         (proc (SigP.samples x))