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

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes

-}
module Synthesizer.State.Filter.Recursive.MovingAverage
(sumsStaticInt,
modulatedFrac,
) where

import qualified Synthesizer.State.Signal  as Sig
import qualified Synthesizer.State.Filter.Recursive.Integration as Integration

import qualified Synthesizer.State.Filter.Delay as Delay

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 PreludeBase
import NumericPrelude

{- |
Like 'Synthesizer.State.Filter.NonRecursive.sums' but in a recursive form.
This needs only linear time (independent of the window size)
but may accumulate rounding errors.

@
ys = xs * (1,0,0,0,-1) \/ (1,-1)
ys * (1,-1) = xs * (1,0,0,0,-1)
ys = xs * (1,0,0,0,-1) + ys * (0,1)
@
-}
{-# INLINE sumsStaticInt #-}
sumsStaticInt :: (Additive.C v) => Int -> Sig.T v -> Sig.T v
sumsStaticInt n xs =
Integration.run (xs - Delay.staticPos n xs)

{-
staticInt :: (Module.C a v, Additive.C v) => Int -> Sig.T v -> Sig.T v
staticInt n xs =
-}

{-
Sum of a part of a vector with negative sign for reverse order.
It adds from @from@ (inclusively) to @to@ (exclusively),
that is, it sums up @abs (to-from)@ values

sumFromTo :: (Additive.C v) => Int -> Int -> Sig.T v -> v
sumFromTo from to =
if from <= to
then          Sig.sum . Sig.take (to-from) . Sig.drop from
else negate . Sig.sum . Sig.take (from-to) . Sig.drop to
-}

{-# INLINE sumFromToFrac #-}
sumFromToFrac :: (RealField.C a, Module.C a v) => a -> a -> Sig.T v -> v
sumFromToFrac from to xs =
let (fromInt, fromFrac) = splitFraction from
(toInt,   toFrac)   = splitFraction to
in  case compare fromInt toInt of
EQ -> (to-from) *> Sig.index fromInt xs
LT ->
Sig.sum \$
Sig.zipWith id
(((1-fromFrac) *>) `Sig.cons`
Sig.replicate (toInt-fromInt-1) id `Sig.append`
Sig.singleton (toFrac *>)) \$
Sig.drop fromInt xs
GT ->
negate \$ Sig.sum \$
Sig.zipWith id
(((1-toFrac) *>) `Sig.cons`
Sig.replicate (fromInt-toInt-1) id `Sig.append`
Sig.singleton (fromFrac *>)) \$
Sig.drop toInt xs

{-
run \$

type Accumulator v a =
WriterT (Dual (Endo v)) (StateT (Sig.T v) Maybe a)

getNext :: Accumulator v a
getNext =
lift \$ StateT \$ viewListL

addAccum x = tell ((x+) \$!)

addNextWeighted :: Module.C a v => a -> Accumulator v ()
addAccum . (w *>) =<< getNext
-}

{-
newtype Accumulator v =
Accumulator ((v, Sig.T v) -> v -> (Sig.T v, v))

Accumulator \$ \(x,xs) s -> (xs, x+s)

addNextWeighted :: Module.C a v => a -> Accumulator v
Accumulator \$ \(x,xs) s -> (xs, a*>x + s)

bindAccum :: Accumulator v -> Accumulator v -> Accumulator v
bindAccum (Accumulator f) (Accumulator g) =
Accumulator \$ \x s0 ->
let (ys,s1) = f x s0
in  maybe s1 () (viewListL ys)
-}

{- |
Sig.T a must contain only non-negative elements.
-}
{-# INLINE sumDiffsModulated #-}
sumDiffsModulated :: (RealField.C a, Module.C a v) =>
a -> Sig.T a -> Sig.T v -> Sig.T v
sumDiffsModulated d ds =
Sig.init .
-- prevent negative d's since 'drop' cannot restore past values
Sig.zipWithTails (uncurry sumFromToFrac)
(Sig.zip (Sig.cons (d+1) ds) (Sig.map (1+) ds)) .
Sig.cons zero
{-
Sig.zipWithTails (uncurry sumFromToFrac)
(Sig.zip (Sig.cons d (Sig.map (subtract 1) ds)) ds)
-}

{-
sumsModulated :: (RealField.C a, Module.C a v) =>
Int -> Sig.T a -> Sig.T v -> Sig.T v
sumsModulated maxDInt ds xs =
let maxD  = fromIntegral maxDInt
posXs = sumDiffsModulated 0 ds xs
negXs = sumDiffsModulated maxD (Sig.map (maxD-) ds) (Delay.static maxDInt xs)
in  Integration.run (posXs - negXs)
-}

{- |
Shift sampling points by a half sample period
in order to preserve signals for window widths below 1.
-}
{-# INLINE sumsModulatedHalf #-}
sumsModulatedHalf :: (RealField.C a, Module.C a v) =>
Int -> Sig.T a -> Sig.T v -> Sig.T v
sumsModulatedHalf maxDInt ds xs =
let maxD  = fromIntegral maxDInt
d0    = maxD+0.5
delXs = Delay.staticPos maxDInt xs
posXs = sumDiffsModulated d0 (Sig.map (d0+) ds) delXs
negXs = sumDiffsModulated d0 (Sig.map (d0-) ds) delXs
in  Integration.run (posXs - negXs)

{-# INLINE modulatedFrac #-}
modulatedFrac :: (RealField.C a, Module.C a v) =>
Int -> Sig.T a -> Sig.T v -> Sig.T v
modulatedFrac maxDInt ds xs =
Sig.zipWith (\d y -> recip (2*d) *> y) ds \$
sumsModulatedHalf maxDInt ds xs

```