//#################################### spats.lib ##########################################
// This library contains a collection of tools for sound spatialization.
// Its official prefix is `sp`.
//
// #### References
// * <https://github.com/grame-cncm/faustlibraries/blob/master/spats.lib>
//########################################################################################

/************************************************************************
 ************************************************************************
FAUST library file
Copyright (C) 2003-2025 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.

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

si = library("signals.lib");
ma = library("maths.lib");
fi = library("filters.lib");
de = library("delays.lib");

declare name "Faust Spatialization Library";
declare version "1.2.0";

//-----------------------`(sp.)panner`------------------------
// A simple linear stereo panner.
// `panner` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : panner(g) : _,_
// ```
//
// Where:
//
// * `g`: the panning (0-1)
//------------------------------------------------------------
panner(g) = _ <: *(1-g), *(g);

// TODO: need demo function for panner here


//---------------`(sp.)constantPowerPan`----------------------
// Apply the constant power pan rule to a stereo signal.
// The channels are not respatialized. Their gains are simply
// adjusted. A pan of 0 preserves the left channel and silences
// the right channel. A pan of 1 has the opposite effect.
// A pan value of 0.5 applies a gain of 0.5 to both channels. 
//
// #### Usage
//
// ```
// _,_ : constantPowerPan(p) : _,_
// ```
//
// Where:
//
// * `p`: the panning (0-1)
//------------------------------------------------------------
declare dryWetMixer author "David Braun";

constantPowerPan(p, x, y) = x*gainLeft, y*gainRight
with {
    theta = ma.PI*p/2.;
    gainLeft = cos(theta)/sqrt(2.);
    gainRight = sin(theta)/sqrt(2.);
};


//-----------------------`(sp.)spat`------------------------
// GMEM SPAT: n-outputs spatializer.
// `spat` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : spat(N,r,d) : si.bus(N)
// ```
//
// Where:
//
// * `N`: number of outputs (a constant numerical expression)
// * `r`: rotation (between 0 et 1)
// * `d`: distance of the source (between 0 et 1)
//------------------------------------------------------
declare spat author "Laurent Pottier, revised by Romain Michon";

spat(N,a,d)= _ <: par(i, N, *(scaler(i, N, a, d) : si.smooth(0.9999)))
with {
    scaler(i,N,a,d) = (d/2.0+0.5)
    * sqrt(max(0.0, 1.0 - abs(fmod(a+0.5+float(N-i)/N, 1.0) - 0.5) * N * d));
};


//-------`(sp.)wfs`-------------------
// Wave Field Synthesis algorithm for multiple sound sources. 
// Implementation generalized starting from Pierre Lecomte version. 
//
// #### Usage
//
// ```
// wfs(xref, yref, zref, speakersDist, nSources, nSpeakers, inGain, xs, ys, zs) : si.bus(nSpeakers)
// ```
//
// Where:
//
// * `xref`: x-coordinate of the reference listening point
// * `yref`: y-coordinate of the reference listening point
// * `zref`: z-coordinate of the reference listening point
// * `speakersDist`: distance between speakers
// * `nSources`: number of sound sources
// * `nSpeakers`: number of speakers
// * `inGain`: input gain (0-1) as a function of the source index
// * `xs`: x-coordinate of the sound source as a function of the source index
// * `ys`: y-coordinate of the sound source as a function of the source index
// * `zs`: z-coordinate of the sound source as a function of the source index
//--------------------------------------
declare wfs author "Pierre Lecomte, revised by GRAME";
wfs(xref, yref, zref, speakersDist, nSources, nSpeakers, inGain, xs, ys, zs) 
    = par(i, nSources, oneSource(xref,yref,zref,speakersDist,nSpeakers,inGain(i),xs(i),ys(i),zs(i))) 
    :> si.bus(nSpeakers)
with {

    // Physical constant: speed of sound in air (m/s)
    c = 340;

    // Filtering stage (pre-defined transfer function approximation)
    filter = fi.tf2(b0,b1,b2,a1,a2)
    with {
        b0 = 1;
        b1 = -0.8715;
        b2 = 0.0412;
        a1 = -0.31134;
        a2 = -0.088955;
    };

    // Function to compute Euclidean distance between two points in 3D space
    norm(x1, y1, z1, x2, y2, z2) = sqrt((x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2);

    // Driving function to compute delay and amplitude compensation
    driving(speakersDist, nSpeakers, pregain, r, rLong) = *(pregain) / r^2 : de.fdelay(rmax, d)
    with {
        rmax = (nSpeakers*speakersDist)/c*ma.SR;
        d = (r-rLong)/c*ma.SR;
    };

    // Function to compute the minimum value of N elements
    mins(1) = _;
    mins(2) = min;
    mins(N) = mins(int(N/2)),mins(N-int(N/2)) : min;

    // List defining the x-position of each speaker in the WFS setup
    speakersXpos(speakersDist, nSpeakers) = par(i,nSpeakers,(-speakersDist*nSpeakers/2 + speakersDist/2 + i*speakersDist));

    // Function to synthesize a single sound source for WFS
    oneSource(xref, yref, zref, speakersDist, nSpeakers, inGain, xs, ys, zs) =
        *(inGain)
        : filter <:
        (speakersXpos(speakersDist,nSpeakers) <: par(i,nSpeakers,norm(xref,yref,zref,_,0,0)*ys),(par(i,nSpeakers,norm(xs,ys,zs,_,0,0))
         <: par(i,nSpeakers,_),(mins(nSpeakers)
         <: par(i,nSpeakers,_)))),par(i,nSpeakers,_)
         : route(nSpeakers*4,nSpeakers*4,par(i,(nSpeakers*4)-1,(((i*nSpeakers))%(nSpeakers*4-1)+1,i+1)),(nSpeakers*4,nSpeakers*4))
         : par(i,nSpeakers,driving(speakersDist,nSpeakers));

};

//-------`(sp.)wfs_ui`-------------------
// Wave Field Synthesis algorithm for multiple sound sources with a built-in UI.
//
// #### Usage
//
// ```
// wfs_ui(xref, yref, zref, speakersDist, nSources, nSpeaker) : si.bus(nSpeakers)
// ```
//
// Where:
//
// * `xref`: x-coordinate of the reference listening point
// * `yref`: y-coordinate of the reference listening point
// * `zref`: z-coordinate of the reference listening point
// * `speakersDist`: distance between speakers
// * `nSources`: number of sound sources
// * `nSpeakers`: number of speakers
//
// #### Example test program
//
// ```
// // Distance between speakers in meters
// speakersDist = 0.0783;  
//
// // Reference listening point (central position for WFS)
// xref = 0;
// yref = 1;
// zref = 0;
//
// Spatialize 4 sound sources on 16 speakers
// process = wfs_ui(xref,yref,zref,speakersDist,4,16);
// ```
//--------------------------------------
declare wfs author "Pierre Lecomte, revised by GRAME";
wfs_ui(xref, yref, zref, speakersDist, nSources, nSpeakers) = wfs(xref, yref, zref, speakersDist, nSources, nSpeakers, inGain, xs, ys, zs)
with {
    xs(i) = hslider("v:Sound%2i/xs", 0, -10, 10, 0.01);
    ys(i) = hslider("v:Sound%2i/ys", 1.5,  1, 10, 0.01);
    zs(i) = hslider("v:Sound%2i/zs", 0, -10, 10, 0.01);
    inGain(i) = hslider("v:Sound%2i/inGain",1,0,10,0.01);
};

//---------------`(sp.)stereoize`-------------
// Transform an arbitrary processor `p` into a stereo processor with 2 inputs
// and 2 outputs.
//
// #### Usage
//
// ```
// _,_ : stereoize(p) : _,_
// ```
//
// Where:
//
// * `p`: the arbitrary processor
//----------------------------------------
stereoize(p) = S(inputs(p), outputs(p))
with {
	// degenerated processor with no outputs
	S(n,0) = !,! : 0,0; 		// just in case, probably a rare case
	// processors with no inputs
	S(0,1) = !,! : p <: _,_; 	// add two fake inputs and split output
	S(0,2) = !,! : p;
	S(0,n) = !,! : p,p :> _,_;	// we are sure this will work if n is odd
	// processors with one input
	S(1,1) = p,p; 				// add two fake inputs and split output
	S(1,n) = p,p :> _,_;		// we are sure this will work if n is odd
	// processors with two inputs
	S(2,1) = p <: _,_; 			// split the output
	S(2,2) = p; 				// nothing to do, p is already stereo
	// processors with inputs > 2 and outputs > 2
	S(n,m) = _,_ <: p,p :> _,_;	// we are sure this works if n or p are odd
};

// TODO: need demo function of spat here

//////////////////////////////////////////////////////////////////////////////////////////
// UNDOCUMENTED/DISMISSED ELEMENTS
//////////////////////////////////////////////////////////////////////////////////////////

// music.lib:
	// The following functions could remain available but they would have to be
	// factorized and reimplemented using the `par` function...
	// bus2 = _,_;
	// bus3 = _,_,_;
	// bus4 = _,_,_,_;
	// bus5 = _,_,_,_,_;
	// bus6 = _,_,_,_,_,_;
	// bus7 = _,_,_,_,_,_,_;
	// bus8 = _,_,_,_,_,_,_,_;
	// gain2(g) = *(g),*(g);
	// gain3(g) = *(g),*(g),*(g);
	// gain4(g) = *(g),*(g),*(g),*(g);
	// gain5(g) = *(g),*(g),*(g),*(g),*(g);
	// gain6(g) = *(g),*(g),*(g),*(g),*(g),*(g);
	// gain7(g) = *(g),*(g),*(g),*(g),*(g),*(g),*(g);
	// gain8(g) = *(g),*(g),*(g),*(g),*(g),*(g),*(g),*(g);