```{-# OPTIONS -O2 -fno-implicit-prelude -fglasgow-exts #-}
{- |
Copyright   :  (c) Henning Thielemann 2006

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

Basic waveforms

If you want to use parametrized waves with two parameters
then zip your parameter signals and apply 'uncurry' to the wave function.
-}
module Synthesizer.Basic.Wave where

import qualified Synthesizer.Plain.ToneModulation as ToneMod
import qualified Synthesizer.Plain.Interpolation  as Interpolation
import Data.Array ((!), listArray)

import qualified Synthesizer.Basic.Phase as Phase

import qualified Algebra.RealTranscendental    as RealTrans
import qualified Algebra.Transcendental        as Trans
import qualified Algebra.RealField             as RealField
import qualified Algebra.Algebraic             as Algebraic
import qualified Algebra.Module                as Module
import qualified Algebra.Field                 as Field
import qualified Algebra.Real                  as Real
import qualified Algebra.Ring                  as Ring

import qualified MathObj.Polynomial as Poly
import qualified Number.Complex     as Complex

import Synthesizer.Utility (swap)

import NumericPrelude.Condition (select, )
import NumericPrelude

-- import qualified Prelude as P
import PreludeBase

{- * Definition and construction -}

newtype T t y = Cons {decons :: Phase.T t -> y}

{-# INLINE fromFunction #-}
fromFunction :: (t -> y) -> (T t y)
fromFunction wave = Cons (wave . Phase.toRepresentative)

{- * Operations on waves -}

{-# INLINE raise #-}
raise :: (Additive.C y) => y -> T t y -> T t y
raise y = distort (y+)

{-# INLINE amplify #-}
amplify :: (Ring.C y) => y -> T t y -> T t y
amplify k = distort (k*)

{-# INLINE distort #-}
distort :: (y -> z) -> T t y -> T t z
distort g (Cons f) = Cons (g . f)

{-# INLINE apply #-}
apply :: T t y -> (Phase.T t -> y)
apply = decons

{-# INLINE zero #-}
{-# INLINE (+) #-}
{-# INLINE (-) #-}
{-# INLINE negate #-}
zero = Cons (const zero)
(+) (Cons f) (Cons g) = Cons (\t -> f t + g t)
(-) (Cons f) (Cons g) = Cons (\t -> f t - g t)
negate = distort negate

instance Module.C a y => Module.C a (T t y) where
{-# INLINE (*>) #-}
s *> w = distort (s*>) w

{- |
Turn an unparametrized waveform into a parametrized one,
where the parameter is a phase offset.
This way you express a phase modulated oscillator
using a shape modulated oscillator.
-}
{-# SPECULATE phaseOffset :: (T Double b) -> (Double -> T Double b) #-}
{-# INLINE phaseOffset #-}
phaseOffset :: (RealField.C a) => T a b -> (a -> T a b)
phaseOffset (Cons wave) offset =
Cons (wave . Phase.increment offset)

{- * Examples -}

{- ** unparameterized -}

{- | map a phase to value of a sine wave -}
{-# SPECULATE sine :: Double -> Double #-}
{-# INLINE sine #-}
sine :: Trans.C a => T a a
sine = fromFunction \$ \x -> sin (2*pi*x)

{-# INLINE cosine #-}
cosine :: Trans.C a => T a a
cosine = fromFunction \$ \x -> cos (2*pi*x)

{-# INLINE helix #-}
helix :: Trans.C a => T a (Complex.T a)
helix = fromFunction \$ \x -> Complex.cis (2*pi*x)

{- |
Approximation of sine by parabolas.
Surprisingly not really faster than 'sine'.
-}
{-# INLINE fastSine2 #-}
fastSine2 :: (Ord a, Ring.C a) => T a a
fastSine2 = fromFunction \$ \x ->
if 2*x<1
then 1 - sqr (4*x-1)
else sqr (4*x-3) - 1

{- |
Approximation of sine by fourth order polynomials.
-}
{-# INLINE fastSine4 #-}
fastSine4 :: (Ord a, Trans.C a) => T a a
fastSine4 = fromFunction \$ \x ->
-- minimal least squares fit
let pi2 = pi*pi
pi3 = pi2*pi
c = 3*((10080/pi2 - 1050) / pi3 + 1) -- 0.2248391014
{-# INLINE bow #-}
bow y = let y2 = y*y in 1-y2*(1+c*(1-y2))
in  if 2*x<1
then   bow (4*x-1)
else - bow (4*x-3)
{-
add a residue to fastSine2 and choose 'c' which minimizes the squared error
in  if 2*x<1
then let y = (4*x-1)^2 in 1-y-c*y*(1-y)
else let y = (4*x-3)^2 in y-1+c*y*(1-y)
-}

{-
GNUPlot.plotFuncs [] (GNUPlot.linearScale 1000 (0,1::Double)) [sine, fastSine2, fastSine4]
-}

{- | saw tooth,
it's a ramp down in order to have a positive coefficient for the first partial sine
-}
{-# SPECULATE saw :: Double -> Double #-}
{-# INLINE saw #-}
saw :: Ring.C a => T a a
saw = fromFunction \$ \x -> 1-2*x

{- |
This wave has the same absolute Fourier coefficients as 'saw'
but the partial waves are shifted by 90 degree.
That is, it is the Hilbert transform of the saw wave.
The formula is derived from 'sawComplex'.
-}
{-# INLINE sawCos #-}
sawCos :: (Real.C a, Trans.C a) => T a a
sawCos = fromFunction \$ \x -> log (2 * sin (pi*x)) * (-2/pi)

{- |
@sawCos + i*saw@

This is an analytic function and thus it may be used for frequency shifting.

The formula can be derived from the power series of the logarithm function.
-}
{-# INLINE sawComplex #-}
sawComplex ::
(Complex.Power a, RealTrans.C a) =>
T a (Complex.T a)
sawComplex = fromFunction \$ \x -> log (1 + Complex.cis (-pi*(1-2*x))) * (-2/pi)
{-
GNUPlot.plotFuncs [] (GNUPlot.linearScale 100 (0,1::Double)) [Complex.real . sawComplex, sawCos]

GNUPlot.plotFuncs [] (GNUPlot.linearScale 100 (0,1::Double)) [sawCos, composedHarmonics (take 20 \$ harmonic 0 0 : map (\n -> harmonic 0.25 ((2/pi) / fromInteger n)) [1..])]
-}

{-
Matching implementation that do not match 'saw' exactly.

sawCos :: (Real.C a, Trans.C a) => T a a
sawCos = fromFunction \$ \x -> log (2 * abs (cos (pi*x)))

sawComplex ::
(Complex.Power a, Trans.C a) =>
T a (Complex.T a)
sawComplex = fromFunction \$ \x -> log (1 + Complex.cis (2*pi*x))
-}

{- | square -}
{-# SPECULATE square :: Double -> Double #-}
{-# INLINE square #-}
square :: (Ord a, Ring.C a) => T a a
square = fromFunction \$ \x -> if 2*x<1 then 1 else -1

{- |
This wave has the same absolute Fourier coefficients as 'square'
but the partial waves are shifted by 90 degree.
That is, it is the Hilbert transform of the saw wave.
-}
{-# INLINE squareCos #-}
squareCos :: (RealField.C a, Trans.C a) => T a a
squareCos = fromFunction \$ \x ->
log (abs (tan (pi*x))) * (-2/pi)
-- sawCos x - sawCos (fraction (0.5-x))

{- |
@squareCos + i*square@

This is an analytic function and thus it may be used for frequency shifting.

The formula can be derived from the power series of the area tangens function.
-}
{-# INLINE squareComplex #-}
squareComplex ::
(Complex.Power a, RealTrans.C a) =>
T a (Complex.T a)
squareComplex = fromFunction \$ \x ->
{- these formulas are equivalent but wrong

log (0 +: 2 * sine x) * (2/pi)

log ((1 - Complex.cis (-2*pi*x)) *
(1 + Complex.cis ( 2*pi*x))) * (2/pi)

sawComplex x + sawComplex (0.5-x)
-}

{-
The Fourier series is equal to the power series of 'atanh'.
-}
atanh (Complex.cis (2*pi*x)) * (4/pi)
{-
GNUPlot.plotFuncs [] (GNUPlot.linearScale 100 (0,1::Double)) [squareCos, composedHarmonics (take 20 \$ zipWith (\b n -> harmonic 0.25 (if b then (4/pi) / fromInteger n else 0)) (cycle [False,True]) [0..])]
-}

{- | triangle -}
{-# SPECULATE triangle :: Double -> Double #-}
{-# INLINE triangle #-}
triangle :: (Ord a, Ring.C a) => T a a
triangle = fromFunction \$ \x ->
let x4 = 4*x
in  select (2-x4)
[(x4<1, x4),
(x4>3, x4-4)]

{-

int(arctan(x)/x,x);

- polylog(2, x*I)*1/2*I + polylog(2, x*(-I))*1/2*I

series(int(arctan(x)/x,x),x,10);

x - 1/9*x^3 + 1/25*x^5 - 1/49*x^7 + 1/81*x^9 + O(x^11)

int(arctan(I*x)/(I*x),x);
int(arctanh(x)/(x),x);

1/2*polylog(2, x) - 1/2*polylog(2, -x)
int(1/x*arctanh(x), x)

polylog(2,x) = dilog(1-x);    -- dilog is implemented in GSL for complex arguments
polylog(2,x) = hypergeom([1,1,1],[2,2],x) * x;

series(int(arctan(I*x)/(I*x),x),x,10);

x + 1/9*x^3 + 1/25*x^5 + 1/49*x^7 + 1/81*x^9 + O(x^11)
-}

sample :: (RealField.C a) =>
Interpolation.T a v -> [v] -> T a v
sample ip wave =
let len = length wave
arr = listArray (0, pred len) wave
in  fromFunction \$ \ phase ->
let (n,q) = splitFraction (phase * fromIntegral len)
xs = map (arr!) (map (flip mod len)
(enumFrom (n - Interpolation.offset ip)))
--                map (arr!) (enumFromTo (n - Interpolation.offset ip)) ++ cycle wave
in  Interpolation.func ip q xs

{- ** discretely parameterized -}

{- |
A truncated cosine. This has rich overtones.
-}
truncOddCosine :: Trans.C a =>
Int -> T a a
truncOddCosine k =
let f = pi * fromIntegral (2*k+1)
in  fromFunction \$ \ x -> cos (f*x)

{- |
For parameter zero this is 'saw'.
-}
truncOddTriangle :: (RealField.C a) =>
Int -> T a a
truncOddTriangle k =
let s = fromIntegral (2*k+1)
in  fromFunction \$ \ x ->
let (n,frac) = splitFraction (s*x)
in  if even (n::Int)
then 1-2*frac
else 2*frac-1

{- ** continuously parameterized -}

{- |
A truncated cosine plus a ramp that guarantees a bump of high 2 at the boundaries.

It is @truncCosine (2 * fromIntegral n + 0.5) == truncOddCosine (2*n)@
-}
truncCosine :: Trans.C a =>
a -> T a a
truncCosine k =
let f = 2 * pi * k
s = 2 * (sin (f*0.5) - 1)
in  fromFunction \$ \ x0 ->
let x = x0-0.5
in  - sin (f*x) + s*x
{-
GNUPlot.plotFuncs [] (GNUPlot.linearScale 1000 (0,1::Double)) (map truncCosine [0.5,0.7..2.5])
-}

truncTriangle :: (RealField.C a) =>
a -> T a a
truncTriangle k =
let tr x =
let (n,frac) = splitFraction (2*k*x+0.5)
in  if even (n::Int)
then 1-2*frac
else 2*frac-1
s = 2 * (1 + tr 0.5)
in  fromFunction \$ \ x0 ->
let x = x0-0.5
in  tr x - s*x
{-
GNUPlot.plotFuncs [] (GNUPlot.linearScale 1000 (0,1::Double)) (map truncTriangle [0,0.25..2.5])
-}

{- |
Power function.
-}

{- |
Roughly the map @\x p -> x**p@
but retains the sign of @x@ and
normalizes the mapping over @[-1,1]@ to L2 norm of 1.
-}
{-# INLINE powerNormed #-}
powerNormed :: (Real.C a, Trans.C a) => a -> T a a
powerNormed p = fromFunction \$ \x -> power01Normed p (2*x-1)

-- | auxiliary
{-# INLINE power01Normed #-}
power01Normed :: (Real.C a, Trans.C a) => a -> a -> a
power01Normed p x = (p+0.5) * powerSigned p x

-- | auxiliary
{-# INLINE powerSigned #-}
powerSigned :: (Real.C a, Trans.C a) => a -> a -> a
powerSigned p x = signum x * abs x ** p

{- |
Tangens hyperbolicus allows interpolation
between some kind of saw tooth and square wave.
In principle it is not necessary
because you can distort a saw tooth oscillation by @map tanh@.
-}
logitSaw :: (Trans.C a) => a -> T a a
logitSaw c = distort tanh \$ amplify c saw

{- |
Tangens hyperbolicus of a sine allows interpolation
between some kind of sine and square wave.
In principle it is not necessary
because you can distort a square oscillation by @map tanh@.
-}
logitSine :: (Trans.C a) => a -> T a a
logitSine c = distort tanh \$ amplify c sine

{- |
Interpolation between 'sine' and 'square'.
-}
{-# INLINE sineSquare #-}
sineSquare :: (Real.C a, Trans.C a) =>
a {- ^ 0 for 'sine', 1 for 'square' -}
-> T a a
sineSquare c =
distort (powerSigned (1-c)) sine

{- |
Interpolation between 'fastSine2' and 'saw'.
We just shrink the parabola towards the borders
and insert a linear curve such that its slope matches the one of the parabola.
-}
{-# INLINE piecewiseParabolaSaw #-}
piecewiseParabolaSaw :: (Algebraic.C a, Ord a) =>
a {- ^ 0 for 'fastSine2', 1 for 'saw' -}
-> T a a
piecewiseParabolaSaw c =
let xb  = (1 - sqrt c) / 2
y x = 1 - ((4*x - (1-c))/(1-c))^2
in  fromFunction \$ \ x ->
select
((2*x - 1)/(2*xb - 1) * y xb)
[(x <   xb,   y x),
(x > 1-xb, - y (1-x))]

{-
equ0 c x =
let y  = 1 - ((4*x - (3+c))/(1-c))^2
secant  = y/(x-1/2)
tangent = - 8 * (4*x - (3+c))/(1-c)^2
in  (tangent, secant)

equ1 c x =
let secant  = (1 - ((4*x - (3+c))/(1-c))^2)/(x-1/2)
tangent = - 8 * (4*x - (3+c))/(1-c)^2
in  (tangent, secant)

equ2 c x =
(1, ((4*x - (3+c))/(1-c))^2
- 8 * (x-1/2) * (4*x - (3+c))/(1-c)^2)

equ3 c x =
((1-c)^2,
(4*x - (3+c) - 4 * (2*x-1)) * (4*x - (3+c)))

equ4 c x =
(4*x - (1-c)) * (4*x - (3+c)) + (1-c)^2

equ5 c x =
(4*x - 2) ^ 2 - (1+c)^2 + (1-c)^2

equ6 c x =
(4*x - 2) ^ 2 - 4*c
-}

{- |
Interpolation between 'sine' and 'saw'.
We just shrink the sine towards the borders
and insert a linear curve such that its slope matches the one of the sine.
-}
{-# INLINE piecewiseSineSaw #-}
piecewiseSineSaw :: (Trans.C a, Ord a) =>
a {- ^ 0 for 'sine', 1 for 'saw' -}
-> T a a
piecewiseSineSaw c =
let {- This simple fix point iteration converges very slow for small 'c',
maybe we should use a Newton iteration. -}
iter z = iterate (\zi -> pi + atan (zi - pi / (1-c))) z !! 10
xb = (1-c)/(2*pi) * iter 0
-- iter (xInit * (2*pi) / (1-c))
-- xb  = (1 - sqrt c) / 2
-- y x = sine (x/(1-c))
y x = sin (2*pi*x/(1-c))
in  fromFunction \$ \ x -> select
((2*x - 1)/(2*xb - 1) * y xb)
[(x <   xb,   y x),
(x > 1-xb, - y (1-x))]

{-
equ0 c x =
let secant  = 2 * sin (2*pi*x/(1-c)) / (2*x - 1)
tangent = 2*pi/(1-c) * cos (2*pi*x/(1-c))
in  (tangent, secant)

iter0 c x =
-- secant / tangent
-- (x - 1/2) = tan (2*pi*x/(1-c)) * (1-c) / (2*pi)
tan (2*pi*x/(1-c)) * (1-c) / (2*pi) + 1/2

iter1 c x =
(1-c)/(2*pi) * (pi + atan ((x - 1/2) * (2*pi) / (1-c)))

iter2 c x =
let iter z = iterate (\zi -> pi + atan (zi - pi / (1-c))) z !! 10
in  (1-c)/(2*pi) * iter (x * (2*pi) / (1-c))
-}

{- |
Interpolation between 'sine' and 'saw'
with smooth intermediate shapes but no perfect saw.
-}
{-# INLINE sineSawSmooth #-}
sineSawSmooth :: (Trans.C a) =>
a {- ^ 0 for 'sine', 1 for 'saw' -}
-> T a a
sineSawSmooth c =
distort (\x -> sin (affineComb c (pi * x, asin x * 2))) saw

{- |
Interpolation between 'sine' and 'saw'
with perfect saw, but sharp intermediate shapes.
-}
{-# INLINE sineSawSharp #-}
sineSawSharp :: (Trans.C a) =>
a {- ^ 0 for 'sine', 1 for 'saw' -}
-> T a a
sineSawSharp c =
distort (\x -> sin (affineComb c (pi * x, asin x))) saw

affineComb :: Ring.C a => a -> (a,a) -> a
affineComb phase (x0,x1) = (1-phase)*x0 + phase*x1

{-
{- |
Smooth saw generated by a quintic polynomial function.
Unfortunately if 'c' approaches the right border,
the function will overshoot the 'y' range (-1,1).
-}
quinticSaw :: Field.C a =>
a  {- ^ position of the right minimum -}
-> a
-> a
quinticSaw c x =
let (s,t) = ToneMod.solveSLE2 ((c^2-1, 3*c^2-1), (c^4-1, 5*c^4-1)) (-1/c,0)
r = - s - t
x2 = x^2
in  x * (r + x2 * (s + x2*t))
{-
r*x + s*  x^3 + t*  x^5
0 = r   + s       + t
-1 = r*c + s*  c^3 + t*  c^5
0 = r   + s*3*c^2 + t*5*c^4

-1/c = r   + s*  c^2 + t*  c^4

-1/c = s*(c^2-1)   + t*(c^4-1)
0 = s*(3*c^2-1) + t*(5*c^4-1)
-}
-}

{- |
saw with space
-}
{-# SPECULATE sawPike :: Double -> Double -> Double #-}
{-# INLINE sawPike #-}
sawPike :: (Ord a, Field.C a) =>
a {- ^ pike width ranging from 0 to 1, 1 yields 'saw' -}
-> T a a
sawPike r = fromFunction \$ \x ->
if x<r
then 1-2/r*x
else 0

{- |
triangle with space
-}
{-# SPECULATE trianglePike :: Double -> Double -> Double #-}
{-# INLINE trianglePike #-}
trianglePike :: (Real.C a, Field.C a) =>
a  {- ^ pike width ranging from 0 to 1, 1 yields 'triangle' -}
-> T a a
trianglePike r = fromFunction \$ \x ->
if x < 1/2
then max 0 (1 - abs (4*x-1) / r)
else min 0 (abs (4*x-3) / r - 1)

{- |
triangle with space and shift
-}
{-# SPECULATE trianglePikeShift :: Double -> Double -> Double -> Double #-}
{-# INLINE trianglePikeShift #-}
trianglePikeShift :: (Real.C a, Field.C a) =>
a  {- ^ pike width ranging from 0 to 1 -}
-> a  {- ^ shift ranges from -1 to 1; 0 yields 'trianglePike' -}
-> T a a
trianglePikeShift r s = fromFunction \$ \x ->
if x < 1/2
then max 0 (1 - abs (4*x-1+s*(r-1)) / r)
else min 0 (abs (4*x-3+s*(1-r)) / r - 1)

{- |
square with space,
can also be generated by mixing square waves with different phases
-}
{-# SPECULATE squarePike :: Double -> Double -> Double #-}
{-# INLINE squarePike #-}
squarePike :: (Real.C a) =>
a  {- ^ pike width ranging from 0 to 1, 1 yields 'square' -}
-> T a a
squarePike r = fromFunction \$ \x ->
if 2*x < 1
then if abs(4*x-1)<r then  1 else 0
else if abs(4*x-3)<r then -1 else 0

{- |
square with space and shift
-}
{-# SPECULATE squarePikeShift :: Double -> Double -> Double -> Double #-}
{-# INLINE squarePikeShift #-}
squarePikeShift :: (Real.C a) =>
a  {- ^ pike width ranging from 0 to 1 -}
-> a  {- ^ shift ranges from -1 to 1; 0 yields 'squarePike' -}
-> T a a
squarePikeShift r s = fromFunction \$ \x ->
if 2*x < 1
then if abs(4*x-1+s*(r-1))<r then  1 else 0
else if abs(4*x-3+s*(1-r))<r then -1 else 0

{- |
square with different times for high and low
-}
{-# SPECULATE squareAsymmetric :: Double -> Double -> Double #-}
{-# INLINE squareAsymmetric #-}
squareAsymmetric :: (Ord a, Ring.C a) =>
a  {- ^ value between -1 and 1 controlling the ratio of high and low time:
-1 turns the high time to zero,
1 makes the low time zero,
0 yields 'square' -}
-> T a a
squareAsymmetric r = fromFunction \$ \x ->
if 2*x < r+1 then 1 else -1

{- | Like 'squareAsymmetric' but with zero average.
It could be simulated by adding two saw oscillations
with 180 degree phase difference and opposite sign.
-}
{-# SPECULATE squareBalanced :: Double -> Double -> Double #-}
{-# INLINE squareBalanced #-}
squareBalanced :: (Ord a, Ring.C a) => a -> T a a
squareBalanced r =
raise (-r) \$ squareAsymmetric r

{- |
triangle
-}
{-# SPECULATE sawPike :: Double -> Double -> Double #-}
{-# INLINE triangleAsymmetric #-}
triangleAsymmetric :: (Ord a, Field.C a) =>
a  {- ^ asymmetry parameter ranging from -1 to 1:
For 0 you obtain the usual triangle.
For -1 you obtain a falling saw tooth starting with its maximum.
For 1 you obtain a rising saw tooth starting with a zero. -}
-> T a a
triangleAsymmetric r = fromFunction \$ \x ->
select ((2-4*x)/(1-r))
[(4*x < 1+r, 4/(1+r)*x),
(4*x > 3-r, 4/(1+r)*(x-1))]

{- |
Mixing 'trapezoid' and 'trianglePike' you can get back a triangle wave form
-}
{-# SPECULATE trapezoid :: Double -> Double -> Double #-}
{-# INLINE trapezoid #-}
trapezoid :: (Real.C a, Field.C a) =>
a  {- ^ width of the plateau ranging from 0 to 1:
0 yields 'triangle', 1 yields 'square' -}
-> T a a
trapezoid w = fromFunction \$ \x ->
if x < 1/2
then min   1  ((1 - abs (4*x-1)) / (1-w))
else max (-1) ((abs (4*x-3) - 1) / (1-w))

{- |
trapezoid with distinct high and low time
-}
{-# SPECULATE trapezoidAsymmetric :: Double -> Double -> Double -> Double #-}
{-# INLINE trapezoidAsymmetric #-}
trapezoidAsymmetric :: (Real.C a, Field.C a) =>
a  {- ^ sum of the plateau widths ranging from 0 to 1:
0 yields 'triangleAsymmetric',
1 yields 'squareAsymmetric' -}
-> a  {- ^ asymmetry of the plateau widths ranging from -1 to 1 -}
-> T a a
trapezoidAsymmetric w r = fromFunction \$ \x ->
let c0 = 1+w*r
c1 = 1-w*r
in  if 2*x < c0
then min   1  ((c0 - abs (4*x-c0)) / (1-w))
else max (-1) ((abs (4*(1-x)-c1) - c1) / (1-w))
{-
let c = w*r+1
in  if 2*x < c
then min   1  ((1 - abs (4*x/c-1))*c/(1-w))
else max (-1) ((abs (4*(1-x)/(2-c)-1) - 1)*(2-c)/(1-w))
-}
{-
let c = (w*r+1)/2
in  if x < c
then min   1  ((1 - abs (2*x/c-1))*2*c/(1-w))
else max (-1) ((abs (2*(1-x)/(1-c)-1) - 1)*2*(1-c)/(1-w))
-}

{- |
trapezoid with distinct high and low time
-}
{-# SPECULATE trapezoidBalanced :: Double -> Double -> Double -> Double #-}
{-# INLINE trapezoidBalanced #-}
trapezoidBalanced :: (Real.C a, Field.C a) => a -> a -> T a a
trapezoidBalanced w r =
raise (-w*r) \$ trapezoidAsymmetric w r

{- |
We assume that a tone was generated by a shape modulated oscillator.
We try to reconstruct the wave function
(with parameters shape control and phase)
from a tone by interpolation.

The unit for the shape control parameter is the sampling period.
That is the shape parameter is a time parameter
pointing to a momentary shape of the prototype signal.
Of course this momentary shape does not exist
and we can only guess it using interpolation.

At the boundaries we repeat the outermost shapes
that can be reconstructed entirely from interpolated data
(that is, no extrapolation is needed).
This way we cannot reproduce the shape at the boundaries
because we have no data for cyclically extending it.
On the other hand this method guarantees a nice wave shape
with the required fractional period.

It must be
@length tone >=
Interpolation.number ipStep +
Interpolation.number ipLeap * ceiling period@.
-}
sampledTone :: (RealField.C a) =>
Interpolation.T a v ->
Interpolation.T a v ->
a -> [v] -> a -> T a v
sampledTone ipLeap ipStep period tone shape = Cons \$ \phase ->
uncurry (ToneMod.interpolateCell ipLeap ipStep) \$
ToneMod.sampledToneCell
(ToneMod.makePrototype ipLeap ipStep period tone)
shape phase

{- |
Interpolate first within waves and then across waves,
which is simpler but maybe less efficient.
-}
sampledToneAlt :: (RealField.C a) =>
Interpolation.T a v ->
Interpolation.T a v ->
a -> [v] -> a -> T a v
sampledToneAlt ipLeap ipStep period tone shape = Cons \$ \phase ->
uncurry (ToneMod.interpolateCell ipStep ipLeap . swap) \$
ToneMod.sampledToneAltCell
(ToneMod.makePrototype ipLeap ipStep period tone)
shape phase

{-
*Synthesizer.Basic.Wave>
GNUPlot.plotFunc [] (GNUPlot.linearScale 1000 (0,12)) (\t -> sampledTone Interpolation.linear Interpolation.linear (6::Double) ([-5,-3,-1,1,3,5,-4,-4,-4,4,4,4]++replicate 20 0) t (t/6))

*Synthesizer.Plain.Oscillator>
let period = 6.3::Double in GNUPlot.plotFunc [] (GNUPlot.linearScale 1000 (-10,20)) (\t -> Wave.sampledTone Interpolation.linear Interpolation.cubic period (take 20 \$ staticSine 0 (1/period)) t (t/period))
-}

{- |
This is similar to Polar coordinates,
but the range of the phase is from @0@ to @1@, @0@ to @2*pi@.
-}
data Harmonic a =
Harmonic {harmonicPhase :: Phase.T a, harmonicAmplitude :: a}

{-# INLINE harmonic #-}
harmonic :: Phase.T a -> a -> Harmonic a
harmonic = Harmonic

{- |
Specify the wave by its harmonics.

The function is implemented quite efficiently
by applying the Horner scheme to a polynomial with complex coefficients
(the harmonic parameters)
using a complex exponential as argument.
-}
{-# INLINE composedHarmonics #-}
composedHarmonics :: Trans.C a => [Harmonic a] -> T a a
composedHarmonics hs =
let p = Poly.fromCoeffs \$
map (\h -> Complex.fromPolar (harmonicAmplitude h)
(2*pi * Phase.toRepresentative (harmonicPhase h))) hs
in  distort (Complex.imag . Poly.evaluate p) helix
{-
GNUPlot.plotFunc [] (GNUPlot.linearScale 1000 (0,1::Double)) (composedHarmonics [harmonic 0 0, harmonic 0 0, harmonic 0 0, harmonic 0.25 1])
-}
```