From 94bb45da7bbb0a1368758a2f4c1512c2b6ab4faf Mon Sep 17 00:00:00 2001 From: Alain Bonardi Date: Mon, 10 Oct 2022 15:45:50 +0200 Subject: [PATCH] Added routes and delays libs --- faustCodes/library/hoa2.lib | 881 ++++++++++++++++++------------------ 1 file changed, 442 insertions(+), 439 deletions(-) diff --git a/faustCodes/library/hoa2.lib b/faustCodes/library/hoa2.lib index 6b6c76fb..8f043d10 100755 --- a/faustCodes/library/hoa2.lib +++ b/faustCodes/library/hoa2.lib @@ -43,8 +43,11 @@ 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 "0.4"; declare author "Pierre Guillot"; declare author "Eliott Paris"; declare author "Julien Colafrancesco"; @@ -77,8 +80,7 @@ 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 +// 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 @@ -90,27 +92,27 @@ encoder(N, x, a) = encoder(N-1, x, a), x*sin(N*a), x*cos(N*a); // Where: // // * `N`: the ambisonic order (constant numerical expression) -// * `sp': the azimuth speed expressed as angular speed (2PI/sec), positive or negative +// * `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. +// * `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; - }; +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) +// Encoding of a stereo pair of channels with symetric angles (a/2, -a/2). // // #### Usage // @@ -121,18 +123,18 @@ rEncoder(N, sp, a, it) = thisEncoder // 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 +// * `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; +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 +// Encoding of a set of P signals distributed on the unit circle according to a list of P speeds and P angles. // // #### Usage // @@ -143,23 +145,23 @@ stereoEncoder(N, a) = (leftEncoder, rightEncoder) :> si.bus(2*N+1) // 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. +// * `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; - }; +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; + }; }; @@ -184,14 +186,14 @@ multiEncoder(N, lspeed, langle, it) = par(i, P, thisEncoder(ba.take(i+1, lspeed) //------------------------------------------------------------------- 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); + 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 +// (30° - 330°): Sound field lose energy around 180°. You should // use `inPhase` optimization with ponctual sources. // #### Usage // @@ -204,21 +206,21 @@ with { // * `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)) :> _) + (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)); - }; + 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)); + }; }; @@ -241,15 +243,15 @@ with { // * `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); +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 circularScaledVBAP function provides a circular scaled VBAP with all loudspeakers and the virtual source on the unit-circle +// The function provides a circular scaled VBAP with all loudspeakers and the virtual source on the unit-circle. // // #### Usage // @@ -259,46 +261,47 @@ iBasicDecoder(N, la, direct, shift) = (par(i, 2*N+1, _) <: par(i, P, speaker(N, // // 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 +// * `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))); +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 +// using 2D VBAP for compensation. // // #### Usage // @@ -314,17 +317,17 @@ circularScaledVBAP(l, t) = thisCircularVbap with { // * `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); +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 +// and to switch between multi-channel and stereo. // // #### Usage // @@ -335,27 +338,27 @@ imlsDecoder(N, la, direct, shift) = si.bus(2*N+1) : iVBAPDecoder // 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 +// * `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)); +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)); }; @@ -399,17 +402,17 @@ optimBasic(N) = par(i, 2*N+1, _); // * `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); - }; +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 +// The inPhase optimization optimizes energy vector and put all loudspeakers signals // in phase. It should be used for an auditory. // // #### Usage @@ -422,20 +425,20 @@ optimMaxRe(N) = par(i, 2*N+1, optim(i, N, _)) // // * `N`: the ambisonic order (constant numerical expression) //----------------------------------------------------- -optimInPhase(N) = par(i, 2*N+1, optim(i, N, _)) +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); - }; + 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, (ho.)optimInPhase +// `(ho).optimBasic`, `(ho).optimMaxRe` and `(ho.)optimInPhase`. // // #### Usage // @@ -446,14 +449,14 @@ with { // Where: // // * `N`: the ambisonic order (constant numerical expression) -// * 'ot' : optimization type (0 for optimBasic, 1 for optimMaxRe, 2 for optimInPhase) +// * `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); +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); }; @@ -473,24 +476,24 @@ optim(N, ot) = thisOptimizer // * `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, _)) +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); - }; - }; - }; + 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 +// Mirroring effect on the sound field. // // #### Usage // @@ -501,7 +504,7 @@ with { // Where: // // * `N`: the ambisonic order (constant numerical expression) -// * 'fa' : mirroring type (1=original sound field; 0=original+mirrored sound field; -1=mirrored sound field) +// * `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)))); @@ -523,10 +526,10 @@ mirror(N, fa) = (*(1), par(i, N, (*(fa), *(1)))); // * `r`: the radius // * `a`: the angle in radian //----------------------------------------------------- -map(N, x, r, a) = encoder(N, x * volume(r), a) : wider(N, ouverture(r)) +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); + volume(r) = 1. / (r * r * (r > 1) + (r <= 1)); + ouverture(r) = r * (r < 1) + (r >= 1); }; @@ -546,18 +549,18 @@ with { //----------------------------------------------------- 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); - }; + 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 +// Produces an XY pair of signals representing the ambisonic sound field. // // #### Usage // @@ -568,60 +571,61 @@ with { // Where: // // * `N`: the ambisonic order (constant numerical expression) -// * 'rt' : refreshment time in milliseconds +// * `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), _, _); +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 +//============================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 +// 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. +// `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) +// +// * 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) +//-(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 @@ -633,117 +637,120 @@ scope(N, rt) = thisScope // 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) +// * `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; +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 +// 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. +// `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, moire and more components are delayed. +// +// * 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) +//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) +//-(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) : _,_, ... +// _,_, ... : 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) +// * `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; +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 +// Spatial ring modulation in syn mode. // -//SYN RINGMOD GENERATES SPATIAL COMPONENTS IN AMBISONICS FROM ONE MONO SIGNAL THANKS TO RING MODULATION -// 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 abmisonic order n. +// `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. +// +// 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 +// +// 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 +//-(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 @@ -755,52 +762,54 @@ synDecorrelation(N, d, wf, fa, fd, tf) = _ <: par(i, 2*N+1, crossFade(d, i, 2*N // 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) +// * `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 ; +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 +// Spatial ring modulation in syn mode. // -// synRingMod generates spatial components in amibsonics from one mono signal thanks to ring modulation. -// The ring modulation is defined for each spatial component among P=2*n+1 at the abmisonic order n. +// `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. +// +// 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. +// +// 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 +//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 +//-(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 @@ -812,35 +821,32 @@ fxRingMod(N, f0, fa, tf) = par(i, 2*N+1, gate_ringmod(f0, i, 2*N+1, 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) +// * `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 ; +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); @@ -849,20 +855,20 @@ 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) ; +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); //======================================================================================== @@ -870,7 +876,6 @@ th(21, i, p) = sqrt(1 - ((i+1) / p - 1)^2) ; //======================================================================================== //======================================================================================== - //----------------------`(ho.)encoder3D`--------------------------------- // Ambisonic encoder. Encodes a signal in the circular harmonics domain // depending on an order of decomposition, an angle and an elevation. @@ -888,73 +893,71 @@ th(21, i, p) = sqrt(1 - ((i+1) / p - 1)^2) ; // * `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); - }; - }; - }; - }; + // 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 +// 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 @@ -966,24 +969,24 @@ with { // 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 +// * `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. +// * `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; - }; + 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; + }; }; @@ -1021,9 +1024,9 @@ optimBasic3D(N) = par(i, (N+1) * (N+1), _); //----------------------------------------------------- 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); + // The degree l of the harmonic[l, m] + degree(index) = int(sqrt(index)); + MaxRe(N, l, _)= _ * cos(l / (2*N+2) * ma.PI); }; @@ -1043,19 +1046,19 @@ with { //----------------------------------------------------- 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); - }; + // 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, (ho.)optimInPhase3D +// `(ho).optimBasic3D`, `(ho).optimMaxRe3D` and `(ho.)optimInPhase3D`. // // #### Usage // @@ -1066,13 +1069,13 @@ with { // Where: // // * `N`: the ambisonic order (constant numerical expression) -// * 'ot' : optimization type (0 for optimBasic, 1 for optimMaxRe, 2 for optimInPhase) +// * `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; +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; }; \ No newline at end of file