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

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

   {- ** Amplification -}
   negate,
   envelope,
   envelopeVector,
   convolveVector,

   {- ** Smooth -}
   mean,
   meanStatic,

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


   {- * Recursive -}

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

   {- ** Allpass -}
   allpassCascade,
   allpassFlangerPhase,

   {- ** Reverb -}
   comb,

   {- * Helper functions -}
   interpolateMultiRelativeZeroPad,
) where

-- import qualified Synthesizer.Dimensional.Abstraction.Linear as Lin
import qualified Synthesizer.Dimensional.Abstraction.Homogeneous as Hom
import qualified Synthesizer.Dimensional.Abstraction.RateIndependent as Ind
import qualified Synthesizer.Dimensional.Abstraction.Flat as Flat

import qualified Synthesizer.Dimensional.RatePhantom as RP

import qualified Synthesizer.Dimensional.Amplitude.Filter       as FiltV
import qualified Synthesizer.Dimensional.Process as Proc
-- import qualified Synthesizer.Dimensional.Rate as Rate

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

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.State.Signal as Sig
import Synthesizer.Plain.Signal (Modifier, )

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

import qualified Synthesizer.Causal.Process       as Causal
import qualified Synthesizer.Causal.Interpolation as Interpolation
import qualified Synthesizer.State.Displacement as Disp
import qualified Synthesizer.State.Filter.Delay as Delay
import qualified Synthesizer.State.Filter.Recursive.MovingAverage as MA
import qualified Synthesizer.State.Filter.NonRecursive as FiltNR

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 FiltRec

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

-- import qualified Synthesizer.Generic.Interpolation as InterpolationG
import qualified Synthesizer.Generic.Filter.Recursive.MovingAverage as MAG
import qualified Synthesizer.Generic.Filter.NonRecursive as FiltG
import qualified Synthesizer.Generic.Filter.Delay as DelayG
import qualified Synthesizer.Generic.Signal  as SigG
import qualified Synthesizer.Generic.Signal2 as SigG2

import qualified Synthesizer.Frame.Stereo as Stereo

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

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 Foreign.Storable (Storable, )

-- import qualified Data.List as List

-- import Control.Monad(liftM2)

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


{-# INLINE negate #-}
negate :: (Hom.C sig, Additive.C yv, Dim.C u) =>
      Proc.T s u t (
        RP.T s sig yv
     -> RP.T s sig yv)
negate = Proc.pure FiltV.negate


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

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

{-# INLINE convolveVector #-}
convolveVector ::
   (Hom.C sig, Module.C q yv, Field.C q, Dim.C u) =>
      Proc.T s u q (
        SigA.R s (Dim.Recip u) q q
                              {- v the filter window -}
     -> RP.T s sig yv         {- v the signal to be enveloped -}
     -> RP.T s sig yv)
convolveVector =
   do toFreq <- Proc.withParam toFrequencyScalar
      return $ \ window ->
         Hom.processSamples
            (FiltNR.generic (SigA.scalarSamples toFreq window))


{- | needs a better handling of boundaries, yet -}
{-# INLINE meanStatic #-}
meanStatic :: (Hom.C sig, Additive.C yv, RealField.C q,
         Module.C q yv, Dim.C u) =>
      DN.T (Dim.Recip u) q    {- ^ cut-off freqeuncy -}
   -> Proc.T s u q (
        RP.T s sig yv
     -> RP.T s sig yv)
meanStatic freq =
   do f <- toFrequencyScalar freq
      return $
         let tInt  = round ((recip f - 1)/2)
             width = tInt*2+1
         in  Hom.processSamples
                ((asTypeOf (recip (fromIntegral width)) f *> ) .
                 Delay.staticNeg tInt .
                 MA.sumsStaticInt width)

{- | needs a better handling of boundaries, yet -}
{-# INLINE mean #-}
mean :: (Hom.C sig, Additive.C yv, RealField.C q,
         Module.C q yv, Dim.C u, Storable q, Storable 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 -}
     -> RP.T s sig yv
     -> RP.T s sig yv)
mean minFreq =
   do mf <- toFrequencyScalar minFreq
      frequencyControl $ \ freqs ->
         let tMax   = ceiling (recip (2*mf))
             err    = error "Filter.mean: frequencies must be positive"
             widths = Sig.map (\f -> if f>0 then recip (2*f) else err) freqs
         in  Hom.processSamples
                (fromStorable .
--                 MAG.sumsStaticInt tMax .
                 MAG.modulatedFrac tMax (toStorable widths) .
                 toStorable)

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


{-# INLINE toStorable #-}
toStorable :: (Storable a) => Sig.T a -> SigSt.T a
toStorable = Sig.toStorableSignal SigSt.defaultChunkSize

{-# INLINE fromStorable #-}
fromStorable :: (Storable a) => SigSt.T a -> Sig.T a
fromStorable = Sig.fromStorableSignal

{-# INLINE phaseModulation #-}
phaseModulation ::
   (Hom.C sig, Additive.C yv, RealField.C q, Dim.C u,
    Storable q, Storable yv) =>
      Interpolation.T q yv
   -> DN.T u q
          {- ^ minimal deviation from current time, usually negative -}
   -> DN.T u q
          {- ^ maximal deviation, it must be @minDev <= maxDev@
               and the modulation must always be
               in the range [minDev,maxDev]. -}
   -> Proc.T s u q (
        SigA.R s u q q
          {- v deviation control,
               positive numbers meanStatic prefetch,
               negative numbers meanStatic delay -}
     -> RP.T s sig yv
     -> RP.T s sig yv)
phaseModulation ip minDev maxDev =
   fmap
      (\f devs ->
         Hom.processSamples
            (Sig.fromStorableSignal .
             f (SigA.processSamples toStorable devs) .
             toStorable))
      (phaseModulationGeneric ip minDev maxDev)

{-# INLINE phaseModulationGeneric #-}
phaseModulationGeneric ::
   (Additive.C yv, RealField.C q, Dim.C u,
    SigG2.Transform sig q yv, SigG.Write sig yv) =>
      Interpolation.T q yv
   -> DN.T u q
          {- ^ minimal deviation from current time, usually negative -}
   -> DN.T u q
          {- ^ maximal deviation, it must be @minDev <= maxDev@
               and the modulation must always be
               in the range [minDev,maxDev]. -}
   -> Proc.T s u q (
        RP.T s (SigA.D u q (SigS.T sig)) q
          {- v deviation control,
               positive numbers meanStatic prefetch,
               negative numbers meanStatic delay -}
     -> sig yv
     -> sig yv)
phaseModulationGeneric ip minDev _maxDev =
   fmap
      (\toTime devs ->
          let t0    = toTime minDev
              tInt0 = floor t0
          in  DelayG.modulated ip tInt0
                 (SigG.map (max t0) (SigA.scalarSamplesGeneric toTime devs)))
      (Proc.withParam toTimeScalar)


{-
FIXME: move to Dimensional.Straight
-}
{-# INLINE frequencyModulation #-}
frequencyModulation ::
   (Hom.C sig, Flat.C flat t,
    Additive.C yv, RealField.C t, Dim.C u) =>
      Interpolation.T t yv
   -> Proc.T s u t (
        RP.T s flat t    {- v frequency factors -}
     -> RP.T s sig yv
     -> RP.T s sig yv)
frequencyModulation ip =
   Proc.pure $
      \ factors ->
          Hom.processSamples
             (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 ::
   (Hom.C sig, Flat.C flat t,
    Additive.C yv, RealField.C t, Dim.C u) =>
      Interpolation.T t yv
   -> SigP.T u t sig yv
                   {- ToDo: We could also allow any signal from Generic.Read class. -}
   -> Proc.T s u t (
        RP.T s flat t {- v frequency factors -}
     -> RP.T s sig yv)
frequencyModulationDecoupled ip y =
   fmap
      (\toFreq factors ->
         RP.fromSignal $
         flip Hom.unwrappedProcessSamples (SigP.signal y) $
         interpolateMultiRelativeZeroPad ip $
         SigA.scalarSamples toFreq $
         SigA.fromSamples (SigP.sampleRate y) $
         Flat.toSamples factors)
      (Proc.withParam Proc.toFrequencyScalar)



{-# INLINE interpolateMultiRelativeZeroPad #-}
interpolateMultiRelativeZeroPad ::
    (RealField.C q, Additive.C yv) =>
    Interpolation.T q yv
    -> Sig.T q
    -> Sig.T yv
    -> Sig.T yv
interpolateMultiRelativeZeroPad ip k x =
    Causal.apply (Interpolation.relativeZeroPad zero ip zero x) k

{- | symmetric phaser -}
{-# INLINE phaser #-}
phaser ::
   (Hom.C sig, Additive.C yv, RealField.C q,
    Module.C q yv, Dim.C u,
    Storable q, Storable yv) =>
      Interpolation.T q yv
   -> DN.T u q  {- ^ maxDev, must be positive -}
   -> Proc.T s u q (
        SigA.R s u q q
                {- v delay control -}
     -> RP.T s sig yv
     -> RP.T s sig yv)
phaser ip maxDev =
   fmap
      (\p devs ->
         Hom.processSamples
            (FiltNR.amplifyVector (SigA.asTypeOfAmplitude 0.5 devs) .
             uncurry Disp.mix . p devs))
      (phaserCore ip maxDev)

{-# INLINE phaserStereo #-}
phaserStereo ::
   (Hom.C sig, Additive.C yv, RealField.C q,
    Module.C q yv, Dim.C u,
    Storable q, Storable yv) =>
      Interpolation.T q yv
   -> DN.T u q   {- ^ maxDev, must be positive -}
   -> Proc.T s u q (
        SigA.R s u q q
                 {- v delay control -}
     -> RP.T s sig yv
     -> RP.T s sig (Stereo.T yv))
phaserStereo ip maxDev =
   fmap
      (\p devs ->
            Hom.processSamples (uncurry (Sig.zipWith Stereo.cons) . p devs))
      (phaserCore ip maxDev)

{-# INLINE phaserCore #-}
phaserCore ::
   (Additive.C yv, RealField.C q,
    Module.C q yv, Dim.C u,
    Storable q, Storable yv) =>
      Interpolation.T q yv
   -> DN.T u q   {- ^ maxDev, must be positive -}
   -> Proc.T s u q (
        SigA.R s u q q
                 {- v delay control -}
     -> Sig.T yv
     -> (Sig.T yv, Sig.T yv))
phaserCore ip maxDev =
   do let minDev  = Additive.negate maxDev
      pm <- phaseModulationGeneric ip minDev maxDev
      return $ \ devs x ->
         let devsPos = SigA.processSamples toStorable devs
             devsNeg = SigA.processSamples FiltG.negate devsPos
             xst     = toStorable x
         in  (fromStorable (pm devsPos xst),
              fromStorable (pm devsNeg xst))


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

{-# INLINE firstOrderGen #-}
firstOrderGen ::
   (Hom.C sig, Trans.C q, Module.C q yv, Dim.C u) =>
      (Modifier yv (Filt1.Parameter q) yv yv)
   -> Proc.T s u q (
        SigA.R s (Dim.Recip u) q q
     -> RP.T s sig yv
     -> RP.T s sig yv)
firstOrderGen modif =
   frequencyControl $ \ freqs ->
      modifyModulated Filt1.parameter modif freqs


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

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

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


{-# INLINE higherOrderNoResoGen #-}
higherOrderNoResoGen ::
   (Hom.C sig, Flat.C flat q, Field.C q, Dim.C u) =>
      (Int -> [q] -> [q] -> [yv] -> [yv])
   -> NonNeg.Int
   -> Proc.T s u q (
        RP.T s flat q
     -> SigA.R s (Dim.Recip u) q q
     -> RP.T s sig yv
     -> RP.T s sig yv)
higherOrderNoResoGen filt order =
   fmap flip $ frequencyControl $ \ freqs ratios ->
      Hom.processSampleList
         (filt (NonNeg.toNumber order) (Sig.toList (Flat.toSamples ratios)) (Sig.toList freqs))



{-# INLINE highpassFromUniversal #-}
{-# INLINE bandpassFromUniversal #-}
{-# INLINE lowpassFromUniversal #-}
{-# INLINE bandlimitFromUniversal #-}
highpassFromUniversal, lowpassFromUniversal,
  bandpassFromUniversal, bandlimitFromUniversal ::
   (Hom.C sig) =>
        RP.T s sig (UniFilter.Result yv)
     -> RP.T s sig yv
{-
   (Hom.C sig, Dim.C u) =>
      Proc.T s u q (
        RP.T s sig (UniFilter.Result yv)
     -> RP.T s sig yv)
-}
highpassFromUniversal  = homogeneousMap UniFilter.highpass
bandpassFromUniversal  = homogeneousMap UniFilter.bandpass
lowpassFromUniversal   = homogeneousMap UniFilter.lowpass
bandlimitFromUniversal = homogeneousMap UniFilter.bandlimit

homogeneousMap ::
   (Hom.C sig, Ind.C w) =>
   (y0 -> y1) ->
   w sig y0 -> w sig y1
homogeneousMap f =
   Ind.processSignal (Hom.unwrappedProcessSamples (Sig.map f))

{-
homogeneousMap0 :: (Hom.C sig) =>
   (y0 -> y1) ->
   RP.T s sig y0 -> RP.T s sig y1
homogeneousMap0 f =
   Hom.processSamples (Sig.map f)

homogeneousMap1 :: (Hom.C sig) =>
   (y0 -> y1) ->
   Proc.T s1 u t (RP.T s sig y0 -> RP.T s sig y1)
homogeneousMap1 f =
   Proc.pure (Hom.processSamples (Sig.map f))
-}


{-# INLINE universal #-}
universal ::
   (Hom.C sig, Flat.C flat q, Trans.C q, Module.C q yv, Dim.C u) =>
      Proc.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 -}
     -> RP.T s sig yv
                    {- v input signal -}
     -> RP.T s sig (UniFilter.Result yv))
                    {- ^ highpass, bandpass, lowpass filter -}
universal =
   fmap flip $ frequencyControl $ \ freqs reso ->
      let resos = Flat.toSamples reso
      in  modifyModulated
             UniFilter.parameter
             UniFilter.modifier
             (Sig.zipWith FiltRec.Pole resos freqs)

{-# INLINE moogLowpass #-}
moogLowpass :: (Hom.C sig, Flat.C flat q, Trans.C q, Module.C q yv, Dim.C u) =>
      NonNeg.Int
   -> Proc.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 -}
     -> RP.T s sig yv
     -> RP.T s sig yv)
moogLowpass order =
   fmap flip $ frequencyControl $ \ freqs reso ->
      let resos = Flat.toSamples reso
          orderInt = NonNeg.toNumber order
      in  modifyModulated
             (Moog.parameter orderInt)
             (Moog.lowpassModifier orderInt)
             (Sig.zipWith FiltRec.Pole resos freqs)


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

{-# INLINE allpassFlangerPhase #-}
allpassFlangerPhase :: Trans.C a => a
allpassFlangerPhase = Allpass.flangerPhase


{- | Infinitely many equi-delayed exponentially decaying echos. -}
{-# INLINE comb #-}
comb :: (Hom.C sig, RealField.C t, Module.C y yv, Dim.C u, Storable yv) =>
   DN.T u t -> y -> Proc.T s u t (RP.T s sig yv -> RP.T s sig yv)
comb time gain =
   do t <- toTimeScalar time
      return $ Hom.processSamples
         (fromStorable . Comb.run (round t) gain . toStorable)


-- * auxiliary functions

{-# INLINE frequencyControl #-}
frequencyControl :: (Dim.C u, Field.C y) =>
      (Sig.T y -> t)
   -> Proc.T s u y (
        SigA.R s (Dim.Recip u) y y
     -> t)
frequencyControl f =
   do toFreq <- Proc.withParam toFrequencyScalar
      return $ \ freq -> f (SigA.scalarSamples toFreq freq)


{-# INLINE modifyModulated #-}
modifyModulated :: Hom.C sig =>
   (param -> ctrl) ->
   Modifier state ctrl y0 y1 ->
   Sig.T param ->
   RP.T s sig y0 ->
   RP.T s sig y1
modifyModulated makeParam modif params =
   Hom.processSamples (Sig.modifyModulated modif (Sig.map makeParam params))