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

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes
-}
module Synthesizer.Plain.Filter.Recursive.MovingAverage (
   sumsStaticInt,
   modulatedFrac,
   ) where

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

import Synthesizer.Plain.Filter.NonRecursive (delay, )

import qualified Algebra.Module                as Module
import qualified Algebra.RealField             as RealField
import qualified Algebra.Additive              as Additive

import Control.Monad.Fix (fix)
import Data.List (tails)

import NumericPrelude.Numeric
import NumericPrelude.Base



{- |
Like 'Synthesizer.Plain.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)
@
-}
sumsStaticInt :: (Additive.C v) => Int -> Sig.T v -> Sig.T v
sumsStaticInt :: forall v. C v => Int -> T v -> T v
sumsStaticInt Int
n T v
xs =
   (T v -> T v) -> T v
forall a. (a -> a) -> a
fix (\T v
ys -> let (T v
xs0,T v
xs1) = Int -> T v -> (T v, T v)
forall a. Int -> [a] -> ([a], [a])
splitAt Int
n T v
xs
               in  (T v
xs0 T v -> T v -> T v
forall a. [a] -> [a] -> [a]
++ (T v
xs1T v -> T v -> T v
forall a. C a => a -> a -> a
-T v
xs)) T v -> T v -> T v
forall a. C a => a -> a -> a
+ (v
forall a. C a => a
zerov -> T v -> T v
forall a. a -> [a] -> [a]
:T v
ys))

{-
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 :: forall v. C v => Int -> Int -> T v -> v
_sumFromTo Int
from Int
to =
   if Int
from Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
to
     then          [v] -> v
forall a. C a => [a] -> a
sum ([v] -> v) -> ([v] -> [v]) -> [v] -> v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [v] -> [v]
forall a. Int -> [a] -> [a]
take (Int
toInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
from) ([v] -> [v]) -> ([v] -> [v]) -> [v] -> [v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [v] -> [v]
forall a. Int -> [a] -> [a]
drop Int
from
     else v -> v
forall a. C a => a -> a
negate (v -> v) -> ([v] -> v) -> [v] -> v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [v] -> v
forall a. C a => [a] -> a
sum ([v] -> v) -> ([v] -> [v]) -> [v] -> v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [v] -> [v]
forall a. Int -> [a] -> [a]
take (Int
fromInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
to) ([v] -> [v]) -> ([v] -> [v]) -> [v] -> [v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [v] -> [v]
forall a. Int -> [a] -> [a]
drop Int
to

{-
It would be a nice approach to interpolate not just linearly at the borders
but in a way that the cut-off frequency is perfectly suppressed.
However suppression depends on the phase shift of the wave.
Actually, we could use a complex factor, but does this help?
-}
sumFromToFrac :: (RealField.C a, Module.C a v) => a -> a -> Sig.T v -> v
sumFromToFrac :: forall a v. (C a, C a v) => a -> a -> T v -> v
sumFromToFrac a
from a
to T v
xs =
   let (Int
fromInt, a
fromFrac) = a -> (Int, a)
forall b. C b => a -> (b, a)
forall a b. (C a, C b) => a -> (b, a)
splitFraction a
from
       (Int
toInt,   a
toFrac)   = a -> (Int, a)
forall b. C b => a -> (b, a)
forall a b. (C a, C b) => a -> (b, a)
splitFraction a
to
   in  case Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Int
fromInt Int
toInt of
          Ordering
EQ -> (a
toa -> a -> a
forall a. C a => a -> a -> a
-a
from) a -> v -> v
forall a v. C a v => a -> v -> v
*> (T v
xs T v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!! Int
fromInt)
          Ordering
LT ->
            T v -> v
forall a. C a => [a] -> a
sum (T v -> v) -> T v -> v
forall a b. (a -> b) -> a -> b
$
            ((v -> v) -> v -> v) -> [v -> v] -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (v -> v) -> v -> v
forall a. a -> a
id
               (((a
1a -> a -> a
forall a. C a => a -> a -> a
-a
fromFrac) a -> v -> v
forall a v. C a v => a -> v -> v
*>) (v -> v) -> [v -> v] -> [v -> v]
forall a. a -> [a] -> [a]
:
                Int -> (v -> v) -> [v -> v]
forall a. Int -> a -> [a]
replicate (Int
toIntInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
fromIntInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1) v -> v
forall a. a -> a
id [v -> v] -> [v -> v] -> [v -> v]
forall a. [a] -> [a] -> [a]
++
                (a
toFrac a -> v -> v
forall a v. C a v => a -> v -> v
*>) (v -> v) -> [v -> v] -> [v -> v]
forall a. a -> [a] -> [a]
:
                []) (T v -> T v) -> T v -> T v
forall a b. (a -> b) -> a -> b
$
            Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
fromInt T v
xs
          Ordering
GT ->
            v -> v
forall a. C a => a -> a
negate (v -> v) -> v -> v
forall a b. (a -> b) -> a -> b
$ T v -> v
forall a. C a => [a] -> a
sum (T v -> v) -> T v -> v
forall a b. (a -> b) -> a -> b
$
            ((v -> v) -> v -> v) -> [v -> v] -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (v -> v) -> v -> v
forall a. a -> a
id
               (((a
1a -> a -> a
forall a. C a => a -> a -> a
-a
toFrac) a -> v -> v
forall a v. C a v => a -> v -> v
*>) (v -> v) -> [v -> v] -> [v -> v]
forall a. a -> [a] -> [a]
:
                Int -> (v -> v) -> [v -> v]
forall a. Int -> a -> [a]
replicate (Int
fromIntInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
toIntInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1) v -> v
forall a. a -> a
id [v -> v] -> [v -> v] -> [v -> v]
forall a. [a] -> [a] -> [a]
++
                (a
fromFrac a -> v -> v
forall a v. C a v => a -> v -> v
*>) (v -> v) -> [v -> v] -> [v -> v]
forall a. a -> [a] -> [a]
:
                []) (T v -> T v) -> T v -> T v
forall a b. (a -> b) -> a -> b
$
            Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
toInt T v
xs



{- |
Sig.T a must contain only non-negative elements.
-}
sumDiffsModulated :: (RealField.C a, Module.C a v) =>
   a -> Sig.T a -> Sig.T v -> Sig.T v
sumDiffsModulated :: forall a v. (C a, C a v) => a -> T a -> T v -> T v
sumDiffsModulated a
d T a
ds =
   -- prevent negative d's since 'drop' cannot restore past values
   (a -> a -> T v -> v) -> T a -> T a -> [T v] -> T v
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 a -> a -> T v -> v
forall a v. (C a, C a v) => a -> a -> T v -> v
sumFromToFrac ((a
da -> a -> a
forall a. C a => a -> a -> a
+a
1) a -> T a -> T a
forall a. a -> [a] -> [a]
: T a
ds) ((a -> a) -> T a -> T a
forall a b. (a -> b) -> [a] -> [b]
map (a
1a -> a -> a
forall a. C a => a -> a -> a
+) T a
ds) ([T v] -> T v) -> (T v -> [T v]) -> T v -> T v
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   [T v] -> [T v]
forall a. HasCallStack => [a] -> [a]
init ([T v] -> [T v]) -> (T v -> [T v]) -> T v -> [T v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [T v] -> [T v]
forall a. HasCallStack => [a] -> [a]
init ([T v] -> [T v]) -> (T v -> [T v]) -> T v -> [T v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T v -> [T v]
forall a. [a] -> [[a]]
tails (T v -> [T v]) -> (T v -> T v) -> T v -> [T v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (v
forall a. C a => a
zerov -> T v -> T v
forall a. a -> [a] -> [a]
:)
{-
   zipWith3 sumFromToFrac (d : map (subtract 1) ds) ds .
   init . tails
-}

_sumsModulated :: (RealField.C a, Module.C a v) =>
   Int -> Sig.T a -> Sig.T v -> Sig.T v
_sumsModulated :: forall a v. (C a, C a v) => Int -> T a -> T v -> T v
_sumsModulated Int
maxDInt T a
ds T v
xs =
   let maxD :: a
maxD  = Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
maxDInt
       posXs :: T v
posXs = a -> T a -> T v -> T v
forall a v. (C a, C a v) => a -> T a -> T v -> T v
sumDiffsModulated a
0 T a
ds T v
xs
       negXs :: T v
negXs = a -> T a -> T v -> T v
forall a v. (C a, C a v) => a -> T a -> T v -> T v
sumDiffsModulated a
maxD ((a -> a) -> T a -> T a
forall a b. (a -> b) -> [a] -> [b]
map (a
maxDa -> a -> a
forall a. C a => a -> a -> a
-) T a
ds) (Int -> T v -> T v
forall v. C v => Int -> T v -> T v
delay Int
maxDInt T v
xs)
   in  T v -> T v
forall v. C v => T v -> T v
Integration.run (T v
posXs T v -> T v -> T v
forall a. C a => a -> a -> a
- T v
negXs)

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

{-
*Synthesizer.Plain.Filter.NonRecursive> movingAverageModulated 10 (replicate 10 (3::Double) ++ [1.1,2.2,2.6,0.7,0.1,0.1]) (repeat (1::Double))
[0.5,0.6666666666666666,0.8333333333333333,1.0,1.0,1.0,1.0,1.0,1.0,1.0,0.9999999999999999,1.0,0.9999999999999998,0.999999999999999,0.9999999999999942,0.9999999999999942]
-}

modulatedFrac :: (RealField.C a, Module.C a v) =>
   Int -> Sig.T a -> Sig.T v -> Sig.T v
modulatedFrac :: forall a v. (C a, C a v) => Int -> T a -> T v -> T v
modulatedFrac Int
maxDInt T a
ds T v
xs =
   (a -> v -> v) -> T a -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
d v
y -> a -> a
forall a. C a => a -> a
recip (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
d) a -> v -> v
forall a v. C a v => a -> v -> v
*> v
y) T a
ds (T v -> T v) -> T v -> T v
forall a b. (a -> b) -> a -> b
$
   Int -> T a -> T v -> T v
forall a v. (C a, C a v) => Int -> T a -> T v -> T v
sumsModulatedHalf Int
maxDInt T a
ds T v
xs