{-# LANGUAGE NoImplicitPrelude #-}
{- |
The distortion functions have slope 1 at zero,
if they are differentiable at that point, at all.
This ensures that signals with low amplitude
are only slightly altered.
Non-differentiable distortions try to have an overall slope of 1.
-}
module Synthesizer.Basic.Distortion (
   clip, logit,
   zigZag, sine,
   oddChebyshev, {- swing, -}
   quantize,
   powerSigned,
   ) where

import qualified Algebra.Transcendental        as Trans
import qualified Algebra.RealField             as RealField
import qualified Algebra.Field                 as Field
import qualified Algebra.RealRing              as RealRing
import qualified Algebra.Absolute              as Absolute
import qualified Algebra.Ring                  as Ring

import Data.List.HT (mapAdjacent, )
import Data.Ord.HT (limit, )

import NumericPrelude.Numeric
import NumericPrelude.Base


{- * Clipping -}

{- |
limit, fuzz booster
-}
clip :: (RealRing.C a) => a -> a
clip :: forall a. C a => a -> a
clip = forall a. Ord a => (a, a) -> a -> a
limit (forall a. C a => a -> a
negate forall a. C a => a
one, forall a. C a => a
one)

{- |
logit, tanh
-}
logit :: (Trans.C a) => a -> a
logit :: forall a. C a => a -> a
logit = forall a. C a => a -> a
tanh

{-
probit, error function
-}



{- * Wrapping -}

{- |
zig-zag
-}
zigZag :: (RealField.C a) => a -> a
zigZag :: forall a. C a => a -> a
zigZag a
x =
   let (Int
n,a
y) = forall a b. (C a, C b) => a -> (b, a)
splitFraction ((a
xforall a. C a => a -> a -> a
+a
1)forall a. C a => a -> a -> a
/a
2)
   in  if forall a. (C a, C a) => a -> Bool
even (Int
n::Int)
         then a
2forall a. C a => a -> a -> a
*a
y forall a. C a => a -> a -> a
- a
1
         else a
1 forall a. C a => a -> a -> a
- a
2forall a. C a => a -> a -> a
*a
y

{- |
sine
-}
sine :: (Trans.C a) => a -> a
sine :: forall a. C a => a -> a
sine = forall a. C a => a -> a
sin


{- |
Odd Chebyshev polynomial

@oddChebyshev n@ is an appropriately scaled Chebyshev polynomial of order @2*n+1@.
The argument @n@ must be non-negative.

> Graphics.Gnuplot.Simple.plotFuncs [Graphics.Gnuplot.Simple.YRange (-1,1)] (Graphics.Gnuplot.Simple.linearScale 1000 (-7,7::Double)) (List.map oddChebyshev [0..5])
-}
oddChebyshev :: (Trans.C a) => (Field.C a) => Int -> a -> a
oddChebyshev :: forall a. (C a, C a) => Int -> a -> a
oddChebyshev Int
n a
xn =
   let order :: Int
order = Int
2forall a. C a => a -> a -> a
*Int
nforall a. C a => a -> a -> a
+Int
1
       {-
       slope of normal Chebyshev polynomials at zero is @order@
       which can be seen when considering slope of @x -> cos (order * arccos x)@
       -}
       x :: a
x  = forall a. C a => Int -> a -> a
parityFlip Int
n (a
xn forall a. C a => a -> a -> a
/ forall a b. (C a, C b) => a -> b
fromIntegral Int
order)
       ys :: [a]
ys = a
1 forall a. a -> [a] -> [a]
: a
x forall a. a -> [a] -> [a]
: forall a b. (a -> a -> b) -> [a] -> [b]
mapAdjacent (\a
x0 a
x1 -> a
2forall a. C a => a -> a -> a
*a
xforall a. C a => a -> a -> a
*a
x1 forall a. C a => a -> a -> a
- a
x0) [a]
ys
   in  [a]
ys forall a. [a] -> Int -> a
!! Int
order

parityFlip :: Ring.C a => Int -> a -> a
parityFlip :: forall a. C a => Int -> a -> a
parityFlip Int
n a
x =
   if forall a. (C a, C a) => a -> Bool
even Int
n then a
x else -a
x


{- |
A polynomial function with zeros at every integral point
weighted in order to equalize the local extreme points.

However, the weighting is difficult enough,
that it might be easier to use just a truncated Taylor series of sine.

We could compute a weighting denominator polynomial
by dividing our equidistant zeros polynomial by the sine series.

equidist / weight = sine
weight = equidist / sine

However we have to normalize the zeros,
thus powers of pi enter the scene
and then power series division becomes inexact.
-}
_swing :: (Trans.C a) => (Field.C a) => Int -> a -> a
_swing :: forall a. (C a, C a) => Int -> a -> a
_swing Int
n a
x =
{-
   foldl (*) x
      (map
         (\ni ->
            let x2 = x^2
                n2 = ni^2
            in  (x2-n2)/sqrt(x2+n2))
         (take n (iterate (1+) 1)))
-}
   forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall a. C a => a -> a -> a
(*) a
x
      (forall a b. (a -> b) -> [a] -> [b]
map
         (\a
ni ->
            let x2 :: a
x2 = a
xforall a. C a => a -> Integer -> a
^Integer
2
                n2 :: a
n2 = a
niforall a. C a => a -> Integer -> a
^Integer
2
            in  (a
x2forall a. C a => a -> a -> a
-a
n2)forall a. C a => a -> a -> a
/(a
x2forall a. C a => a -> a -> a
+a
n2))
         (forall a. Int -> [a] -> [a]
take Int
n (forall a. (a -> a) -> a -> [a]
iterate (a
1forall a. C a => a -> a -> a
+) a
1)))
{-
   foldl (*) x
      (map (\ni -> (x/ni)^2-1)
         (take n (iterate (1+) 1)))
-}
{-
   let xu = iterate (1+) x
       xl = iterate (subtract 1) x
   in  foldl (*) x (take n (tail (zipWith (*) xu xl)))
-}
--   in  product (x : take n (tail xu) ++ take n (tail xl))



{- * Quantization -}

quantize :: (RealField.C a) => a -> a
quantize :: forall a. C a => a -> a
quantize a
x = forall a b. (C a, C b) => a -> b
fromIntegral (forall a b. (C a, C b) => a -> b
round a
x :: Int)


{- * other -}

{- |
Power function.
Roughly the map @\p x -> x**p@ but retains the sign of @x@.
-}
{-# INLINE powerSigned #-}
powerSigned :: (Absolute.C a, Trans.C a) => a -> a -> a
powerSigned :: forall a. (C a, C a) => a -> a -> a
powerSigned a
p a
x = forall a. C a => a -> a
signum a
x forall a. C a => a -> a -> a
* forall a. C a => a -> a
abs a
x forall a. C a => a -> a -> a
** a
p