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

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes
-}
module Synthesizer.Generic.Filter.NonRecursive where

import qualified Synthesizer.Generic.Signal as SigG
import qualified Synthesizer.Generic.Signal2 as SigG2

import qualified Synthesizer.Generic.Control as Ctrl
import qualified Synthesizer.State.Signal as SigS
import qualified Synthesizer.Plain.Filter.NonRecursive as Filt

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

import Algebra.Module( {- linearComb, -} (*>), )

import Data.Function.HT (nest, )
import Data.Tuple.HT (mapSnd, mapPair, )

import PreludeBase
import NumericPrelude as NP


{- * Envelope application -}

{-# INLINE negate #-}
negate ::
   (Additive.C a, SigG.Transform sig a) =>
   sig a -> sig a
negate = SigG.map Additive.negate

{-# INLINE amplify #-}
amplify ::
   (Ring.C a, SigG.Transform sig a) =>
   a -> sig a -> sig a
amplify v = SigG.map (v*)

{-# INLINE amplifyVector #-}
amplifyVector ::
   (Module.C a v, SigG.Transform sig v) =>
   a -> sig v -> sig v
amplifyVector v = SigG.map (v*>)

{-# INLINE envelope #-}
envelope ::
   (Ring.C a, SigG.Transform sig a) =>
      sig a  {-^ the envelope -}
   -> sig a  {-^ the signal to be enveloped -}
   -> sig a
envelope = SigG.zipWith (*)

{-# INLINE envelopeVector #-}
envelopeVector ::
   (Module.C a v, SigG.Read sig a, SigG.Transform sig v) =>
      sig a  {-^ the envelope -}
   -> sig v  {-^ the signal to be enveloped -}
   -> sig v
envelopeVector = SigG.zipWith (*>)



{-# INLINE fadeInOut #-}
fadeInOut ::
   (Field.C a, SigG.Write sig a) =>
   Int -> Int -> Int -> sig a -> sig a
fadeInOut tIn tHold tOut xs =
   let slopeIn  =                  recip (fromIntegral tIn)
       slopeOut = Additive.negate (recip (fromIntegral tOut))
       {-
       Since we use the size only for the internal envelope
       no laziness effect can be observed outside the function.
       We could also create the envelope as State.Signal.
       But I assume that concatenating chunks of an envelope
       is more efficient than concatenating generator loops.
       However, our intermediate envelope is still observable,
       because we have to use SigG.Write class.
       -}
       leadIn  = SigG.take tIn  $ Ctrl.linear SigG.defaultLazySize slopeIn  0
       leadOut = SigG.take tOut $ Ctrl.linear SigG.defaultLazySize slopeOut 1
       (partIn, partHoldOut) = SigG.splitAt tIn xs
       (partHold, partOut)   = SigG.splitAt tHold partHoldOut
   in  envelope leadIn partIn `SigG.append`
       partHold `SigG.append`
       envelope leadOut partOut


{- * Smoothing -}

{-# INLINE delay #-}
delay :: (Additive.C y, SigG.Write sig y) =>
   Int -> sig y -> sig y
delay =
   delayPad zero

{-# INLINE delayPad #-}
delayPad :: (SigG.Write sig y) =>
   y -> Int -> sig y -> sig y
delayPad z n =
   if n<0
     then SigG.drop (Additive.negate n)
     else SigG.append (SigG.replicate SigG.defaultLazySize n z)

{-# INLINE delayPos #-}
delayPos :: (Additive.C y, SigG.Write sig y) =>
   Int -> sig y -> sig y
delayPos n =
   SigG.append (SigG.replicate SigG.defaultLazySize n zero)

{-# INLINE delayNeg #-}
delayNeg :: (SigG.Transform sig y) =>
   Int -> sig y -> sig y
delayNeg = SigG.drop



{-# INLINE delayLazySize #-}
delayLazySize :: (Additive.C y, SigG.Write sig y) =>
   SigG.LazySize -> Int -> sig y -> sig y
delayLazySize size =
   delayPadLazySize size zero

{- |
The pad value @y@ must be defined,
otherwise the chunk size of the padding can be observed.
-}
{-# INLINE delayPadLazySize #-}
delayPadLazySize :: (SigG.Write sig y) =>
   SigG.LazySize -> y -> Int -> sig y -> sig y
delayPadLazySize size z n =
   if n<0
     then SigG.drop (Additive.negate n)
     else SigG.append (SigG.replicate size n z)

{-# INLINE delayPosLazySize #-}
delayPosLazySize :: (Additive.C y, SigG.Write sig y) =>
   SigG.LazySize -> Int -> sig y -> sig y
delayPosLazySize size n =
   SigG.append (SigG.replicate size n zero)



{-| Unmodulated non-recursive filter -}
{-# INLINE generic #-}
generic ::
   (Module.C a v, SigG.Transform sig a, SigG.Write sig v) =>
   sig a -> sig v -> sig v
generic m x =
   let mr = SigG.reverse m
       xp = delayPos (pred (SigG.length m)) x
   in  SigG.mapTails (SigG.linearComb mr) xp

{-
{- |
@eps@ is the threshold relatively to the maximum.
That is, if the gaussian falls below @eps * gaussian 0@,
then the function truncated.
-}
gaussian ::
   (Trans.C a, RealField.C a, Module.C a v) =>
   a -> a -> a -> sig v -> sig v
gaussian eps ratio freq =
   let var    = ratioFreqToVariance ratio freq
       area   = var * sqrt (2*pi)
       gau t  = exp (-(t/var)^2/2) / area
       width  = ceiling (var * sqrt (-2 * log eps))  -- inverse gau
       gauSmp = map (gau . fromIntegral) [-width .. width]
   in  drop width . generic gauSmp
-}

{-
GNUPlot.plotList [] (take 1000 $ gaussian 0.001 0.5 0.04 (Filter.Test.chirp 5000) :: [Double])

The filtered chirp must have amplitude 0.5 at 400 (0.04*10000).
-}

{-
  We want to approximate a Gaussian by a binomial filter.
  The latter one can be implemented by a convolutional power.
  However we still require a number of operations per sample
  which is proportional to the variance.
-}
{-# INLINE binomial #-}
binomial ::
   (Trans.C a, RealField.C a, Module.C a v, SigG.Transform sig v) =>
   a -> a -> sig v -> sig v
binomial ratio freq =
   let width = ceiling (2 * ratioFreqToVariance ratio freq ^ 2)
   in  SigG.drop width .
       nest (2*width) (amplifyVector (asTypeOf 0.5 freq) . binomial1)

{-
exp (-(t/var)^2/2) / area *> cis (2*pi*f*t)
  == exp (-(t/var)^2/2 +: 2*pi*f*t) / area
  == exp ((-t^2 +: 2*var^2*2*pi*f*t) / (2*var^2)) / area
  == exp ((t^2 - i*2*var^2*2*pi*f*t) / (-2*var^2)) / area
  == exp (((t^2 - i*var^2*2*pi*f)^2 + (var^2*2*pi*f)^2) / (-2*var^2)) / area
  == exp (((t^2 - i*var^2*2*pi*f)^2 / (-2*var^2) - (var*2*pi*f)^2/2)) / area

sumMap (\t -> exp (-(t/var)^2/2) / area *> cis (2*pi*f*t))
       [-infinity..infinity]
  ~ sumMap (\t -> exp (-(t/var)^2/2)) [-infinity..infinity]
       * exp (-(var*2*pi*f)^2/2) / area
  = exp (-(var*2*pi*f)^2/2)
-}
{- |
  Compute the variance of the Gaussian
  such that its Fourier transform has value @ratio@ at frequency @freq@.
-}
{-# INLINE ratioFreqToVariance #-}
ratioFreqToVariance :: (Trans.C a) => a -> a -> a
ratioFreqToVariance ratio freq =
   sqrt (Additive.negate (2 * log ratio)) / (2*pi*freq)
           -- inverse of the fourier transformed gaussian

{-# INLINE binomial1 #-}
binomial1 ::
   (Additive.C v, SigG.Transform sig v) => sig v -> sig v
binomial1 = SigG.mapAdjacent (+)





{- |
Moving (uniformly weighted) average in the most trivial form.
This is very slow and needs about @n * length x@ operations.
-}
{-# INLINE sums #-}
sums ::
   (Additive.C v, SigG.Transform sig v) =>
   Int -> sig v -> sig v
sums n = SigG.mapTails (SigG.sum . SigG.take n)


sumsDownsample2 ::
   (Additive.C v, SigG.Write sig v) =>
   SigG.LazySize -> sig v -> sig v
sumsDownsample2 cs =
   SigG.unfoldR cs (\xs ->
      flip fmap (SigG.viewL xs) $ \xxs0@(x0,xs0) ->
         SigG.switchL xxs0 {- xs0 is empty -}
            (\ x1 xs1 -> (x0+x1, xs1))
            xs0)

downsample2 ::
   (SigG.Write sig v) =>
   SigG.LazySize -> sig v -> sig v
downsample2 cs =
   SigG.unfoldR cs
      (fmap (mapSnd SigG.laxTail) . SigG.viewL)


{-
{- |
Given a list of numbers
and a list of sums of (2*k) of successive summands,
compute a list of the sums of (2*k+1) or (2*k+2) summands.

Eample for 2*k+1

@
 [0+1+2+3, 2+3+4+5, 4+5+6+7, ...] ->
    [0+1+2+3+4, 1+2+3+4+5, 2+3+4+5+6, 3+4+5+6+7, 4+5+6+7+8, ...]
@

Example for 2*k+2

@
 [0+1+2+3, 2+3+4+5, 4+5+6+7, ...] ->
    [0+1+2+3+4+5, 1+2+3+4+5+6, 2+3+4+5+6+7, 3+4+5+6+7+8, 4+5+6+7+8+9, ...]
@
-}
sumsUpsampleOdd :: (Additive.C v) => Int -> sig v -> sig v -> sig v
sumsUpsampleOdd n {- 2*k -} xs ss =
   let xs2k = drop n xs
   in  (head ss + head xs2k) :
          concat (zipWith3 (\s x0 x2k -> [x0+s, s+x2k])
                           (tail ss)
                           (downsample2 (tail xs))
                           (tail (downsample2 xs2k)))

sumsUpsampleEven :: (Additive.C v) => Int -> sig v -> sig v -> sig v
sumsUpsampleEven n {- 2*k -} xs ss =
   sumsUpsampleOdd (n+1) xs (zipWith (+) ss (downsample2 (drop n xs)))

sumsPyramid :: (Additive.C v) => Int -> sig v -> sig v
sumsPyramid n xs =
   let aux 1 ys = ys
       aux 2 ys = ys + tail ys
       aux m ys =
          let ysd = sumsDownsample2 ys
          in  if even m
                then sumsUpsampleEven (m-2) ys (aux (div (m-2) 2) ysd)
                else sumsUpsampleOdd  (m-1) ys (aux (div (m-1) 2) ysd)
   in  aux n xs


propSums :: Bool
propSums =
   let n  = 1000
       xs = [0::Double ..]
       naive   =              sums        n xs
       rec     = drop (n-1) $ sumsRec     n xs
       pyramid =              sumsPyramid n xs
   in  and $ take 1000 $
         zipWith3 (\x y z -> x==y && y==z) naive rec pyramid

-}

sumRange ::
   (Additive.C v, SigG.Transform sig v) =>
   sig v -> (Int,Int) -> v
sumRange =
   Filt.sumRangePrepare $ \ (l,r) ->
   SigG.sum . SigG.take (r-l) . SigG.drop l

pyramid ::
   (Additive.C v, SigG.Write sig v) =>
   Int -> sig v -> ([Int], [sig v])
pyramid height sig =
   let sizes =
          reverse $ take (1+height) $ iterate (2*) 1
   in  (sizes,
        scanl (flip sumsDownsample2) sig (map SigG.LazySize $ tail sizes))

sumRangeFromPyramid ::
   (Additive.C v, SigG.Write sig v) =>
   [sig v] -> (Int,Int) -> v
sumRangeFromPyramid =
   Filt.sumRangePrepare $ \(l0,r0) pyr0 ->
   case pyr0 of
      [] -> error "empty pyramid"
      (ps0:pss) ->
         foldr
            (\psNext k (l,r) ps s ->
               case r-l of
                  0 -> s
                  1 -> s + SigG.index ps l
                  _ ->
                     let (lh,ll) = NP.negate $ divMod (NP.negate l) 2
                         (rh,rl) = divMod r 2
                         {-# INLINE inc #-}
                         inc b x = if b==0 then id else (x+)
                     in  k (lh,rh) psNext $
                         inc ll (SigG.index ps l) $
                         inc rl (SigG.index ps (r-1)) $
                         s)
            (\(l,r) ps s ->
               s + (SigG.sum $ SigG.take (r-l) $ SigG.drop l ps))
            pss (l0,r0) ps0 zero

sumsPosModulated ::
   (Additive.C v, SigG2.Transform sig (Int,Int) v) =>
   sig (Int,Int) -> sig v -> sig v
sumsPosModulated ctrl xs =
   SigG2.zipWithTails (flip sumRange) ctrl xs


{- |
Moving average, where window bounds must be always non-negative.

The laziness granularity is @2^height@.
-}
sumsPosModulatedPyramid ::
   (Additive.C v, SigG.Transform sig (Int,Int), SigG.Write sig v) =>
   Int -> sig (Int,Int) -> sig v -> sig v
sumsPosModulatedPyramid height ctrl xs =
   let (sizes,pyr0) = pyramid height xs
       blockSize = head sizes
       pyrStarts =
          iterate (zipWith SigG.drop sizes) pyr0
       ctrlBlocks =
          SigS.toList $
          SigG.sliceVertical blockSize ctrl
   in  SigG.concat $
       zipWith
          (\pyr ->
              SigG.fromState (SigG.LazySize blockSize) .
              SigS.map (sumRangeFromPyramid pyr) .
              SigS.zipWith (\d -> mapPair ((d+), (d+))) (SigS.iterate (1+) 0) .
              SigG.toState)
          pyrStarts ctrlBlocks

{- |
The first argument is the amplification.
The main reason to introduce it,
was to have only a Module constraint instead of Field.
This way we can also filter stereo signals.
-}
movingAverageModulatedPyramid ::
   (Field.C a, Module.C a v,
    SigG2.Transform sig Int (Int,Int), SigG.Write sig v) =>
   a -> Int -> Int -> sig Int -> sig v -> sig v
movingAverageModulatedPyramid amp height maxC ctrl xs =
   SigG.zipWith (\c x -> (amp / fromIntegral (2*c+1)) *> x) ctrl $
   sumsPosModulatedPyramid height
      (SigG2.map (\c -> (maxC - c, maxC + c)) ctrl)
      (delay maxC xs)



{- * Filter operators from calculus -}

{- |
Forward difference quotient.
Shortens the signal by one.
Inverts 'Synthesizer.Generic.Filter.Recursive.Integration.run' in the sense that
@differentiate (zero : integrate x) == x@.
The signal is shifted by a half time unit.
-}
{-# INLINE differentiate #-}
differentiate ::
   (Additive.C v, SigG.Transform sig v) =>
   sig v -> sig v
differentiate x = SigG.mapAdjacent subtract x

{- |
Central difference quotient.
Shortens the signal by two elements,
and shifts the signal by one element.
(Which can be fixed by prepending an appropriate value.)
For linear functions this will yield
essentially the same result as 'differentiate'.
You obtain the result of 'differentiateCenter'
if you smooth the one of 'differentiate'
by averaging pairs of adjacent values.

ToDo: Vector variant
-}
{-
This implementation is a bit cumbersome,
but it fits both StorableVector and State.Signal
(since it avoids recomputation).
-}
{-# INLINE differentiateCenter #-}
differentiateCenter ::
   (Field.C v, SigG.Transform sig v) =>
   sig v -> sig v
differentiateCenter =
   SigG.drop 2 .
   SigG.crochetL
      (\x0 (x1,x2) -> Just ((x2-x0)/2, (x0,x1)))
      (zero,zero)

{- |
Second derivative.
It is @differentiate2 == differentiate . differentiate@
but 'differentiate2' should be faster.
-}
{-# INLINE differentiate2 #-}
differentiate2 ::
   (Additive.C v, SigG.Transform sig v) =>
   sig v -> sig v
differentiate2 = differentiate . differentiate