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

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes
-}
module Synthesizer.Plain.Filter.NonRecursive (
   amplify,
   amplifyVector,
   binomial,
   binomial1,
   delay,
   delayPad,
   differentiate,
   differentiate2,
   differentiateCenter,
   downsample2,
   envelope,
   envelopeVector,
   fadeInOut,
   fadeInOutAlt,
   gaussian,
   generic,
   genericAlt,
   minRange,
   movingAverageModulatedPyramid,
   sumRange,
   sumRangeFromPyramid,
   sums,
   sumsDownsample2,
   sumsPosModulated,
   sumsPosModulatedPyramid,
   sumsPyramid,

   -- for testing
   propGeneric,
   sumRangeFromPyramidFoldr,
   sumRangeFromPyramidRec,
   getRangeFromPyramid,
   pyramid,
   ) where

import Synthesizer.Basic.Filter.NonRecursive (
   unitSizesFromPyramid,
   sumRangePrepare,
   symmetricRangePrepare,
   ratioFreqToVariance,
   )

import qualified Synthesizer.Plain.Control as Ctrl
import qualified Synthesizer.Plain.Signal  as Sig

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 (mapPair, )
import Data.List.HT (sliceVertical, )
import Data.List (tails, )

import NumericPrelude.Numeric
import NumericPrelude.Base


{- * Envelope application -}

amplify :: (Ring.C a) => a -> Sig.T a -> Sig.T a
amplify :: forall a. C a => a -> T a -> T a
amplify a
v = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
va -> a -> a
forall a. C a => a -> a -> a
*)

amplifyVector :: (Module.C a v) => a -> Sig.T v -> Sig.T v
amplifyVector :: forall a v. C a v => a -> T v -> T v
amplifyVector = a -> T v -> T v
forall a v. C a v => a -> v -> v
(*>)


envelope :: (Ring.C a) =>
      Sig.T a  {-^ the envelope -}
   -> Sig.T a  {-^ the signal to be enveloped -}
   -> Sig.T a
envelope :: forall a. C a => T a -> T a -> T a
envelope = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. C a => a -> a -> a
(*)

envelopeVector :: (Module.C a v) =>
      Sig.T a  {-^ the envelope -}
   -> Sig.T v  {-^ the signal to be enveloped -}
   -> Sig.T v
envelopeVector :: forall a v. C a v => T a -> T v -> T v
envelopeVector = (a -> v -> v) -> [a] -> [v] -> [v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> v -> v
forall a v. C a v => a -> v -> v
(*>)



fadeInOut :: (Field.C a) => Int -> Int -> Int -> Sig.T a -> Sig.T a
fadeInOut :: forall a. C a => Int -> Int -> Int -> T a -> T a
fadeInOut Int
tIn Int
tHold Int
tOut T a
xs =
   let leadIn :: T a
leadIn  = Int -> T a -> T a
forall a. Int -> [a] -> [a]
take Int
tIn  (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$ a -> a -> T a
forall y. C y => y -> y -> T y
Ctrl.linearMultiscale (  a -> a
forall a. C a => a -> a
recip (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
tIn))  a
0
       leadOut :: T a
leadOut = Int -> T a -> T a
forall a. Int -> [a] -> [a]
take Int
tOut (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$ a -> a -> T a
forall y. C y => y -> y -> T y
Ctrl.linearMultiscale (- a -> a
forall a. C a => a -> a
recip (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
tOut)) a
1
       (T a
partIn, T a
partHoldOut) = Int -> T a -> (T a, T a)
forall a. Int -> [a] -> ([a], [a])
splitAt Int
tIn T a
xs
       (T a
partHold, T a
partOut) = Int -> T a -> (T a, T a)
forall a. Int -> [a] -> ([a], [a])
splitAt Int
tHold T a
partHoldOut
   in  T a -> T a -> T a
forall a. C a => T a -> T a -> T a
envelope T a
leadIn T a
partIn T a -> T a -> T a
forall a. [a] -> [a] -> [a]
++
       T a
partHold T a -> T a -> T a
forall a. [a] -> [a] -> [a]
++
       T a -> T a -> T a
forall a. C a => T a -> T a -> T a
envelope T a
leadOut T a
partOut


fadeInOutAlt :: (Field.C a) => Int -> Int -> Int -> Sig.T a -> Sig.T a
fadeInOutAlt :: forall a. C a => Int -> Int -> Int -> T a -> T a
fadeInOutAlt Int
tIn Int
tHold Int
tOut =
   ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a. a -> a
id
      (((Int -> a -> a) -> [Int] -> [a -> a]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
x a
y -> a
y a -> a -> a
forall a. C a => a -> a -> a
* Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
x a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
tIn) [Int
0..Int
tInInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1]) [a -> a] -> [a -> a] -> [a -> a]
forall a. [a] -> [a] -> [a]
++
       (Int -> (a -> a) -> [a -> a]
forall a. Int -> a -> [a]
replicate Int
tHold a -> a
forall a. a -> a
id) [a -> a] -> [a -> a] -> [a -> a]
forall a. [a] -> [a] -> [a]
++
       ((Int -> a -> a) -> [Int] -> [a -> a]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
x a
y -> a
y a -> a -> a
forall a. C a => a -> a -> a
* Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
x a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
tOut) [Int
tOutInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1,Int
tOutInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
2..Int
0]))



{- * Shift -}

delay :: Additive.C y => Int -> Sig.T y -> Sig.T y
delay :: forall y. C y => Int -> T y -> T y
delay = y -> Int -> T y -> T y
forall y. y -> Int -> T y -> T y
delayPad y
forall a. C a => a
zero

delayPad :: y -> Int -> Sig.T y -> Sig.T y
delayPad :: forall y. y -> Int -> T y -> T y
delayPad y
z Int
n =
   if Int
nInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
0
     then Int -> [y] -> [y]
forall a. Int -> [a] -> [a]
drop (Int -> Int
forall a. C a => a -> a
negate Int
n)
     else (Int -> y -> [y]
forall a. Int -> a -> [a]
replicate Int
n y
z [y] -> [y] -> [y]
forall a. [a] -> [a] -> [a]
++)


{- * Smoothing -}


{-| Unmodulated non-recursive filter -}
generic :: Module.C a v =>
   Sig.T a -> Sig.T v -> Sig.T v
generic :: forall a v. C a v => T a -> T v -> T v
generic T a
m T v
x =
   let mr :: T a
mr = T a -> T a
forall a. [a] -> [a]
reverse T a
m
       xp :: T v
xp = Int -> T v -> T v
forall y. C y => Int -> T y -> T y
delay (Int -> Int
forall a. Enum a => a -> a
pred (T a -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length T a
m)) T v
x
   in  (T v -> v) -> [T v] -> T v
forall a b. (a -> b) -> [a] -> [b]
map (T a -> T v -> v
forall a v. C a v => [a] -> [v] -> v
linearComb T a
mr) ([T v] -> [T v]
forall a. HasCallStack => [a] -> [a]
init (T v -> [T v]
forall a. [a] -> [[a]]
tails T v
xp))

{-|
Unmodulated non-recursive filter
Output has same length as the input.

It is elegant but leaks memory.
 -}
genericAlt :: Module.C a v =>
   Sig.T a -> Sig.T v -> Sig.T v
genericAlt :: forall a v. C a v => T a -> T v -> T v
genericAlt T a
m T v
x =
   (T v -> v) -> [T v] -> T v
forall a b. (a -> b) -> [a] -> [b]
map (T a -> T v -> v
forall a v. C a v => [a] -> [v] -> v
linearComb T a
m)
      ([T v] -> [T v]
forall a. HasCallStack => [a] -> [a]
tail ((T v -> v -> T v) -> T v -> T v -> [T v]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl ((v -> T v -> T v) -> T v -> v -> T v
forall a b c. (a -> b -> c) -> b -> a -> c
flip (:)) [] T v
x))


propGeneric :: (Module.C a v, Eq v) =>
   Sig.T a -> Sig.T v -> Bool
propGeneric :: forall a v. (C a v, Eq v) => T a -> T v -> Bool
propGeneric T a
m T v
x =
--   generic m x == genericAlt m x
   [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ([Bool] -> Bool) -> [Bool] -> Bool
forall a b. (a -> b) -> a -> b
$ (v -> v -> Bool) -> T v -> T v -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> v -> Bool
forall a. Eq a => a -> a -> Bool
(==) (T a -> T v -> T v
forall a v. C a v => T a -> T v -> T v
generic T a
m T v
x) (T a -> T v -> T v
forall a v. C a v => T a -> T v -> T v
genericAlt T a
m T v
x)



{- |
@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.T v -> Sig.T v
gaussian :: forall a v. (C a, C a, C a v) => a -> a -> a -> T v -> T v
gaussian a
eps a
ratio a
freq =
   let var :: a
var    = a -> a -> a
forall a. C a => a -> a -> a
ratioFreqToVariance a
ratio a
freq
       area :: a
area   = a
var a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
sqrt (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
forall a. C a => a
pi)
       gau :: a -> a
gau a
t  = a -> a
forall a. C a => a -> a
exp (-(a
ta -> a -> a
forall a. C a => a -> a -> a
/a
var)a -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2a -> a -> a
forall a. C a => a -> a -> a
/a
2) a -> a -> a
forall a. C a => a -> a -> a
/ a
area
       width :: Int
width  = a -> Int
forall b. C b => a -> b
forall a b. (C a, C b) => a -> b
ceiling (a
var a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
sqrt (-a
2 a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
log a
eps))  -- inverse gau
       gauSmp :: [a]
gauSmp = (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a
gau (a -> a) -> (Int -> a) -> Int -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral) [-Int
width .. Int
width]
   in  Int -> [v] -> [v]
forall a. Int -> [a] -> [a]
drop Int
width ([v] -> [v]) -> ([v] -> [v]) -> [v] -> [v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [v] -> [v]
forall a v. C a v => T a -> T v -> T v
generic [a]
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.
-}
binomial :: (Trans.C a, RealField.C a, Module.C a v) => a -> a -> Sig.T v -> Sig.T v
binomial :: forall a v. (C a, C a, C a v) => a -> a -> T v -> T v
binomial a
ratio a
freq =
   let width :: Int
width = a -> Int
forall b. C b => a -> b
forall a b. (C a, C b) => a -> b
ceiling (a
2 a -> a -> a
forall a. C a => a -> a -> a
* a -> a -> a
forall a. C a => a -> a -> a
ratioFreqToVariance a
ratio a
freq a -> Integer -> a
forall a. C a => a -> Integer -> a
^ Integer
2)
   in  Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
width (T v -> T v) -> (T v -> T v) -> T v -> T v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> (T v -> T v) -> T v -> T v
forall a. Int -> (a -> a) -> a -> a
nest (Int
2Int -> Int -> Int
forall a. C a => a -> a -> a
*Int
width) ((a -> a -> a
forall a. a -> a -> a
asTypeOf a
0.5 a
freq a -> T v -> T v
forall a v. C a v => a -> v -> v
*>) (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 v. C v => T v -> T v
binomial1)

binomial1 :: (Additive.C v) => Sig.T v -> Sig.T v
binomial1 :: forall v. C v => T v -> T v
binomial1 xt :: [v]
xt@(v
x:[v]
xs) = v
x v -> [v] -> [v]
forall a. a -> [a] -> [a]
: ([v]
xs [v] -> [v] -> [v]
forall a. C a => a -> a -> a
+ [v]
xt)
binomial1 [] = []






{- |
Moving (uniformly weighted) average in the most trivial form.
This is very slow and needs about @n * length x@ operations.
-}
sums :: (Additive.C v) => Int -> Sig.T v -> Sig.T v
sums :: forall y. C y => Int -> T y -> T y
sums Int
n = ([v] -> v) -> [[v]] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map ([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
n) ([[v]] -> [v]) -> ([v] -> [[v]]) -> [v] -> [v]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[v]] -> [[v]]
forall a. HasCallStack => [a] -> [a]
init ([[v]] -> [[v]]) -> ([v] -> [[v]]) -> [v] -> [[v]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [v] -> [[v]]
forall a. [a] -> [[a]]
tails



sumsDownsample2 :: (Additive.C v) => Sig.T v -> Sig.T v
sumsDownsample2 :: forall v. C v => T v -> T v
sumsDownsample2 (v
x0:v
x1:[v]
xs) = (v
x0v -> v -> v
forall a. C a => a -> a -> a
+v
x1) v -> [v] -> [v]
forall a. a -> [a] -> [a]
: [v] -> [v]
forall v. C v => T v -> T v
sumsDownsample2 [v]
xs
sumsDownsample2 [v]
xs         = [v]
xs

downsample2 :: Sig.T a -> Sig.T a
downsample2 :: forall a. [a] -> [a]
downsample2 (a
x0:a
_:[a]
xs) = a
x0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
forall a. [a] -> [a]
downsample2 [a]
xs
downsample2 [a]
xs        = [a]
xs


{- |
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.

Example 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.T v -> Sig.T v -> Sig.T v
sumsUpsampleOdd :: forall v. C v => Int -> T v -> T v -> T v
sumsUpsampleOdd Int
n {- 2*k -} T v
xs T v
ss =
   let xs2k :: T v
xs2k = Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
n T v
xs
   in  (T v -> v
forall a. HasCallStack => [a] -> a
head T v
ss v -> v -> v
forall a. C a => a -> a -> a
+ T v -> v
forall a. HasCallStack => [a] -> a
head T v
xs2k) v -> T v -> T v
forall a. a -> [a] -> [a]
:
          [T v] -> T v
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ((v -> v -> v -> T v) -> T v -> T v -> T v -> [T v]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\v
s v
x0 v
x2k -> [v
x0v -> v -> v
forall a. C a => a -> a -> a
+v
s, v
sv -> v -> v
forall a. C a => a -> a -> a
+v
x2k])
                           (T v -> T v
forall a. HasCallStack => [a] -> [a]
tail T v
ss)
                           (T v -> T v
forall a. [a] -> [a]
downsample2 (T v -> T v
forall a. HasCallStack => [a] -> [a]
tail T v
xs))
                           (T v -> T v
forall a. HasCallStack => [a] -> [a]
tail (T v -> T v
forall a. [a] -> [a]
downsample2 T v
xs2k)))

sumsUpsampleEven :: (Additive.C v) => Int -> Sig.T v -> Sig.T v -> Sig.T v
sumsUpsampleEven :: forall v. C v => Int -> T v -> T v -> T v
sumsUpsampleEven Int
n {- 2*k -} T v
xs T v
ss =
   Int -> T v -> T v -> T v
forall v. C v => Int -> T v -> T v -> T v
sumsUpsampleOdd (Int
nInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1) T v
xs ((v -> v -> v) -> T v -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> v -> v
forall a. C a => a -> a -> a
(+) T v
ss (T v -> T v
forall a. [a] -> [a]
downsample2 (Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
n T v
xs)))

sumsPyramid :: (Additive.C v) => Int -> Sig.T v -> Sig.T v
sumsPyramid :: forall y. C y => Int -> T y -> T y
sumsPyramid =
   let aux :: Int -> T a -> T a
aux Int
1 T a
ys = T a
ys
       aux Int
2 T a
ys = T a
ys T a -> T a -> T a
forall a. C a => a -> a -> a
+ T a -> T a
forall a. HasCallStack => [a] -> [a]
tail T a
ys  -- binomial1
       aux Int
m T a
ys =
          let ysd :: T a
ysd = T a -> T a
forall v. C v => T v -> T v
sumsDownsample2 T a
ys
          in  if Int -> Bool
forall a. (C a, C a) => a -> Bool
even Int
m
                then Int -> T a -> T a -> T a
forall v. C v => Int -> T v -> T v -> T v
sumsUpsampleEven (Int
mInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
2) T a
ys (Int -> T a -> T a
aux (Int -> Int -> Int
forall a. C a => a -> a -> a
div (Int
mInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
2) Int
2) T a
ysd)
                else Int -> T a -> T a -> T a
forall v. C v => Int -> T v -> T v -> T v
sumsUpsampleOdd  (Int
mInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1) T a
ys (Int -> T a -> T a
aux (Int -> Int -> Int
forall a. C a => a -> a -> a
div (Int
mInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1) Int
2) T a
ysd)
   in  Int -> T v -> T v
forall y. C y => Int -> T y -> T y
aux


{- |
Compute the sum of the values from index l to (r-1).
(I.e. somehow a right open interval.)
This can be used for implementation of a moving average filter.
However, its counterpart 'sumRangeFromPyramid'
is much faster for large windows.
-}
sumRange :: (Additive.C v) => Sig.T v -> (Int,Int) -> v
sumRange :: forall v. C v => T v -> (Int, Int) -> v
sumRange =
   ((Int, Int) -> T v -> v) -> T v -> (Int, Int) -> v
forall v source.
C v =>
((Int, Int) -> source -> v) -> source -> (Int, Int) -> v
sumRangePrepare (((Int, Int) -> T v -> v) -> T v -> (Int, Int) -> v)
-> ((Int, Int) -> T v -> v) -> T v -> (Int, Int) -> v
forall a b. (a -> b) -> a -> b
$ \ (Int
l,Int
r) ->
   T v -> v
forall a. C a => [a] -> a
sum (T v -> v) -> (T v -> T v) -> T v -> v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> T v -> T v
forall a. Int -> [a] -> [a]
take (Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l) (T v -> T v) -> (T v -> T v) -> T v -> T v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
l

pyramid :: (Additive.C v) => Sig.T v -> [Sig.T v]
pyramid :: forall v. C v => T v -> [T v]
pyramid = (T v -> T v) -> T v -> [T v]
forall a. (a -> a) -> a -> [a]
iterate T v -> T v
forall v. C v => T v -> T v
sumsDownsample2

{- |
This function should be much faster than 'sumRange'
but slower than the recursively implemented @movingAverage@.
However in contrast to @movingAverage@
it should not suffer from cancelation.
-}
sumRangeFromPyramid :: (Additive.C v) => [Sig.T v] -> (Int,Int) -> v
sumRangeFromPyramid :: forall v. C v => [T v] -> (Int, Int) -> v
sumRangeFromPyramid =
   ((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v
forall v source.
C v =>
((Int, Int) -> source -> v) -> source -> (Int, Int) -> v
sumRangePrepare (((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v)
-> ((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v
forall a b. (a -> b) -> a -> b
$ \(Int, Int)
lr0 [T v]
pyr0 ->
   T v -> v
forall a. C a => [a] -> a
sum (T v -> v) -> T v -> v
forall a b. (a -> b) -> a -> b
$ [T v] -> (Int, Int) -> T v
forall v. [T v] -> (Int, Int) -> T v
getRangeFromPyramid [T v]
pyr0 (Int, Int)
lr0

minRange :: (Ord v) => Sig.T v -> (Int,Int) -> v
minRange :: forall v. Ord v => T v -> (Int, Int) -> v
minRange =
   ((Int, Int) -> T v -> v) -> T v -> (Int, Int) -> v
forall source v.
((Int, Int) -> source -> v) -> source -> (Int, Int) -> v
symmetricRangePrepare (((Int, Int) -> T v -> v) -> T v -> (Int, Int) -> v)
-> ((Int, Int) -> T v -> v) -> T v -> (Int, Int) -> v
forall a b. (a -> b) -> a -> b
$ \ (Int
l,Int
r) ->
   T v -> v
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (T v -> v) -> (T v -> T v) -> T v -> v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> T v -> T v
forall a. Int -> [a] -> [a]
take (Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l) (T v -> T v) -> (T v -> T v) -> T v -> T v
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop Int
l

getRangeFromPyramid :: [Sig.T v] -> (Int,Int) -> [v]
getRangeFromPyramid :: forall v. [T v] -> (Int, Int) -> T v
getRangeFromPyramid [T v]
pyr0 (Int, Int)
lr0 =
   case [T v]
pyr0 of
      [] -> [Char] -> T v
forall a. HasCallStack => [Char] -> a
error [Char]
"empty pyramid"
      (T v
ps0:[T v]
pss) ->
         (T v -> (((Int, Int), T v) -> T v) -> ((Int, Int), T v) -> T v)
-> (((Int, Int), T v) -> T v) -> [T v] -> ((Int, Int), T v) -> T v
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
            (\T v
psNext ((Int, Int), T v) -> T v
k ((Int
l,Int
r), T v
ps) ->
               let (Int
lh,Int
ll) = - Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod (-Int
l) Int
2
                   (Int
rh,Int
rl) =   Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod   Int
r  Int
2
                   ls :: T v
ls = if Int
llInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 then [] else [T v
psT v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!!Int
l]
                   rs :: T v
rs = if Int
rlInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 then [] else [T v
psT v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!!(Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1)]
               in  case Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l of
                     Int
0 -> []
                     Int
1 -> [T v
psT v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!!Int
l]
                     Int
_ -> T v
ls T v -> T v -> T v
forall a. [a] -> [a] -> [a]
++ T v
rs T v -> T v -> T v
forall a. [a] -> [a] -> [a]
++ ((Int, Int), T v) -> T v
k ((Int
lh,Int
rh), T v
psNext))
            (\((Int
l,Int
r), T v
ps) -> Int -> T v -> T v
forall a. Int -> [a] -> [a]
take (Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l) (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
l T v
ps)
            [T v]
pss ((Int, Int)
lr0, T v
ps0)

{- mapAccumL cannot work since the pyramid might be infinitely high
   sum $ takeWhileJust $ snd $
   mapAccumL () lr0 $
   zip pyr0 $
   tail (Match.replicate pyr0 False ++ [True])
-}

sumRangeFromPyramidRec ::
   (Additive.C v) =>
   [Sig.T v] -> (Int,Int) -> v
sumRangeFromPyramidRec :: forall v. C v => [T v] -> (Int, Int) -> v
sumRangeFromPyramidRec =
   let recourse :: t -> (Int, Int) -> [[t]] -> t
recourse t
s (Int
l,Int
r) [[t]]
pyr =
          case [[t]]
pyr of
             ([t]
ps:[]) ->
                t
s t -> t -> t
forall a. C a => a -> a -> a
+ ([t] -> t
forall a. C a => [a] -> a
sum ([t] -> t) -> [t] -> t
forall a b. (a -> b) -> a -> b
$ Int -> [t] -> [t]
forall a. Int -> [a] -> [a]
take (Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l) ([t] -> [t]) -> [t] -> [t]
forall a b. (a -> b) -> a -> b
$ Int -> [t] -> [t]
forall a. Int -> [a] -> [a]
drop Int
l [t]
ps)
             ([t]
ps:[[t]]
pss) ->
                let (Int
lh,Int
ll) = - Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod (-Int
l) Int
2
                    (Int
rh,Int
rl) =   Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod   Int
r  Int
2
                    ls :: t
ls = if Int
llInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 then t
forall a. C a => a
zero else [t]
ps[t] -> Int -> t
forall a. HasCallStack => [a] -> Int -> a
!!Int
l
                    rs :: t
rs = if Int
rlInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 then t
forall a. C a => a
zero else [t]
ps[t] -> Int -> t
forall a. HasCallStack => [a] -> Int -> a
!!(Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1)
                in  case Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l of
                      Int
0 -> t
s
                      Int
1 -> t
st -> t -> t
forall a. C a => a -> a -> a
+[t]
ps[t] -> Int -> t
forall a. HasCallStack => [a] -> Int -> a
!!Int
l
                      Int
_ -> t -> (Int, Int) -> [[t]] -> t
recourse (t
s t -> t -> t
forall a. C a => a -> a -> a
+ t
ls t -> t -> t
forall a. C a => a -> a -> a
+ t
rs) (Int
lh,Int
rh) [[t]]
pss
             [] -> [Char] -> t
forall a. HasCallStack => [Char] -> a
error [Char]
"empty pyramid"
   in  ((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v
forall v source.
C v =>
((Int, Int) -> source -> v) -> source -> (Int, Int) -> v
sumRangePrepare (v -> (Int, Int) -> [T v] -> v
forall {t}. C t => t -> (Int, Int) -> [[t]] -> t
recourse v
forall a. C a => a
zero)

sumRangeFromPyramidFoldr ::
   (Additive.C v) =>
   [Sig.T v] -> (Int,Int) -> v
sumRangeFromPyramidFoldr :: forall v. C v => [T v] -> (Int, Int) -> v
sumRangeFromPyramidFoldr =
   ((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v
forall v source.
C v =>
((Int, Int) -> source -> v) -> source -> (Int, Int) -> v
sumRangePrepare (((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v)
-> ((Int, Int) -> [T v] -> v) -> [T v] -> (Int, Int) -> v
forall a b. (a -> b) -> a -> b
$ \(Int, Int)
lr0 [T v]
pyr0 ->
   case [T v]
pyr0 of
      [] -> [Char] -> v
forall a. HasCallStack => [Char] -> a
error [Char]
"empty pyramid"
      (T v
ps0:[T v]
pss) ->
         (T v
 -> ((Int, Int) -> T v -> v -> v) -> (Int, Int) -> T v -> v -> v)
-> ((Int, Int) -> T v -> v -> v)
-> [T v]
-> (Int, Int)
-> T v
-> v
-> v
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
            (\T v
psNext (Int, Int) -> T v -> v -> v
k (Int
l,Int
r) T v
ps v
s ->
               case Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l of
                  Int
0 -> v
s
                  Int
1 -> v
s v -> v -> v
forall a. C a => a -> a -> a
+ T v
psT v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!!Int
l
                  Int
_ ->
                     let (Int
lh,Int
ll) = - Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod (-Int
l) Int
2
                         (Int
rh,Int
rl) =   Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod   Int
r  Int
2
                         {-# INLINE inc #-}
                         inc :: a -> a -> a -> a
inc a
b a
x = if a
ba -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0 then a -> a
forall a. a -> a
id else (a
xa -> a -> a
forall a. C a => a -> a -> a
+)
                     in  (Int, Int) -> T v -> v -> v
k (Int
lh,Int
rh) T v
psNext (v -> v) -> v -> v
forall a b. (a -> b) -> a -> b
$
                         Int -> v -> v -> v
forall {a} {a}. (Eq a, C a, C a) => a -> a -> a -> a
inc Int
ll (T v
psT v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!!Int
l) (v -> v) -> v -> v
forall a b. (a -> b) -> a -> b
$
                         Int -> v -> v -> v
forall {a} {a}. (Eq a, C a, C a) => a -> a -> a -> a
inc Int
rl (T v
psT v -> Int -> v
forall a. HasCallStack => [a] -> Int -> a
!!(Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1)) (v -> v) -> v -> v
forall a b. (a -> b) -> a -> b
$
                         v
s)
            (\(Int
l,Int
r) T v
ps v
s ->
               v
s v -> v -> v
forall a. C a => a -> a -> a
+ (T v -> v
forall a. C a => [a] -> a
sum (T v -> v) -> T v -> v
forall a b. (a -> b) -> a -> b
$ Int -> T v -> T v
forall a. Int -> [a] -> [a]
take (Int
rInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
l) (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
l T v
ps))
            [T v]
pss (Int, Int)
lr0 T v
ps0 v
forall a. C a => a
zero

sumsPosModulated :: (Additive.C v) =>
   Sig.T (Int,Int) -> Sig.T v -> Sig.T v
sumsPosModulated :: forall v. C v => T (Int, Int) -> T v -> T v
sumsPosModulated T (Int, Int)
ctrl T v
xs =
   (T v -> (Int, Int) -> v) -> [T v] -> T (Int, Int) -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith T v -> (Int, Int) -> v
forall v. C v => T v -> (Int, Int) -> v
sumRange (T v -> [T v]
forall a. [a] -> [[a]]
tails T v
xs) T (Int, Int)
ctrl

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

The laziness granularity is @2^height@.
-}
sumsPosModulatedPyramid :: (Additive.C v) =>
   Int -> Sig.T (Int,Int) -> Sig.T v -> Sig.T v
sumsPosModulatedPyramid :: forall v. C v => Int -> T (Int, Int) -> T v -> T v
sumsPosModulatedPyramid Int
height T (Int, Int)
ctrl T v
xs =
   let blockSize :: Int
blockSize = Int
2 Int -> Integer -> Int
forall a. C a => a -> Integer -> a
^ Int -> Integer
forall a b. (C a, C b) => a -> b
fromIntegral Int
height
       pyr0 :: [T v]
pyr0 = Int -> [T v] -> [T v]
forall a. Int -> [a] -> [a]
take (Int
1Int -> Int -> Int
forall a. C a => a -> a -> a
+Int
height) ([T v] -> [T v]) -> [T v] -> [T v]
forall a b. (a -> b) -> a -> b
$ T v -> [T v]
forall v. C v => T v -> [T v]
pyramid T v
xs
       sizes :: [Int]
sizes = [T v] -> [Int]
forall signal. [signal] -> [Int]
unitSizesFromPyramid [T v]
pyr0
       pyrStarts :: [[T v]]
pyrStarts =
          ([T v] -> [T v]) -> [T v] -> [[T v]]
forall a. (a -> a) -> a -> [a]
iterate ((Int -> T v -> T v) -> [Int] -> [T v] -> [T v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> T v -> T v
forall a. Int -> [a] -> [a]
drop [Int]
sizes) [T v]
pyr0
       ctrlBlocks :: [T (Int, Int)]
ctrlBlocks =
          (T (Int, Int) -> T (Int, Int)) -> [T (Int, Int)] -> [T (Int, Int)]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> (Int, Int) -> (Int, Int))
-> [Int] -> T (Int, Int) -> T (Int, Int)
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
d -> (Int -> Int, Int -> Int) -> (Int, Int) -> (Int, Int)
forall a c b d. (a -> c, b -> d) -> (a, b) -> (c, d)
mapPair ((Int
dInt -> Int -> Int
forall a. C a => a -> a -> a
+), (Int
dInt -> Int -> Int
forall a. C a => a -> a -> a
+))) ([Int] -> T (Int, Int) -> T (Int, Int))
-> [Int] -> T (Int, Int) -> T (Int, Int)
forall a b. (a -> b) -> a -> b
$ (Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate (Int
1Int -> Int -> Int
forall a. C a => a -> a -> a
+) Int
0) ([T (Int, Int)] -> [T (Int, Int)])
-> [T (Int, Int)] -> [T (Int, Int)]
forall a b. (a -> b) -> a -> b
$
          Int -> T (Int, Int) -> [T (Int, Int)]
forall a. Int -> [a] -> [[a]]
sliceVertical Int
blockSize T (Int, Int)
ctrl
   in  [T v] -> T v
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([T v] -> T v) -> [T v] -> T v
forall a b. (a -> b) -> a -> b
$
       ([T v] -> T (Int, Int) -> T v)
-> [[T v]] -> [T (Int, Int)] -> [T v]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
          (\[T v]
pyr -> ((Int, Int) -> v) -> T (Int, Int) -> T v
forall a b. (a -> b) -> [a] -> [b]
map ([T v] -> (Int, Int) -> v
forall v. C v => [T v] -> (Int, Int) -> v
sumRangeFromPyramid [T v]
pyr))
          [[T v]]
pyrStarts [T (Int, Int)]
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.

A control value @n@ corresponds to filter window size @2*n+1@.
-}
movingAverageModulatedPyramid ::
   (Field.C a, Module.C a v) =>
   a -> Int -> Int -> Sig.T Int -> Sig.T v -> Sig.T v
movingAverageModulatedPyramid :: forall a v. (C a, C a v) => a -> Int -> Int -> [Int] -> T v -> T v
movingAverageModulatedPyramid a
amp Int
height Int
maxC [Int]
ctrl T v
xs =
   (Int -> v -> v) -> [Int] -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
c v
x -> (a
amp a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (Int
2Int -> Int -> Int
forall a. C a => a -> a -> a
*Int
cInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
1)) a -> v -> v
forall a v. C a v => a -> v -> v
*> v
x) [Int]
ctrl (T v -> T v) -> T v -> T v
forall a b. (a -> b) -> a -> b
$
   Int -> T (Int, Int) -> T v -> T v
forall v. C v => Int -> T (Int, Int) -> T v -> T v
sumsPosModulatedPyramid Int
height
      ((Int -> (Int, Int)) -> [Int] -> T (Int, Int)
forall a b. (a -> b) -> [a] -> [b]
map (\Int
c -> (Int
maxC Int -> Int -> Int
forall a. C a => a -> a -> a
- Int
c, Int
maxC Int -> Int -> Int
forall a. C a => a -> a -> a
+ Int
c Int -> Int -> Int
forall a. C a => a -> a -> a
+ Int
1)) [Int]
ctrl)
      (Int -> T v -> T v
forall y. C y => Int -> T y -> T y
delay Int
maxC T v
xs)



{- * Filter operators from calculus -}

{- |
Forward difference quotient.
Shortens the signal by one.
Inverts 'Synthesizer.Plain.Filter.Recursive.Integration.run' in the sense that
@differentiate (zero : integrate x) == x@.
The signal is shifted by a half time unit.
-}
differentiate :: Additive.C v => Sig.T v -> Sig.T v
differentiate :: forall v. C v => T v -> T v
differentiate T v
x = (v -> v -> v) -> T v -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> v -> v
forall a. C a => a -> a -> a
subtract T v
x (T v -> T v
forall a. HasCallStack => [a] -> [a]
tail T v
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
-}
differentiateCenter :: Field.C v => Sig.T v -> Sig.T v
differentiateCenter :: forall v. C v => T v -> T v
differentiateCenter T v
x =
   (v -> v) -> T v -> T v
forall a b. (a -> b) -> [a] -> [b]
map ((v
1v -> v -> v
forall a. C a => a -> a -> a
/v
2)v -> v -> v
forall a. C a => a -> a -> a
*) (T v -> T v) -> T v -> T v
forall a b. (a -> b) -> a -> b
$
   (v -> v -> v) -> T v -> T v -> T v
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith v -> v -> v
forall a. C a => a -> a -> a
subtract T v
x (T v -> T v
forall a. HasCallStack => [a] -> [a]
tail (T v -> T v
forall a. HasCallStack => [a] -> [a]
tail T v
x))

{- |
Second derivative.
It is @differentiate2 == differentiate . differentiate@
but 'differentiate2' should be faster.
-}
differentiate2 :: Additive.C v => Sig.T v -> Sig.T v
differentiate2 :: forall v. C v => T v -> T v
differentiate2 T v
xs0 =
   let xs1 :: T v
xs1 = T v -> T v
forall a. HasCallStack => [a] -> [a]
tail T v
xs0
       xs2 :: T v
xs2 = T v -> T v
forall a. HasCallStack => [a] -> [a]
tail T v
xs1
   in  (v -> v -> v -> v) -> T v -> T v -> T v -> T v
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\v
x0 v
x1 v
x2 -> v
x0v -> v -> v
forall a. C a => a -> a -> a
+v
x2v -> v -> v
forall a. C a => a -> a -> a
-(v
x1v -> v -> v
forall a. C a => a -> a -> a
+v
x1)) T v
xs0 T v
xs1 T v
xs2