//################################### hoa.lib ############################################
// Faust library for high order ambisonic. Its official prefix is `ho`.
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/hoa.lib>
//########################################################################################

/************************************************************************
 ************************************************************************
FAUST library file
Copyright (C) 2003-2012 GRAME, Centre National de Creation Musicale
----------------------------------------------------------------------
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2.1 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA.

EXCEPTION TO THE LGPL LICENSE : As a special exception, you may create a
larger FAUST program which directly or indirectly imports this library
file and still distribute the compiled code generated by the FAUST
compiler, or a modified version of this compiled code, under your own
copyright and license. This EXCEPTION TO THE LGPL LICENSE explicitly
grants you the right to freely choose the license for the resulting
compiled code. In particular the resulting compiled code has no obligation
to be LGPL or GPL. For example you are free to choose a commercial or
closed source license or any other license if you decide so.

 ************************************************************************
 ************************************************************************/

ma = library("maths.lib");
si = library("signals.lib");
ba = library("basics.lib");
os = library("oscillators.lib");
ho = library("hoa.lib");
ro = library("routes.lib");
de = library("delays.lib");

declare name "High Order Ambisonics library";
declare version "1.4.0";
declare author "Pierre Guillot";
declare author "Eliott Paris";
declare author "Julien Colafrancesco";
declare author "Wargreen";
declare author "Alain Bonardi";
declare author "Paul Goutmann";
declare copyright "2012-2013 Guillot, Paris, Colafrancesco, CICM labex art H2H, U. Paris 8, 2019 Wargreen, 2022 Bonardi, Goutmann";

//============================Encoding/decoding Functions=================================
//========================================================================================

//----------------------`(ho.)encoder`---------------------------------
// Ambisonic encoder. Encodes a signal in the circular harmonics domain
// depending on an order of decomposition and an angle.
//
// #### Usage
//
// ```
// encoder(N, x, a) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `x`: the signal
// * `a`: the angle
//----------------------------------------------------------------
encoder(0, x, a) = x;
encoder(N, x, a) = encoder(N-1, x, a), x*sin(N*a), x*cos(N*a);


//-------`(ho.)rEncoder`----------
// Ambisonic encoder in 2D including source rotation. A mono signal is encoded at a certain ambisonic order
// with two possible modes: either rotation with an angular speed, or static with a fixed angle (when speed is zero).
//
// #### Usage
//
// ```
// _ : rEncoder(N, sp, a, it) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `sp`: the azimuth speed expressed as angular speed (2PI/sec), positive or negative
// * `a`: the fixed azimuth when the rotation stops (sp = 0) in radians
// * `it` : interpolation time (in milliseconds) between the rotation and the fixed modes
//-----------------------------
rEncoder(N, sp, a, it) = thisEncoder
with {
    basicEncoder(sig, angle) = ho.encoder(N, sig, angle);
    thisEncoder = (_, rotationOrStaticAngle) : basicEncoder
    with {
        //converting the static angle from radians to [0; 1]
        an = (a / (2 * ma.PI), 1) : fmod;
        rotationOrStaticAngle = ((1-vn) * x + vn * an) * 2 * ma.PI;
        //to manage the case where frequency is zero, smoothly switches from one mode to another//
        vn = (sp == 0) : si.smooth(ba.tau2pole(it));
        x = (os.phasor(1, sp), an, 1) : (+, _) : fmod;
    };
};


//-------`(ho.)stereoEncoder`----------
// Encoding of a stereo pair of channels with symetric angles (a/2, -a/2).
//
// #### Usage
//
// ```
// _,_ : stereoEncoder(N, a) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `a` : opening angle in radians, left channel at a/2 angle, right channel at -a/2 angle
//-----------------------------
stereoEncoder(N, a) = (leftEncoder, rightEncoder) :> si.bus(2*N+1)
with {
    basicEncoder(sig, angle) = ho.encoder(N, sig, angle);
    leftEncoder = (_, a / 2) : basicEncoder;
    rightEncoder = (_, -a /2) : basicEncoder;
};


//-------`(ho.)multiEncoder`----------
// Encoding of a set of P signals distributed on the unit circle according to a list of P speeds and P angles.
//
// #### Usage
//
// ```
// _,_, ... : multiEncoder(N, lspeed, langle, it) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `lspeed` : a list of P speeds in turns by second (one speed per input signal, positive or negative)
// * `langle` : a list of P angles in radians on the unit circle to localize the sources (one angle per input signal)
// * `it` : interpolation time (in milliseconds) between the rotation and the fixed modes.
//-----------------------------
multiEncoder(N, lspeed, langle, it) = par(i, P, thisEncoder(ba.take(i+1, lspeed), ba.take(i+1, langle), it)) :> si.bus(2*N+1)
with {
    P = outputs(langle); //supposed to be the same as outputs(lspeed)
    basicEncoder(sig, angle) = ho.encoder(N, sig, angle);
    thisEncoder(sp, a, it) = (_, rotationOrStaticAngle) : basicEncoder
    with {
        //converting the static angle from radians to [0; 1]
        an = (a / (2 * ma.PI), 1) : fmod;
        rotationOrStaticAngle = ((1-vn) * x + vn * an) * 2 * ma.PI;
        //to manage the case where frequency is zero, smoothly switches from one mode to another//
        vn = (sp == 0) : si.smooth(ba.tau2pole(it));
        x = (os.phasor(1, sp), an, 1) : (+, _) : fmod;
    };
};


//--------------------------`(ho.)decoder`--------------------------------
// Decodes an ambisonics sound field for a circular array of loudspeakers.
//
// #### Usage
//
// ```
// _ : decoder(N, P) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `P`: the number of speakers (constant numerical expression)
//
// #### Note
//
// The number of loudspeakers must be greater or equal to 2n+1.
// It's preferable to use 2n+2 loudspeakers.
//-------------------------------------------------------------------
decoder(N, P) = par(i, 2*N+1, _) <: par(i, P, speaker(N, 2 * ma.PI*i/P))
with {
    speaker(N,a) = /(2), par(i, 2*N, _), encoder(N, 2/P, a) : si.dot(2*N+1);
};


//-----------------------`(ho.)decoderStereo`------------------------
// Decodes an ambisonic sound field for stereophonic configuration.
// An "home made" ambisonic decoder for stereophonic restitution
// (30° - 330°): Sound field lose energy around 180°. You should
// use `inPhase` optimization with ponctual sources.
// #### Usage
//
// ```
// _ : decoderStereo(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//--------------------------------------------------------------
decoderStereo(N) = decoder(N, P) <: (par(i, 2*N+2, gainLeft(360 * i / P)) :> _),
    (par(i, 2*N+2, gainRight(360 * i / P)) :> _)
with {
    P = 2*N+2;

    gainLeft(a) = _ * sin(ratio_minus + ratio_cortex)
    with {
        ratio_minus = ma.PI*.5 * abs((30 + a) / 60 * ((a <= 30)) + (a - 330) / 60 * (a >= 330));
        ratio_cortex= ma.PI*.5 * abs((120 + a) / 150 * (a > 30) * (a <= 180));
    };

    gainRight(a) = _ * sin(ratio_minus + ratio_cortex)
    with {
        ratio_minus = ma.PI*.5 * abs((390 - a) / 60 * (a >= 330) + (30 - a) / 60 * (a <= 30));
        ratio_cortex= ma.PI*.5 * abs((180 - a) / 150 * (a < 330) * (a >= 180));
    };
};


//-------`(ho.)iBasicDecoder`----------
// The irregular basic decoder is a simple decoder that projects the incoming ambisonic situation
// to the loudspeaker situation (P loudspeakers) whatever it is, without compensation.
// When there is a strong irregularity, there can be some discontinuity in the sound field.
//
// #### Usage
//
// ```
// _,_, ... : iBasicDecoder(N,la, direct, shift) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (there are 2*N+1 inputs to this function)
// * `la` : the list of P angles in degrees, for instance (0, 85, 182, 263) for four loudspeakers
// * `direct`: 1 for direct mode, -1 for the indirect mode (changes the rotation direction)
// * `shift` : angular shift in degrees to easily adjust angles
//-----------------------------
iBasicDecoder(N, la, direct, shift) = (par(i, 2*N+1, _) <: par(i, P, speaker(N, ang(i))))
with {
    P = outputs(la);
    ang(i) = (ba.take(i+1, la)  - direct * shift) * direct * ma.PI / 180.;
    speaker(N,alpha) = /(2), par(i, 2*N, _), ho.encoder(N,2/P,alpha) : si.dot(2*N+1);
};


//-------`(ho.)circularScaledVBAP`----------
// The function provides a circular scaled VBAP with all loudspeakers and the virtual source on the unit-circle.
//
// #### Usage
//
// ```
// _ : circularScaledVBAP(l, t) : _,_, ...
// ```
//
// Where:
//
// * `l` : the list of angles of the loudspeakers in degrees, for instance (0, 85, 182, 263) for four loudspeakers
// * `t` : the current angle of the virtual source in degrees
//-----------------------------
circularScaledVBAP(l, t) = thisCircularVbap
with {
    //modulo indexes between 1 and the number of elements of the list
    modIndex(i, l) = ma.modulo(i, outputs(l)) + 1;
    //
    //pick up the ith angle with a 360 degree modulo
    getElt(i, l) = ma.modulo(ba.take(modIndex(i, l), l), 360);
    //
    //function to compute the sinus of the difference between angles expressed in degrees
    diffSin(u, v) = sin((v - u) * ma.PI / 180.);
    //
    //permutations to be used to compute scaledVBAPGain
    p1(a, b, c, d) = (b, c, d, a);
    p2(a, b, c, d) = (a, c, b, d);
    //
    //computation of the scaled VBAP gain of a pair
    scaledVBAPGain(t1, t2, t) = ((diffSin(t2, t) <:(_, _, _)), (ma.signum(diffSin(t2, t1)) <: (_, _)), (diffSin(t, t1) <:(_, _, _))) : (*, *, *, *) : p1 : (_, _, (+ : sqrt <: (_, _))) : p2 : (/, /);
    sVBAPGain(i, l, t) = scaledVBAPGain(getElt(i, l), getElt(i+1, l), t);
    //
    //computes the left and the right gains using the matrix inversion (VBAP)
    leftGain(i, l, t) = sVBAPGain(i, l, t) : (_, !);
    rightGain(i, l, t) = sVBAPGain(i, l, t) : (!, _);
    //computation of boolean activePair that determines whether the pair of LS is active or not
    //we have to distinguish leftGain >0 and rightGain >= 0
    //if we put >=0 for both, two pairs will be simultaneously active when theta is one of the loudspeaker angles in the list
    //if we put > 0 for both, all the pairs will be inactive when theta is one of the loudspeaker angles in the list
    activePair(i, l, t) = (leftGain(i, l, t) > 0) * (rightGain(i, l, t) >= 0);
    //
    //computes the total gain for each loudspeaker
    cumulatedGain(i, l, t) = rightGain(outputs(l)+i-1, l, t) * activePair(outputs(l)+i-1, l, t) + leftGain(i, l, t) * activePair(i, l, t);
    //
    thisCircularVbap = _ <: par(i, outputs(l), *(cumulatedGain(i, l, t)));
};


//-------`(ho.)imlsDecoder`----------
// Irregular decoder in 2D for an irregular configuration of P loudspeakers
// using 2D VBAP for compensation.
//
// #### Usage
//
// ```
// _,_, ... : imlsDecoder(N,la, direct, shift) : _,_, ...
// ``` 
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `la` : the list of P angles in degrees, for instance (0, 85, 182, 263) for four loudspeakers
// * `direct`: 1 for direct mode, -1 for the indirect mode (changes the rotation direction)
// * `shift` : angular shift in degrees to easily adjust angles
//-----------------------------
imlsDecoder(N, la, direct, shift) = si.bus(2*N+1) : iVBAPDecoder
with {
    P = outputs(la);
    //The VBAP decoder uses VBAP compensation: it balances the regular decoder output enabling to use irregular angular setup.
    Q = max(2*N+2, P);
    iVBAPDecoder = ho.decoder(N, Q) : par(i, Q, circularScaledVBAP(la, (i * 360 / Q - direct * shift) * direct)) :> si.bus(P);
};


//-------`(ho.)iDecoder`----------
// General decoder in 2D enabling an irregular multi-loudspeaker configuration
// and to switch between multi-channel and stereo.
//
// #### Usage
//
// ```
// _,_, ... : iDecoder(N, la, direct, st, g) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `la`: the list of angles in degrees
// * `direct`: 1 for direct mode, -1 for the indirect mode (changes the rotation direction)
// * `shift` : angular shift in degrees to easily adjust angles
// * `st`: 1 for stereo, 0 for multi-loudspeaker configuration. When 1, stereo sounds goes through the first two channels
// * `g` : gain between 0 and 1
//-----------------------------
iDecoder(N, la, direct, shift, st, g) = thisDecoder
with {
    //p is the number of outputs
    P = outputs(la);
    ambi = 1 - st;
    //
    //for stereo decoding
    paddedStereoDecoder(N, P) = (gDecoderStereo, (0 <: si.bus(P-2)))
    with {
        leftDispatcher = _<:(*(1-direct), *(direct));
        rightDispatcher = _<:(*(direct), *(1-direct));
        gDecoderStereo = ho.decoderStereo(N) : (*(g), *(g)) : (leftDispatcher, rightDispatcher) :> (_,_);
    };
    //
    thisDecoder = si.bus(2*N+1) <: (si.bus(2*N+1), si.bus(2*N+1)) : (imlsDecoder(N, la, direct, shift), paddedStereoDecoder(N, P)) : (par(i, P, *(ambi)), *(st), *(st), si.bus(P-2))  :> si.bus(P) : par(i, P, *(g));
};


//============================Optimization Functions======================================
// Functions to weight the circular harmonics signals depending to the
// ambisonics optimization.
// It can be `basic` for no optimization, `maxRe` or `inPhase`.
//========================================================================================


//----------------`(ho.)optimBasic`-------------------------
// The basic optimization has no effect and should be used for a perfect
// circle of loudspeakers with one listener at the perfect center loudspeakers
// array.
//
// #### Usage
//
// ```
// _ : optimBasic(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//-----------------------------------------------------
optimBasic(N) = par(i, 2*N+1, _);


//----------------`(ho.)optimMaxRe`-------------------------
// The maxRe optimization optimizes energy vector. It should be used for an
// auditory confined in the center of the loudspeakers array.
//
// #### Usage
//
// ```
// _ : optimMaxRe(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//-----------------------------------------------------
optimMaxRe(N) = par(i, 2*N+1, optim(i, N, _))
with {
    optim(i, N, _)= _ * cos(indexabs / (2*N+1) * ma.PI)
    with {
        numberOfharmonics = 2 * N + 1;
        indexabs = (int)((i - 1) / 2 + 1);
    };
};


//----------------`(ho.)optimInPhase`-------------------------
//  The inPhase optimization optimizes energy vector and put all loudspeakers signals
// in phase. It should be used for an auditory.
//
// #### Usage
//
// ```
// _ : optimInPhase(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//-----------------------------------------------------
optimInPhase(N) = par(i, 2*N+1, optim(i, N, _))
with {
    optim(i, N, _)= _ * (fact(N)^2.) / (fact(N+indexabs) * fact(N-indexabs))
    with {
        indexabs = (int)((i - 1) / 2 + 1);
        fact(0) = 1;
        fact(n) = n * fact(n-1);
    };
};


//-------`(ho.)optim`----------
// Ambisonic optimizer including the three elementary optimizers:
// `(ho).optimBasic`, `(ho).optimMaxRe` and `(ho.)optimInPhase`.
//
// #### Usage
//
// ```
// _,_, ... : optim(N, ot) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `ot` : optimization type (0 for `optimBasic`, 1 for `optimMaxRe`, 2 for `optimInPhase`)
//-----------------------------
optim(N, ot) = thisOptimizer
with {
    optb = (ot == 0) : si.smoo;
    optm = (ot == 1) : si.smoo;
    opti = (ot == 2) : si.smoo;
    thisOptimizer = ((si.bus(2*N+1) <: ((si.bus(2*N+1):ho.optimBasic(N)), (si.bus(2*N+1):ho.optimMaxRe(N)), (si.bus(2*N+1):ho.optimInPhase(N)))), ((optb <: si.bus(2*N+1)), (optm <: si.bus(2*N+1)), (opti <: si.bus(2*N+1)))) : ro.interleave(6*N+3, 2) : par(i, 6*N+3, *) :> si.bus(2*N+1);
};


//----------------`(ho.)wider`-------------------------
// Can be used to wide the diffusion of a localized sound. The order
// depending signals are weighted and appear in a logarithmic way to
// have linear changes.
//
// #### Usage
//
// ```
// _ : wider(N,w) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `w`: the width value between 0 - 1
//-----------------------------------------------------
wider(N, w) = par(i, 2*N+1, perform(N, w, i, _))
with {
    perform(N, w, i, _) = _ * (log(N+1) * (1 - w) + 1) * clipweight
    with {
        clipweight = weighter(N, w, i) * (weighter(N, w, i) > 0) * (weighter(N, w, i) <= 1) + (weighter(N, w, i) > 1)
        with {
            weighter(N, w, 0) = 1.;
            weighter(N, w, i) = (((w * log(N+1)) - log(indexabs)) / (log(indexabs+1) - log(indexabs)))
            with {
                indexabs = (int)((i - 1) / 2 + 1);
            };
        };
    };
};


//-------`(ho.)mirror`----------
// Mirroring effect on the sound field.
//
// #### Usage
//
// ```
// _,_, ... : mirror(N, fa) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `fa` : mirroring type (1 = original sound field, 0 = original+mirrored sound field, -1 = mirrored sound field)
//-----------------------------
mirror(N, fa) = (*(1), par(i, N, (*(fa), *(1))));


//----------------`(ho.)map`-------------------------
// It simulates the distance of the source by applying a gain
// on the signal and a wider processing on the soundfield.
//
// #### Usage
//
// ```
// map(N, x, r, a)
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `x`: the signal
// * `r`: the radius
// * `a`: the angle in radian
//-----------------------------------------------------
map(N, x, r, a) = encoder(N, x * volume(r), a) : wider(N, ouverture(r))
with {
    volume(r) = 1. / (r * r * (r > 1) + (r <= 1));
    ouverture(r) = r * (r < 1) + (r >= 1);
};


//----------------`(ho.)rotate`-------------------------
// Rotates the sound field.
//
// #### Usage
//
// ```
// _ : rotate(N, a) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `a`: the angle in radian
//-----------------------------------------------------
rotate(N, a) = par(i, 2*N+1, _) <: par(i, 2*N+1, rotation(i, a))
with {
    rotation(i, a) = (par(j, 2*N+1, gain1(i, j, a)), par(j, 2*N+1, gain2(i, j, a)), par(j, 2*N+1, gain3(i, j, a)) :> _)
    with {
        indexabs = (int)((i - 1) / 2 + 1);
        gain1(i, j, a) = _ * cos(a * indexabs) * (j == i);
        gain2(i, j, a) = _ * sin(a * indexabs) * (j-1 == i) * (j != 0) * (i%2 == 1);
        gain3(i, j, a) = (_ * sin(a * indexabs)) * (j+1 == i) * (j != 0) * (i%2 == 0) * (-1);
    };
};


//-------`(ho.)scope`----------
// Produces an XY pair of signals representing the ambisonic sound field.
//
// #### Usage
//
// ```
// _,_, ... : scope(N, rt) : _,_
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `rt` : refreshment time in milliseconds
//-----------------------------
scope(N, rt) = thisScope
with {
    //Angle sweeping at a speed corresponding to refresh period between 0 and 2*PI
    theta = os.phasor(1, 1/rt) * 2 * ma.PI;
    //we get the vector of harmonic functions thanks to the encoding function//
    harmonicsVector = ho.encoder(N, 1, theta);
    //
    normalizedVector(N) = si.bus(N) <: (si.bus(N), norm) : ro.interleave(N, 2) : par(i, N, /)
    with {
        norm = par(i, N, _ <:(_,_) : *) :> _ : sqrt <: ((_ == 0), (_ > 0), _) : (_,*) : + <: si.bus(N);
    };
    //building (2N+1) normalized vectors
    inputVector = (*(0.5), par(i, (2*N), _)) : normalizedVector(2*N+1);
    normalizedHarmonics = harmonicsVector : normalizedVector(2*N+1);
    //
    rho = (inputVector, normalizedHarmonics) : si.dot(2*N+1) ;
    thisScope = (rho <: (ma.fabs, (_ >= 0))) : ((_ <: (_,_)), _) : (*(sin(theta)), *(cos(theta)), _) : (*(-1), _,_);
};


//============================Spatial Sound Processes ====================================
// We propose implementations of processes intricated to the ambisonic model.
// The process is implemented using as many instances as the number of harmonics at at certain order.
// The key control parameters of these instances are computed thanks to distribution functions
// (th functions below) and to a global driving factor.
//========================================================================================

//-------`(ho.).fxDecorrelation`----------
// Spatial ambisonic decorrelation in fx mode.
//
// `fxDecorrelation` applies decorrelations to spatial components already created.
// The decorrelation is defined for each #i spatial component among P=2\*N+1 at the ambisonic order `N`
// as a delay of 0 if factor `fa` is under a certain value 1-(i+1)/P and d\*F((i+1)/p) in the contrary case,
// where `d` is the maximum delay applied (in samples) and F is a distribution function for durations.
// The user can choose this delay time distribution among 22 different ones.
// The delay increases according to the index of ambisonic components.
// But it increases at each step and it is modulated by a threshold.
// Therefore, delays are progressively revealed when the factor increases:
//
// * when the factor is close to 0, only upper components are delayed;
// * when the factor increases, more and more components are delayed.
//
//H                 THRESHOLD            DELAY
//0                 1-1/P                0 OR DELAY*F(1/P)
//-1                1-2/P                0 OR DELAY*F(2/P)
//1                 1-3/P                0 OR DELAY*F(3/P)
//-2                1-4/P                0 OR DELAY*F(4/P)
//2                 1-5/P                0 OR DELAY*F(5/P)
//...
//-(N-1)            1-(P-3)/P            0 OR DELAY*F((P-3)/P)
//(N-1)             1-(P-2)/P            0 OR DELAY*F((P-2)/P)
//-N                1-(P-1)/P            0 OR DELAY*F((P-1)/P)
//N                 1-P/P                0 OR DELAY*F(P/P)
//
//
// #### Usage
//
// ```
// _,_, ... : fxDecorrelation(N, d, wf, fa, fd, tf) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `d`: the maximum delay applied (in samples)
// * `wf`: window frequency (in Hz) for the overlapped delay
// * `fa`: decorrelation factor (between 0 and 1)
// * `fd`: feedback / level of reinjection (between 0 and 1)
// * `tf`: type of function of delay distribution (integer, between 0 and 21)
//-----------------------------
fxDecorrelation(N, d, wf, fa, fd, tf) = par(i, 2*N+1, gate(d, i, 2*N+1, fa, tf, wf, fd))
with {
    gate(d, i, N, fa, tf, wf, fd) = _ <: fdOverlappedDelay(dur(d, i, N, fa, tf), 262144, wf, fd) * env1(fa, i, N), _ * env1c(fa, i, N) : +;
    //
    fdOverlappedDelay(nsamp, nmax, freq, fdbk) = (+ : de.sdelay(nmax, int(ma.SR / freq), nsamp)) ~ (*(fdbk));
    //
    env1(fa, i, N) = (fa > ((N-i-1)/N)) : si.smooth(ba.tau2pole(0.005));
    env1c(fa, i, N) = 1 - env1(fa, i, N);
    //
    //computes the ith duration of the ith delay in samples with twenty two possibilities of distribution
    elemdur(d, i, p, fa, tf, ind) = (tf == ind) * (fa > (1 - x)) * d * x * fa
    with {
        x = th(ind, i, p);
    };
    //duration in samples computed as a sum of the 22 cases//
    dur(d, i, p, fa, tf) = sum(ind, 22, elemdur(d, i, p, fa, tf, ind)) : int;
};

//-------`(ho.).synDecorrelation`----------
// Spatial ambisonic decorrelation in syn mode.
//
// `synDecorrelation` generates spatial decorrelated components in ambisonics from one mono signal.
// The decorrelation is defined for each #i spatial component among P=2\*N+1 at the ambisonic order `N`
// as a delay of 0 if factor `fa` is under a certain value 1-(i+1)/P and d\*F((i+1)/p) in the contrary case,
// where `d` is the maximum delay applied (in samples) and F is a distribution function for durations.
// The user can choose this delay time distribution among 22 different ones.
// The delay increases according to the index of ambisonic components.
// But it increases at each step and it is modulated by a threshold.
// Therefore, delays are progressively revealed when the factor increases:
//
// * when the factor is close to 0, only upper components are delayed;
// * when the factor increases, more and more components are delayed.
//
// When the factor is between [0; 1/P], upper harmonics are progressively faded and the level of the H0 component is compensated
// to avoid source localization and to produce a large mono.
//
//H                THRESHOLD            DELAY
//0                1-1/P                0 OR DELAY*F(1/P)
//-1               1-2/P                0 OR DELAY*F(2/P)
//1                1-3/P                0 OR DELAY*F(3/P)
//-2               1-4/P                0 OR DELAY*F(4/P)
//2                1-5/P                0 OR DELAY*F(5/P)
//...
//-(N-1)           1-(P-3)/P            0 OR DELAY*F((P-3)/P)
//(N-1)            1-(P-2)/P            0 OR DELAY*F((P-2)/P)
//-N               1-(P-1)/P            0 OR DELAY*F((P-1)/P)
//N                1-P/P                0 OR DELAY*F(P/P)
//
//
// #### Usage
//
// ```
// _,_, ... : synDecorrelation(N, d, wf, fa, fd, tf) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `d`: the maximum delay applied (in samples)
// * `wf`: window frequency (in Hz) for the overlapped delay
// * `fa`: decorrelation factor (between 0 and 1)
// * `fd`: feedback / level of reinjection (between 0 and 1)
// * `tf`: type of function of delay distribution (integer, between 0 and 21)
//-----------------------------
synDecorrelation(N, d, wf, fa, fd, tf) =  _ <: par(i, 2*N+1, crossFade(d, i, 2*N+1, fa, tf, wf, fd))
with {
    crossFade(d, i, N, fa, tf, wf, fd) = _ <: fdOverlappedDelay(dur(d, i, N, fa, tf), 262144, wf, fd) * env1(fa, i, N), _ * env1c(fa, i, N) :> _ * env2(fa, i, N);
    //
    fdOverlappedDelay(nsamp, nmax, freq, fdbk) = (+ : de.sdelay(nmax, int(ma.SR / freq), nsamp)) ~ (*(fdbk));
    //
    env1(fa, i, N) = (fa > ((N-i-1)/N)) : si.smooth(ba.tau2pole(0.005));
    env1c(fa, i, N) = 1 - env1(fa, i, N) ;
    env2(fa, i, N) = ((i > 0) * N * min(fa, 1/N)) + ((i == 0) * (sqrt(N) * (1 - (N - sqrt(N)) * min(fa, 1/N)))) : si.smooth(ba.tau2pole(0.005));
    //
    //computes the ith duration of the ith delay in samples with twenty two possibilities of distribution
    elemdur(d, i, p, fa, tf, ind) = (tf == ind) * fa * d * x
    with {
        x = th(ind, i, p);
    };
    //duration in samples computed as a sum of the 22 cases//
    dur(d, i, p, fa, tf) = sum(ind, 22, elemdur(d, i, p, fa, tf, ind)) : int;
};

//-------`(ho.).fxRingMod`----------
// Spatial ring modulation in syn mode.
//
// `fxRingMod` applies ring modulation to spatial components already created.
// The ring modulation is defined for each spatial component among P=2\*n+1 at the ambisonic order `N`.
// For each spatial component #i, the result is either the original signal or a ring modulated signal
// according to a threshold that is i/P.
//
// The general process is drive by a factor `fa` between 0 and 1 and a modulation frequency `f0`.
// If `fa` is greater than theshold (P-i-1)/P, the ith ring modulator is on with carrier frequency of f0\*(i+1)/P.
// On the contrary, it provides the original signal.
//
// Therefore ring modulators are progressively revealed when `fa` increases.
//
//H                THRESHOLD                OUTPUT
//0                (P-1)/P                  ORIGINAL OR RING MODULATION BY F0*1/P
//-1               (P-2)/P                  ORIGINAL OR RING MODULATION BY F0*2/P
//1                (P-3)/P                  ORIGINAL OR RING MODULATION BY F0*3/P
//-2               (P-4)/P                  ORIGINAL OR RING MODULATION BY F0*4/P
//2                (P-5)/P                  ORIGINAL OR RING MODULATION BY F0*5/P
//...
//-(N-1)           3/P                      ORIGINAL OR RING MODULATION BY F0*(P-3)/P
//(N-1)            2/P                      ORIGINAL OR RING MODULATION BY F0*(P-2)/P
//-N               1/P                      ORIGINAL OR RING MODULATION BY F0*(P-1)/P
//N                0                        ORIGINAL OR RING MODULATION BY F0*P/P=F0
//
//
// #### Usage
//
// ```
// _,_, ... : fxRingMod(N, f0, fa, tf) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `f0`: the maximum delay applied (in samples)
// * `fa`: decorrelation factor (between 0 and 1)
// * `tf`: type of function of delay distribution (integer, between 0 and 21)
//-----------------------------
fxRingMod(N, f0, fa, tf) = par(i, 2*N+1, gate_ringmod(f0, i, 2*N+1, fa, tf))
with {
    //
    env1(fa, i, N) = (fa > ((N-i-1)/N)) : si.smooth(ba.tau2pole(0.005));
    env1c(fa, i, N) = 1 - env1(fa, i, N);
    //
    gate_ringmod(f, i, N, fa, tf) = _ <: _ * os.osccos(freq(f, i, N, tf)) * env1(fa, i, N), _ * env1c(fa, i, N) : +;
    //
    ringmodfreq(f, i, N, tf, ind) = (tf == ind) * f * x * coef
    with {
        x = th(ind, i, N);
        coef = min(1, max(N * (fa - (N - i - 1) / N), 0));
    };
    //
    freq(f, i, N, tf) = sum(ind, 22, ringmodfreq(f, i, N, tf, ind)) : int;
};

//-------`(ho.).synRingMod`----------
// Spatial ring modulation in syn mode.
//
// `synRingMod` generates spatial components in ambisonics from one mono signal thanks to ring modulation.
// The ring modulation is defined for each spatial component among P=2\*n+1 at the ambisonic order `N`.
// For each spatial component #i, the result is either the original signal or a ring modulated signal
// according to a threshold that is i/P.
//
// The general process is drive by a factor `fa` between 0 and 1 and a modulation frequency `f0`.
// If `fa` is greater than theshold (P-i-1)/P, the ith ring modulator is on with carrier frequency of f0\*(i+1)/P.
// On the contrary, it provides the original signal.
//
// Therefore ring modulators are progressively revealed when `fa` increases.
// When the factor is between [0; 1/P], upper harmonics are progressively faded and the level of the H0 component is compensated
// to avoid source localization and to produce a large mono.
//
//H                 THRESHOLD                OUTPUT
//0                 (P-1)/P                  ORIGINAL OR RING MODULATION BY F0*1/P
//-1                (P-2)/P                  ORIGINAL OR RING MODULATION BY F0*2/P
//1                 (P-3)/P                  ORIGINAL OR RING MODULATION BY F0*3/P
//-2                (P-4)/P                  ORIGINAL OR RING MODULATION BY F0*4/P
//2                 (P-5)/P                  ORIGINAL OR RING MODULATION BY F0*5/P
//...
//-(N-1)            3/P                      ORIGINAL OR RING MODULATION BY F0*(P-3)/P
//(N-1)             2/P                      ORIGINAL OR RING MODULATION BY F0*(P-2)/P
//-N                1/P                      ORIGINAL OR RING MODULATION BY F0*(P-1)/P
//N                 0                        ORIGINAL OR RING MODULATION BY F0*P/P=F0
//
//
// #### Usage
//
// ```
// _,_, ... : synRingMod(N, f0, fa, tf) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `f0`: the maximum delay applied (in samples)
// * `fa`: decorrelation factor (between 0 and 1)
// * `tf`: type of function of delay distribution (integer, between 0 and 21)
//-----------------------------
synRingMod(N, f0, fa, tf) = _ <: par(i, 2*N+1, crossfade_ringmod(f0, i, 2*N+1, fa, tf))
with {
    //
    env1(fa, i, N) = (fa > ((N-i-1)/N)) : si.smooth(ba.tau2pole(0.005));
    env1c(fa, i, N) = 1 - env1(fa, i, N);
    env2(fa, i, N) = ((i > 0) * N * min(fa, 1/N)) + ((i == 0) * (sqrt(N) * (1 - (N - sqrt(N)) * min(fa, 1/N)))) : si.smooth(ba.tau2pole(0.005));
    //
    crossfade_ringmod(f, i, N, fa, tf) = _ <: _ * os.osccos(freq(f, i, N, tf)) * env1(fa, i, N), _ * env1c(fa, i, N) :> _ * env2(fa, i, N);
    //
    ringmodfreq(f, i, N, tf, ind) = (tf == ind) * f * x * coef
    with {
        x = th(ind, i, N);
        coef = min(1, max(N * (fa - (N - i - 1) / N), 0));
    };
    //
    freq(f, i, N, tf) = sum(ind, 22, ringmodfreq(f, i, N, tf, ind)) : int;
};


//TYPES OF DISTRIBUTIONS: 22 EASING FUNCTIONS FROM [0, 1] to [0,1]
//(i+1)/p belongs to [0, 1] and its image by any function in the list also belongs to the interval

th(0, i, p) = (i+1) / p;
th(1, i, p) = ((i+1) / p)^2;
th(2, i, p) = sin(ma.PI * 0.5 * (i+1) / p);
th(3, i, p) = log10(1 + (i+1) / p) / log10(2);
th(4, i, p) = sqrt((i+1) / p);
th(5, i, p) = 1 - cos(ma.PI * 0.5 * (i+1) / p);
th(6, i, p) = (1 - cos(ma.PI * (i+1) / p)) * 0.5;
th(7, i, p) = 1 - (1 - (i+1) / p )^2;
th(8, i, p) = ((i+1) / p < 0.5) * 2 * ((i+1) / p)^2 + ((i+1) / p >= 0.5) * (1 - (-2 * (i+1) / p + 2)^2 * 0.5);
th(9, i, p) = ((i+1) / p)^3;
th(10, i, p) = 1 - (1 - (i+1) / p)^3;
th(11, i, p) = ((i+1) / p < 0.5) * 4 * ((i+1) / p)^3 + ((i+1) / p >= 0.5) * (1 - (-2 * (i+1) / p + 2)^3 * 0.5);
th(12, i, p) = ((i+1) / p)^4;
th(13, i, p) = 1 - (1 - (i+1) / p)^4;
th(14, i, p) = ((i+1) / p < 0.5) * 8 * ((i+1) / p)^4 + ((i+1) / p >= 0.5) * (1 - (-2 * (i+1) / p + 2)^4 * 0.5);
th(15, i, p) = ((i+1) / p)^5;
th(16, i, p) = 1 - (1 - (i+1) / p)^5;
th(17, i, p) = ((i+1) / p < 0.5) * 16 * ((i+1) / p)^5 + ((i+1) / p >= 0.5) * (1 - (-2 * (i+1) / p + 2)^5 * 0.5);
th(18, i, p) = 2^(10 * (i+1) / p - 10);
th(19, i, p) = ((i+1) / p < 1) * (1 - 2^(-10 * (i+1) / p)) + ((i+1) / p == 1);
th(20, i, p) = 1 - sqrt(1 - ((i+1) / p)^2);
th(21, i, p) = sqrt(1 - ((i+1) / p - 1)^2);


//========================================================================================
//============================3D Functions================================================
//========================================================================================
//========================================================================================

//----------------------`(ho.)encoder3D`---------------------------------
// Ambisonic encoder. Encodes a signal in the circular harmonics domain
// depending on an order of decomposition, an angle and an elevation.
//
// #### Usage
//
// ```
// encoder3D(N, x, a, e) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `x`: the signal
// * `a`: the angle
// * `e`: the elevation
//----------------------------------------------------------------
encoder3D(N, x, theta, phi) = par(i, (N+1) * (N+1), x * y(degree(i), order(i), theta, phi))
with {
    // The degree l of the harmonic[l, m]
    degree(index) = int(sqrt(index));
    // The order m of the harmonic[l, m]
    order(index) = int(index - int(degree(index) * int(degree(index) + 1)));

    // The spherical harmonics
    y(l, m, theta, phi) =  e(m, theta2) * k(l, m) * p(l, m, cos(phi + ma.PI * 0.5))
    with {
        //theta2 enables a continuous movement of elevation (when phi becomes greater than Pi/2)
        theta2 = theta + (1 - int(fmod(fmod(phi / ma.PI - 0.5, 2) + 2, 2))) * ma.PI;
        //
        // The associated Legendre polynomial
        // If l = 0   => p = 1
        // If l = m   => p = -1 * (2 * (l-1) + 1) * sqrt(1 - cphi*cphi) * p(l-1, l-1, cphi)
        // If l = m+1 => p = phi * (2 * (l-1) + 1) * p(l-1, l-1, cphi)
        // Else => p = (cphi * (2 * (l-1) + 1) * p(l-1, abs(m), cphi) - ((l-1) + abs(m)) * p(l-2, abs(m), cphi)) / ((l-1) - abs(m) + 1)
        p(l, m, cphi) = pcalcul(((l != 0) & (l == abs(m))) + ((l != 0) & (l == abs(m)+1)) * 2 + ((l != 0) & (l != abs(m)) & (l != abs(m)+1)) * 3, l, m, cphi)
        with {
            pcalcul(0, l, m, cphi) = 1;
            pcalcul(1, l, m, cphi) = -1 * (2 * (l-1) + 1) * sqrt(1 - cphi*cphi) * p(l-1, l-1, cphi);
            pcalcul(2, l, m, cphi) = cphi * (2 * (l-1) + 1) * p(l-1, l-1, cphi);
            pcalcul(s, l, m, cphi) = (cphi * (2 * (l-1) + 1) * p(l-1, abs(m), cphi) - ((l-1) + abs(m)) * p(l-2, abs(m), cphi)) / ((l-1) - abs(m) + 1);
        };

        // The exponential imaginary
        // If m > 0 => e^i*m*theta = cos(m * theta)
        // If m < 0 => e^i*m*theta = sin(-m * theta)
        // If m = 0 => e^i*m*theta = 1
        e(m, theta) = ecalcul((m > 0) * 2 + (m < 0), m, theta)
        with {
            ecalcul(2, m, theta) = cos(m * theta);
            ecalcul(1, m, theta) = sin(abs(m) * theta);
            ecalcul(s, m, theta) = 1;
        }; 
        
        // The normalization
        // If m  = 0 => k(l, m) = 1
        // If m != 0 => k(l, m) = sqrt((l - abs(m))! / l + abs(m))!) * sqrt(2)
        k(l, m) = kcalcul((m != 0), l, m)
        with {
            kcalcul(0, l, m) = 1;
            kcalcul(1, l, m) = sqrt(2) / sqrtFactQuotient(l+abs(m), l-abs(m))
            with {
                //factorial quotient fq(n, p)=n! / p! = n(n-1)...(p+1) when n > p
                //enables factor simplification
                //and considering the square root of a product as a product of square roots
                sqrtFactQuotient(n, p) = sqrtProd(n-p, p)
                with {
                    //sqrtProd(n, p) computes the product sqrt(p+1) x sqrt(p+2) x ... x sqrt(n)
                    //to enable factorial quotient simplification
                    sqrtProd(1, p) = sqrt(p+1);
                    sqrtProd(n, p) = sqrt(p+n) * sqrtProd(n-1, p);
                };
            };
        };
    };
};


//-------`(ho.)rEncoder3D`----------
// Ambisonic encoder in 3D including source rotation. A mono signal is encoded at at certain ambisonic order
// with two possible modes: either rotation with 2 angular speeds (azimuth and elevation), or static with a fixed pair of angles.
//
// `rEncoder3D` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : rEncoder3D(N, azsp, elsp, az, el, it) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `azsp`: the azimuth speed expressed as angular speed (2PI/sec), positive or negative
// * `elsp`: the elevation speed expressed as angular speed (2PI/sec), positive or negative
// * `az`: the fixed azimuth when the azimuth rotation stops (azsp = 0) in radians
// * `el`: the fixed elevation when the elevation rotation stops (elsp = 0) in radians
// * `it` : interpolation time (in milliseconds) between the rotation and the fixed modes
//-----------------------------
rEncoder3D(N, azsp, elsp, az, el, it) = this3DEncoder
with {
    basic3DEncoder(sig, ang1, ang2) = encoder3D(N, sig, ang1, ang2);
    this3DEncoder = (_, rotationOrStaticAzim, rotationOrStaticElev) : basic3DEncoder
    with {
        x1 = (os.phasor(1, azsp), az, 1) : (+, _) : fmod : *(2 * ma.PI);
        vn1 = (azsp == 0) : si.smooth(ba.tau2pole(it));
        rotationOrStaticAzim = (1-vn1) * x1 + vn1 * az;
        x2 = (os.phasor(1, elsp), el, 1) : (+, _) : fmod : *(2 * ma.PI);
        vn2 = (elsp == 0) : si.smooth(ba.tau2pole(it));
        rotationOrStaticElev =  (1-vn2) * x2 + vn2 * el;
    };
};


//----------------`(ho.)optimBasic3D`-------------------------
// The basic optimization has no effect and should be used for a perfect
// sphere of loudspeakers with one listener at the perfect center loudspeakers
// array.
//
// #### Usage
//
// ```
// _ : optimBasic3D(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//-----------------------------------------------------
optimBasic3D(N) = par(i, (N+1) * (N+1), _);


//----------------`(ho.)optimMaxRe3D`-------------------------
// The maxRe optimization optimize energy vector. It should be used for an
// auditory confined in the center of the loudspeakers array.
//
// #### Usage
//
// ```
// _ : optimMaxRe3D(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//-----------------------------------------------------
optimMaxRe3D(N) = par(i, (N+1) * (N+1), MaxRe(N, degree(i), _))
with {
    // The degree l of the harmonic[l, m]
    degree(index)  = int(sqrt(index));
    MaxRe(N, l, _)= _ * cos(l / (2*N+2) * ma.PI);
};


//----------------`(ho.)optimInPhase3D`-------------------------
// The inPhase Optimization optimizes energy vector and put all loudspeakers signals
// in phase. It should be used for an auditory.
//
// #### Usage
//
// ```
// _ : optimInPhase3D(N) : _
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
//-----------------------------------------------------
optimInPhase3D(N) = par(i, (N+1) * (N+1), InPhase(N, degree(i), _))
with {
    // The degree l of the harmonic[l, m]
    degree(index)  = int(sqrt(index));
    InPhase(N, l, _)= _ * (fact(N) * fact(N)) / (fact(N - l) * fact(N + l))
    with {
        fact(0) = 1;
        fact(n) = n * fact(n-1);
    };
};


//-------`(ho.)optim3D`----------
// Ambisonic optimizer including the three elementary optimizers:
// `(ho).optimBasic3D`, `(ho).optimMaxRe3D` and `(ho.)optimInPhase3D`.
//
// #### Usage
//
// ```
// _,_, ... : optim3D(N, ot) : _,_, ...
// ```
//
// Where:
//
// * `N`: the ambisonic order (constant numerical expression)
// * `ot` : optimization type (0 for optimBasic, 1 for optimMaxRe, 2 for optimInPhase)
//-----------------------------
optim3D(N, ot) = thisOptimizer
with {
    optb = (ot == 0) : si.smoo;
    optm = (ot == 1) : si.smoo;
    opti = (ot == 2) : si.smoo;
    bus3D = si.bus((N+1)*(N+1));
    thisOptimizer = ((bus3D  <: ((bus3D:ho.optimBasic3D(N)), (bus3D:ho.optimMaxRe3D(N)), (bus3D:ho.optimInPhase3D(N)))), ((optb <: bus3D), (optm <: bus3D), (opti <: bus3D))) : ro.interleave(3*(N+1)*(N+1), 2) : par(i, 3*(N+1)*(N+1), *) :> bus3D;
};