Various collection of (transpiled) FFT algorithms

 avatar
unknown
javascript
a month ago
59 kB
10
Indexable
/*
(15 Jan 2025 19:◯◯) I'm kinda obsessed with spectral sounds lately…
(2025-01-16 12:◯◯) There's so much I don't know, most of these are just small snippets converted (possibly even erroneously)…

[KVR] Forum index > Photosounder > Audio Resynthesis Algorithm (397392; 15 Jan:mmtgb, Ghost Archive:DskEq):
	The process involves resynthesizing sound using a bank of bandpass filters.
	These filters are spaced logarithmically across frequencies with varying widths to balance time and frequency resolution.
	The analysis uses a hybrid linear-logarithmic scale (p=a+b*exp(k*f) to convert frequency f to bin p; where k → 0 achieves linearity) for frequency bands.
	The transform is just a variant of the S-transform, except for using cosine-windows in place of Gaussian-windows.
	Equivalently, both forward transforms are instances of the wavelet transform in which the wavelets are complex exponentials enveloped with the given windows; i.e. psi(x) = g(x) exp(ikx) where psi(x) is the mother wavelet and g(x) is the windowing function.
	In fact, the whole reason for the S-transform was that the inverse transform for wavelets couldn't be done with Gaussian-windowed complex exponential wavelets.
[Wiki] constant-Q transform:
	In mathematics and signal processing, the constant-Q transform and variable-Q transform, simply known as CQT and VQT, transforms a data series to the frequency domain. It is related to the Fourier transform and very closely related to the complex Morlet wavelet transform. Its design is suited for musical representation.
	The transform can be thought of as a series of filters f_k, logarithmically spaced in frequency, with the k-th filter having a spectral width δf_k equal to a multiple of the previous filter's width.
*/

//■■■■ [GitHub] jnalon / fast-fourier-transform / c++ ■■■■ Converted to JavaScript ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
class Complex{
	constructor(re=0,im=0){this.r=re;this.i=im;}
	add(c){return new Complex(this.r+c.r,this.i+c.i);}
	sub(c){return new Complex(this.r-c.r,this.i-c.i);}
	mul(c){
		return new Complex(
			this.r*c.r-this.i*c.i,
			this.r*c.i+this.i*c.r
		);
	}
	mulScalar(a){return new Complex(this.r*a,this.i*a);}
	toString(){return`(${this.r}, ${this.i})`;}
}
function cexpn(a){return new Complex(Math.cos(a),Math.sin(a));}
function bitReverse(k,r){
	let l=0;
	for(let i=0;i<r;i++){
		l=(l<<1)+(k&1);
		k>>=1;
	}
	return l;
}
function factor(n){ //Smallest prime factor of a given number.
	const rn=Math.ceil(Math.sqrt(n));
	for(let i=2;i<=rn;i++)
		if(n%i===0)return i;
	return n;
}
function directFT(x){ //Discrete Fourier Transform (DFT) directly from the definition, an algorithm that has O(N²) complexity.
	const N=x.length, //There is no restriction on the size of the array.
	      X=new Array(N).fill(new Complex(0,0)),
	      W=cexpn(-2*Math.PI/N);
	let Wk=new Complex(1,0);
	for(let k=0;k<N;k++){
		X[k]=new Complex(0,0);
		let Wkn=new Complex(1,0);
		for(let n=0;n<N;n++){
			X[k]=X[k].add(Wkn.mul(x[n]));
			Wkn=Wkn.mul(Wk);
		}
		Wk=Wk.mul(W);
	}
	return X;
}
function recursiveFFT(x){ //Fast Fourier Transform (FFT) using a recursive decimation in time algorithm. This has O(N log_2(N)) complexity.
	const N=x.length; //This should always be called with an array of a power of two length, or it will fail.
	if(N===1){
		return x;
	}else{
		const N2=N>>1,
		      xe=new Array(N2),
		      xo=new Array(N2),
		      X=new Array(N);
		for(let k=0;k<N2;k++){
			xe[k]=x[k<<1];
			xo[k]=x[(k<<1)+1];
		}
		const Xe=recursiveFFT(xe),
					Xo=recursiveFFT(xo),
					W=cexpn(-2*Math.PI/N);
		let Wk=new Complex(1,0);
		for(let k=0;k<N2;k++){
			const WkXo=Wk.mul(Xo[k]);
			X[k]=Xe[k].add(WkXo);
			X[k+N2]=Xe[k].sub(WkXo);
			Wk=Wk.mul(W);
		}
		return X;
	}
}
function iterativeFFT(x){ //An iterative in-place decimation in time algorithm. This has O(N log_2(N)) complexity, and since there are less function calls, it will probably be marginally faster than the recursive versions.
	const N=x.length, //Should be a power of two.
	      r=Math.floor(Math.log2(N)),
	      X=new Array(N);
	for(let k=0;k<N;k++){
		const l=bitReverse(k,r);
		X[l]=x[k];
	}
	let step=1;
	for(let k=0;k<r;k++){
		for(let l=0;l<N;l+=2*step){
			const W=cexpn(-Math.PI/step);
			let Wkn=new Complex(1,0);
			for(let n=0;n<step;n++){
				const p=l+n,
							q=p+step,
							temp=Wkn.mul(X[q]);
				X[q]=X[p].sub(temp);
				X[p]=X[p].add(temp);
				Wkn=Wkn.mul(W);
			}
		}
		step<<=1;
	}
	return X;
}
function recursiveNFFT(x){ //Fast Fourier Transform using a recursive decimation in time algorithm. This has smaller complexity than the direct FT, though the exact value is difficult to compute.
	const N=x.length, //Its length must be a composite number, or else the computation will be defered to the direct FT, and there will be no efficiency gain.
	      N1=factor(N);
	if(N1===N){
		return directFT(x);
	}else{
		const N2=N/N1,
		      xj=new Array(N2),
		      X=new Array(N).fill(new Complex(0,0)),
		      W=cexpn(-2*Math.PI/N);
		let Wj=new Complex(1,0);
		for(let j=0;j<N1;j++){
			for(let n=0;n<N2;n++){
				xj[n]=x[n*N1+j];
			}
			const Xj=recursiveNFFT(xj);
			let Wkj=new Complex(1,0);
			for(let k=0;k<N;k++){
				X[k]=X[k].add(Xj[k%N2].mul(Wkj));
				Wkj=Wkj.mul(Wj);
			}
			Wj=Wj.mul(W);
		}
		return X;
	}
}

//Example usage: console.log(recursiveFFT(complexArray).map(c=>c.toString()).join(","));

//■■■■ [GitHub] jtfell / c-fft ■■■■ Converted to JavaScript ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
const complex=(re,im)=>({re,im}),
      convFromPolar=(r,radians)=>complex(r*Math.cos(radians),r*Math.sin(radians)),
      add=(left,right)=>complex(left.re+right.re,left.im+right.im),
      multiply=(left,right)=>complex(
      	left.re*right.re-left.im*right.im,
      	left.re*right.im+left.im*right.re
      ),
      PI=Math.PI;
const DFT_naive=(x,N)=>{ //DFT (Discrete Fourier Transform) using the naïve method.
	const X=new Array(N).fill(null).map(()=>complex(0,0));
	for(let k=0;k<N;k++){
		for(let n=0;n<N;n++){
			const twiddle=convFromPolar(1,-2*PI*n*k/N);
			X[k]=add(X[k],multiply(x[n],twiddle));
		}
	}
	return X;
};
const FFT_GoodThomas=(input,N,N1,N2)=>{ //Good–Thomas FFT algorithm.
	const columns=new Array(N1).fill(null).map(()=>new Array(N2).fill(null)),
	      rows   =new Array(N2).fill(null).map(()=>new Array(N1).fill(null));
	//Reshape input into N1 columns (using Good-Thomas indexing)
	for(let z=0;z<N;z++)
		columns[z%N1][z%N2]=input[z];
	//Compute N1 DFTs of length N2 using naïve method
	for(let k1=0;k1<N1;k1++)
		columns[k1]=DFT_naive(columns[k1],N2);
	//Transpose
	for(let k1=0;k1<N1;k1++){
		for(let k2=0;k2<N2;k2++){
			rows[k2][k1]=columns[k1][k2];
		}
	}
	//Compute N2 DFTs of length N1 using naïve method
	for(let k2=0;k2<N2;k2++)
		rows[k2]=DFT_naive(rows[k2],N1);
	//Flatten into single output (using Chinese Remainder Theorem)
	const output=new Array(N).fill(null);
	for(let k1=0;k1<N1;k1++){
		for(let k2=0;k2<N2;k2++){
			output[(N1*k2+N2*k1)%N]=rows[k2][k1];
		}
	}
	return output;
};
const FFT_CooleyTukey=(input,N,N1,N2)=>{ //Cooley–Tukey FFT algorithm.
	const columns=new Array(N1).fill(null).map(()=>new Array(N2).fill(null)),
	      rows   =new Array(N2).fill(null).map(()=>new Array(N1).fill(null));
	//Reshape input into N1 columns
	for(let k1=0;k1<N1;k1++){
		for(let k2=0;k2<N2;k2++){
			columns[k1][k2]=input[N1*k2+k1];
		}
	}
	//Compute N1 DFTs of length N2 using naïve method
	for(let k1=0;k1<N1;k1++)
		columns[k1]=DFT_naive(columns[k1],N2);
	//Multiply by the twiddle factors (e^(-2*pi*j/N * k1*k2)) and transpose
	for(let k1=0;k1<N1;k1++){
		for(let k2=0;k2<N2;k2++){
			const twiddle=convFromPolar(1,-2*PI*k1*k2/N);
			rows[k2][k1]=multiply(twiddle,columns[k1][k2]);
		}
	}
	//Compute N2 DFTs of length N1 using naïve method
	for(let k2=0;k2<N2;k2++)
		rows[k2]=DFT_naive(rows[k2],N1);
	//Flatten into single output
	const output=new Array(N).fill(null);
	for(let k1=0;k1<N1;k1++){
		for(let k2=0;k2<N2;k2++){
			output[N2*k1+k2]=rows[k2][k1];
		}
	}
	return output;
};

//■■■■ FFTW (Fastest Fourier Transform in the West) ■■■■ Converted from OCaml to JavaScript ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
//==== fftw-3.3.10\genfft\expr.ml =======================================================================================================================================================
/*****************************************
 * Symbolic Arithmetic Expressions
 *****************************************/
//Transcendent type
const Transcendent={
	I:"I",
	MULTI_A:"MULTI_A",
	MULTI_B:"MULTI_B",
	CONJ:"CONJ",
};
//Expression type
class Expr{
	static Num(n){return{type:"Num",value:n};}
	static NaN(t){return{type:"NaN",value:t};}
	static Plus(exprs){return{type:"Plus",exprs};}
	static Times(a,b){return{type:"Times",a,b};}
	static CTimes(a,b){return{type:"CTimes",a,b};}
	static CTimesJ(a,b){return{type:"CTimesJ",a,b};}
	static Uminus(expr){return{type:"Uminus",expr};}
	static Load(v){return{type:"Load",variable:v};}
	static Store(v,expr){return{type:"Store",variable:v,expr};}
}
//Assignment type
class Assignment{
	static Assign(v,expr){
		return{type:"Assign",variable:v,expr};
	}
}
/*****************************************
 * Hash Functions
 *****************************************/
//Hash function for floats
function hashFloat(x){
	const mantissa=x%1;
	const exponent=Math.floor(Math.log2(x));
	return Math.trunc(exponent*1234.567+mantissa*10000);
}
//Convert transcendent to float
function transcendentToFloat(t){
	switch(t){
		case Transcendent.I:return 2.718281828459045235360287471;
		case Transcendent.MULTI_A:return 0.6931471805599453094172321214;
		case Transcendent.MULTI_B:return-0.3665129205816643270124391582;
		case Transcendent.CONJ:return 0.6019072301972345747375400015;
		default:throw new Error("Unknown transcendent");
	}
}
//Hash function for expressions
function hash(expr){
	switch(expr.type){
		case"Num":return hashFloat(Number.toFloat(expr.value));
		case"NaN":return hashFloat(transcendentToFloat(expr.value));
		case"Load":return 1+1237*Variable.hash(expr.variable);
		case"Store":return 2*Variable.hash(expr.variable)-2345*hash(expr.expr);
		case"Plus":return 5+23451*expr.exprs.reduce((sum,e)=>sum+hash(e),0);
		case"Times":return 41+31415*(hash(expr.a)+hash(expr.b));
		case"CTimes":return 49+3245*(hash(expr.a)+hash(expr.b));
		case"CTimesJ":return 31+3471*(hash(expr.a)+hash(expr.b));
		case"Uminus":return 42+12345*hash(expr.expr);
		default:throw new Error("Unknown expression type");
	}
}
/*****************************************
 * Variable Handling
 *****************************************/
//Find all variables in an expression
function findVars(expr){
	switch(expr.type){
		case"Load":return[expr.variable];
		case"Plus":return expr.exprs.flatMap(findVars);
		case"Times":case"CTimes":case"CTimesJ":return[...findVars(expr.a),...findVars(expr.b)];
		case"Uminus":return findVars(expr.expr);
		default:return[];
	}
}
//Check if an expression is a constant
function isConstant(expr){
	switch(expr.type){
		case"Num":case"NaN":return true;
		case"Load":return Variable.isConstant(expr.variable);
		default:return false;
	}
}
//Check if an expression is a known constant
function isKnownConstant(expr){return expr.type==="Num"||expr.type==="NaN";}
/*****************************************
 * String Conversion
 *****************************************/
//Convert a list of strings to a concatenated string
function foldrStringConcat(list){return list.join(" ");}
//Convert transcendent to string
function stringOfTranscendent(t){return t;}
//Convert expression to string
function toString(expr,depth=10){
	if(depth===0)return"...";
	switch(expr.type){
		case"Load":return Variable.unparse(expr.variable);
		case"Num":return Number.toFloat(expr.value).toString();
		case"NaN":return stringOfTranscendent(expr.value);
		case"Plus":return`(+ ${foldrStringConcat(expr.exprs.map((e)=>toString(e,depth-1)))})`;
		case"Times":return`(* ${toString(expr.a,depth-1)} ${toString(expr.b,depth-1)})`;
		case"CTimes":return`(c* ${toString(expr.a,depth-1)} ${toString(expr.b,depth-1)})`;
		case"CTimesJ":return`(cj* ${toString(expr.a,depth-1)} ${toString(expr.b,depth-1)})`;
		case"Uminus":return`(- ${toString(expr.expr,depth-1)})`;
		case"Store":return`(:= ${Variable.unparse(expr.variable)} ${toString(expr.expr,depth-1)})`;
		default:throw new Error("Unknown expression type");
	}
}
//Convert assignment to string
function assignmentToString(assign){return`(:= ${Variable.unparse(assign.variable)} ${toString(assign.expr)})`;}
// Dump assignments to a print function
function dump(print,assignments){assignments.forEach((assign)=>print(`${assignmentToString(assign)}\n`));}
/*****************************************
 * Constants Handling
 *****************************************/
//Extract all constants from an expression
function exprToConstants(expr){
	switch(expr.type){
		case"Num":return[expr.value];
		case"Plus":return expr.exprs.flatMap(exprToConstants);
		case"Times":case"CTimes":case"CTimesJ":return [...exprToConstants(expr.a), ...exprToConstants(expr.b)];
		case"Uminus":return exprToConstants(expr.expr);
		default:return[];
	}
}
//Add a constant to a list if it doesn't already exist
function addFloatKeyValue(listSoFar,k){
	if(listSoFar.some((k2)=>Number.equal(k,k2))){
		return listSoFar;
	}else{
		return[k,...listSoFar];
	}
}
//Ensure constants in a list are unique
function uniqueConstants(constants){return constants.reduce(addFloatKeyValue,[]);}
//==== fftw-3.3.10\genfft\complex.ml ====================================================================================================================================================
/*****************************************
 * Complex Number Operations
 *****************************************/
class Complex{
	constructor(real,imag){
		this.real=real;
		this.imag=imag;
	}
	static zero=new Complex(Expr.makeNum(Number.zero),Expr.makeNum(Number.zero));
	static one=new Complex(Expr.makeNum(Number.one),Expr.makeNum(Number.zero));
	static i=new Complex(Expr.makeNum(Number.zero),Expr.makeNum(Number.one));
	static make(r,i){return new Complex(r,i);}
	static uminus(ce){return new Complex(Expr.makeUminus(ce.real),Expr.makeUminus(ce.imag));}
	static inverseInt(n){
		return new Complex(
			Expr.makeNum(Number.div(Number.one,Number.ofInt(n))),
			Expr.makeNum(Number.zero)
		);
	}
	static inverseIntSqrt(n){
		return new Complex(
			Expr.makeNum(Number.div(Number.one,Number.sqrt(Number.ofInt(n)))),
			Expr.makeNum(Number.zero)
		);
	}
	static intSqrt(n){
		return new Complex(
			Expr.makeNum(Number.sqrt(Number.ofInt(n))),
			Expr.makeNum(Number.zero)
		);
	}
	static nan(x){return new Complex(Expr.NaN(x),Expr.makeNum(Number.zero));}
	static half=Complex.inverseInt(2);
	static times3x3(ce1,ce2){
		const a=ce1.real,b=ce1.imag,c=ce2.real,d=ce2.imag;
		const realPart=Expr.makePlus([
			Expr.makeTimes(c,Expr.makePlus([a,Expr.makeUminus(b)])),
			Expr.makeTimes(b,Expr.makePlus([c,Expr.makeUminus(d)]))
		]);
		const imagPart=Expr.makePlus([
			Expr.makeTimes(a,Expr.makePlus([c,d])),
			Expr.makeUminus(Expr.makeTimes(c,Expr.makePlus([a,Expr.makeUminus(b)])))
		]);
		return new Complex(realPart,imagPart);
	}
	static times(ce1,ce2){
		if(!Magic.threemult){
			const a=ce1.real,b=ce1.imag,c=ce2.real,d=ce2.imag;
			const realPart=Expr.makePlus([
				Expr.makeTimes(a,c),
				Expr.makeUminus(Expr.makeTimes(b,d))
			]);
			const imagPart=Expr.makePlus([
				Expr.makeTimes(a,d),
				Expr.makeTimes(b,c)
			]);
			return new Complex(realPart,imagPart);
		}else if(Expr.isConstant(ce2.real)&&Expr.isConstant(ce2.imag)){
			return Complex.times3x3(ce1,ce2);
		}else{
			return Complex.times3x3(ce2,ce1);
		}
	}
	static ctimes(ce1,ce2){return new Complex(Expr.CTimes(ce1.real,ce2.real),Expr.makeNum(Number.zero));}
	static ctimesj(ce1,ce2){return new Complex(Expr.CTimesJ(ce1.real,ce2.real),Expr.makeNum(Number.zero));}
	static exp(n,i){
		const[c,s]=Number.cexp(n,i);
		return new Complex(Expr.makeNum(c),Expr.makeNum(s));
	}
	static sec(n,m){
		const[c]=Number.cexp(n,m);
		return new Complex(Expr.makeNum(Number.div(Number.one,c)),Expr.makeNum(Number.zero));
	}
	static csc(n,m){
		const[,s]=Number.cexp(n,m);
		return new Complex(Expr.makeNum(Number.div(Number.one,s)),Expr.makeNum(Number.zero));
	}
	static tan(n,m){
		const[c,s]=Number.cexp(n,m);
		return new Complex(Expr.makeNum(Number.div(s,c)),Expr.makeNum(Number.zero));
	}
	static cot(n,m){
		const[c,s]=Number.cexp(n,m);
		return new Complex(Expr.makeNum(Number.div(c,s)),Expr.makeNum(Number.zero));
	}
	static plus(ces){
		const[realParts,imagParts]=ces.reduce(
			([r,i],ce)=>[r.concat([ce.real]),i.concat([ce.imag])],
			[[],[]]
		);
		return new Complex(Expr.makePlus(realParts),Expr.makePlus(imagParts));
	}
	static real(ce){return new Complex(ce.real,Expr.makeNum(Number.zero));}
	static imag(ce){return new Complex(ce.imag,Expr.makeNum(Number.zero));}
	static iimag(ce){return new Complex(Expr.makeNum(Number.zero),ce.imag);}
	static conj(ce){return new Complex(ce.real,Expr.makeUminus(ce.imag));}
	static sigma(a,b,f){return Complex.plus(Util.interval(a,b).map(f));}
	static storeReal(v,ce){return Expr.Store(v,ce.real);}
	static storeImag(v,ce){return Expr.Store(v,ce.imag);}
	static store([vr,vi],ce){return[Complex.storeReal(vr,ce),Complex.storeImag(vi,ce)];}
	static assignReal(v,ce){return Expr.Assign(v,ce.real);}
	static assignImag(v,ce){return Expr.Assign(v,ce.imag);}
	static assign([vr,vi],ce){return[Complex.assignReal(vr,ce),Complex.assignImag(vi,ce)];}
}
/*****************************************
 * Shortcuts and Signal Handling
 *****************************************/
const times=Complex.times,
      plus=Complex.plus,
      uminus=(a,b)=>Complex.plus([a,Complex.uminus(b)]),
      mul=(a,b)=>Complex.times(a,b),
      add=(a,b)=>Complex.plus([a,b]),
      sub=(a,b)=>Complex.plus([a,Complex.uminus(b)]);
//Type of complex signals
class Signal{
	constructor(fn){
		this.fn=fn;
	}
	//Evaluate the signal at index i
	eval(i){
		return this.fn(i);
	}
	//Make a finite signal infinite
	static infinite(n,signal){
		return new Signal((i)=>{
			if(i>=0&&i<n)return signal.eval(i);
			return Complex.zero;
		});
	}
	//Hermitian symmetry for signals
	static hermitian(n,signal){
		return new Signal((i)=>{
			if(i===0)return Complex.real(signal.eval(0));
			if(i<n-i)return signal.eval(i);
			if(i>n-i)return Complex.conj(signal.eval(n-i));
			return Complex.real(signal.eval(i));
		});
	}
	//Anti-Hermitian symmetry for signals
	static antihermitian(n,signal){
		return new Signal((i)=>{
			if(i===0)return Complex.iimag(signal.eval(0));
			if(i<n-i)return signal.eval(i);
			if(i>n-i)return Complex.uminus(Complex.conj(signal.eval(n-i)));
			return Complex.iimag(signal.eval(i));
		});
	}
}
/*****************************************
 * Utility Functions
 *****************************************/
//Abstraction of sum_{i=0}^{n-1}
const sigma=(a,b,f)=>{
	const interval=Array.from({length:b-a},(_,i)=>a+i);
	return Complex.plus(interval.map(f));
};
//Store real and imaginary parts of a complex number
const storeReal=(v,ce)=>Expr.Store(v,ce.real),
      storeImag=(v,ce)=>Expr.Store(v,ce.imag),
      store=([vr,vi],ce)=>[storeReal(vr,ce),storeImag(vi,ce)];
//Assign real and imaginary parts of a complex number
const assignReal=(v,ce)=>Expr.Assign(v,ce.real),
      assignImag=(v,ce)=>Expr.Assign(v,ce.imag),
      assign=([vr,vi],ce)=>[assignReal(vr,ce),assignImag(vi,ce)];
//==== fftw-3.3.10\genfft\util.ml =======================================================================================================================================================
/*****************************************
 *Integer operations
 *****************************************/
//Find the inverse of n modulo m
function invmod(n,m){
	function loop(i){
		if((i*n)%m===1)return i;
		return loop(i+1);
	}
	return loop(1);
}
//Euclidean algorithm for GCD
function gcd(n,m){
	if(n>m)return gcd(m,n);
	const r=m%n;
	if(r===0)return n;
	return gcd(r,n);
}
//Reduce the fraction m/n to lowest terms,modulo factors of n/n
function lowestTerms(n,m){
	if(m%n===0)return[1,0];
	const nn=Math.abs(n);
	const mm=m*(n/nn);
	const mpos=mm>0?mm%nn:(mm+(1+Math.abs(mm)/nn)*nn)%nn;
	const d=gcd(nn,Math.abs(mm));
	return[nn/d,mpos/d];
}
//Find a generator for the multiplicative group mod p (p must be prime)
class NoGenerator extends Error{}
function findGenerator(p){
	function period(x,prod){
		if(prod===1)return 1;
		return 1+period(x,(prod*x)%p);
	}
	function findgen(x){
		if(x===0)throw new NoGenerator();
		if(period(x,x)===p-1)return x;
		return findgen((x+1)%p);
	}
	return findgen(1);
}
//Raise x to a power n modulo p (requires n>0)
class NegativePower extends Error{}
function powMod(x,n,p){
	if(n===0)return 1;
	if(n<0)throw new NegativePower();
	if(n%2===0)return powMod((x*x)%p,n/2,p);
	return(x*powMod(x,n-1,p))%p;
}
/******************************************
 *Auxiliary functions
 ******************************************/
//Functional forall (like a fold)
function forall(id,combiner,a,b,f){
	if(a>=b)return id;
	return combiner(f(a),forall(id,combiner,a+1,b,f));
}
//Sum of a list
function sumList(l){
	return l.reduce((acc,val)=>acc+val,0);
}
//Maximum of a list
function maxList(l){
	return l.reduce((acc,val)=>Math.max(acc,val),-Infinity);
}
//Minimum of a list
function minList(l){
	return l.reduce((acc,val)=>Math.min(acc,val),Infinity);
}
//Count elements satisfying a predicate
function count(pred,l){
	return l.reduce((acc,elem)=>pred(elem)?acc+1:acc,0);
}
//Remove an element from a list
function remove(elem,l){
	return l.filter(e=>e!==elem);
}
//Functional composition
function compose(f,g){
	return x=>f(g(x));
}
//Flatten a list
function forallFlat(a,b,f){
	return forall([],(a,b)=>a.concat(b),a,b,f);
}
//Identity function
function identity(x){
	return x;
}
//Find the minimum element in a list based on a function
function minimize(f,l){
	if(l.length===0)return null;
	const[elem,...rest]=l;
	const minRest=minimize(f,rest);
	if(minRest===null)return elem;
	return f(minRest)>=f(elem)?elem:minRest;
}
//Find an element satisfying a condition
function findElem(condition,l){
	for(const elem of l){
		if(condition(elem))return elem;
	}
	return null;
}
//Find x>=a such that pred(x) is true
function suchthat(a,pred){
	if(pred(a))return a;
	return suchthat(a+1,pred);
}
//Print an information message
function info(string){
	if(Magic.verbose){
		const now=Date.now();
		const pid=process.pid;
		console.error(`${pid}:at t=${now}:${string}`);
	}
}
//Generate a list [0,1,...,n-1]
function iota(n){
	return forall([],(a,b)=>a.concat([b]),0,n,identity);
}
//Generate a list [a,a+1,...,b-1]
function interval(a,b){
	return iota(b-a).map(x=>x+a);
}
//Freeze a function into an array (lazy evaluation)
function array(n,f){
	const a=Array.from({length:n},(_,i)=>()=>f(i));
	return i=>a[i]();
}
//Take the first n elements of a list
function take(n,l){
	if(n===0)return[];
	if(l.length===0)throw new Error("take");
	const[a,...b]=l;
	return[a,...take(n-1,b)];
}
//Drop the first n elements of a list
function drop(n,l){
	if(n===0)return l;
	if(l.length===0)throw new Error("drop");
	const[,...b]=l;
	return drop(n-1,b);
}
//Return the first non-null value
function either(a,b){
	return a!==null&&a!==undefined?a:b;
}
//==== fftw-3.3.10\genfft\fft.ml ========================================================================================================================================================
class Complex{
	constructor(real,imag){
		this.real=real;
		this.imag=imag;
	}
	add(other){
		return new Complex(this.real+other.real,this.imag+other.imag);
	}
	sub(other){
		return new Complex(this.real-other.real,this.imag-other.imag);
	}
	mul(other){
		return new Complex(
			this.real*other.real-this.imag*other.imag,
			this.real*other.imag+this.imag*other.real
		);
	}
	conj(){
		return new Complex(this.real,-this.imag);
	}
	static zero=new Complex(0,0);
	static one=new Complex(1,0);
	static i=new Complex(0,1);
}
function gcd(a,b){
	while(b!==0){
		let temp=b;
		b=a%b;
		a=temp;
	}
	return a;
}
function chooseFactor(n){
	function choose1(n){
		let f=1;
		for(let i=1;i*i<=n;i++){
			if(n%i===0&&gcd(i,n/i)===1){
				f=i;
			}
		}
		return f;
	}
	function choose2(n){
		let f=1;
		for(let i=1;i*i<=n;i++){
			if(n%i===0){
				f=i;
			}
		}
		return f;
	}
	const i=choose1(n);
	return i>1?i:choose2(n);
}
function isPowerOfTwo(n){
	return n>0&&(n&(n-1))===0;
}
function dftPrime(sign,n,input){
	function sum(filter,i){
		let result=Complex.zero;
		for(let j=0;j<n;j++){
			const coeff=filter(exp(n,sign*i*j));
			result=result.add(coeff.mul(input(j)));
		}
		return result;
	}
	const computationEven=Array.from({length:n},(_,i)=>sum(x=>x,i)),
	      computationOdd=(()=>{
		const sumr=Array.from({length:n},(_,i)=>sum(x=>x.real,i)),
		      sumi=Array.from({length:n},(_,i)=>sum(x=>x.imag.mul(Complex.i),i));
		return Array.from({length:n},(i)=>{
			if(i===0){
				return input(0).add(Array.from({length:(n+1)/2},(j)=>input(j).add(input(n-j))).reduce((a,b)=>a.add(b),Complex.zero));
			}else{
				const iPrime=Math.min(i,n-i);
				return i<n-i?sumr[iPrime].add(sumi[iPrime]):sumr[iPrime].sub(sumi[iPrime]);
			}
		});
	})();
	if(n>=Magic.raderMin){
		return dftRader(sign,n,input);
	}else if(n===2){
		return computationEven;
	}else{
		return computationOdd;
	}
}
function dftRader(sign,p,input){
	const half=new Complex(0.5,0);
	function makeProduct(n,a,b){
		const scaleFactor=inverseInt(n);
		return Array.from({length:n},(i)=>a(i).mul(b(i).mul(scaleFactor)));
	}
	function genConvolutionByFFT(n,a,b,addToAll){
		const fftA=dft(1,n,a),
		      fftB=dft(1,n,b),
		      fftAB=makeProduct(n,fftA,fftB),
		      dcTerm=i=>i===0?addToAll:Complex.zero,
		      fftAB1=Array.from({length:n},(i)=>fftAB(i).add(dcTerm(i))),
		      sum=fftA(0),
		      conv=dft(-1,n,fftAB1);
		return[sum,conv];
	}
	function genConvolutionByFFTAlt(n,a,b,addToAll){
		const ap=Array.from({length:n},(i)=>half.mul(a(i).add(a((n-i)%n)))),
		      am=Array.from({length:n},(i)=>half.mul(a(i).sub(a((n-i)%n)))),
		      bp=Array.from({length:n},(i)=>half.mul(b(i).add(b((n-i)%n)))),
		      bm=Array.from({length:n},(i)=>half.mul(b(i).sub(b((n-i)%n)))),
		      fftAp=dft(1,n,ap),
		      fftAm=dft(1,n,am),
		      fftBp=dft(1,n,bp),
		      fftBm=dft(1,n,bm),
		      fftABpp=makeProduct(n,fftAp,fftBp),
		      fftABpm=makeProduct(n,fftAp,fftBm),
		      fftABmp=makeProduct(n,fftAm,fftBp),
		      fftABmm=makeProduct(n,fftAm,fftBm),
		      sum=fftAp(0).add(fftAm(0)),
		      dcTerm=i=>i===0?addToAll:Complex.zero,
		      fftAB1=Array.from({length:n},(i)=>fftABpp(i).add(fftABmm(i)).add(dcTerm(i))),
		      fftAB2=Array.from({length:n},(i)=>fftABpm(i).add(fftABmp(i))),
		      conv1=dft(-1,n,fftAB1),
		      conv2=dft(-1,n,fftAB2),
		      conv=Array.from({length:n},(i)=>conv1(i).add(conv2(i)));
		return[sum,conv];
	}
	const genConvolution=p<=Magic.alternateConvolution?genConvolutionByFFTAlt:genConvolutionByFFT,
	      g=findGenerator(p),
	      gInv=powMod(g,p-2,p),
	      inputPerm=Array.from({length:p},(i)=>input(powMod(g,i,p))),
	      omegaPerm=Array.from({length:p},(i)=>exp(p,sign*powMod(gInv,i,p))),
	      outputPerm=Array.from({length:p},(i)=>powMod(gInv,i,p)),
	      [sum,conv]=genConvolution(p-1,inputPerm,omegaPerm,input(0));
	return Array.from({length:p},(i)=>{
		if(i===0){
			return input(0).add(sum);
		}else{
			const iPrime=outputPerm.findIndex((val)=>val===i);
			return conv(iPrime);
		}
	});
}
function newsplit(sign,n,input){
	function s(n,k){
		if(n<=4)return Complex.one;
		const k4=Math.abs(k)%(n/4),
		      k4Prime=k4<=n/8?k4:n/4-k4;
		return s(n/4,k4Prime).mul(new Complex(Math.cos((2*Math.PI*k4Prime)/n),0));
	}
	function sinv(n,k){
		if(n<=4)return Complex.one;
		const k4=Math.abs(k)%(n/4),
		      k4Prime=k4<=n/8?k4:n/4-k4;
		return sinv(n/4,k4Prime).mul(new Complex(1/Math.cos((2*Math.PI*k4Prime)/n),0));
	}
	function sdiv2(n,k){
		return s(n,k).mul(sinv(2*n,k));
	}
	function sdiv4(n,k){
		const k4=Math.abs(k)%n,
		      angle=k4<=n/2?k4:n-k4;
		return new Complex(1/Math.cos((2*Math.PI*angle)/(4*n)),0);
	}
	function t(n,k){
		return exp(n,k).mul(sdiv4(n/4,k));
	}
	function dft1(input){
		return input;
	}
	function dft2(input){
		return Array.from({length:2},(k)=>input(0).add(input(1).mul(exp(2,k))));
	}
	function newsplit0(sign,n,input){
		if(n===1)return dft1(input);
		if(n===2)return dft2(input);
		const u=newsplit0(sign,n/2,(i)=>input(i*2)),
		      z=newsplitS(sign,n/4,(i)=>input(i*4+1)),
		      zPrime=newsplitS(sign,n/4,(i)=>input((n+i*4-1)%n)),
		      twid=Array.from({length:n},(k)=>s(n/4,k).mul(exp(n,sign*k))),
		      w=Array.from({length:n},(k)=>twid[k].mul(z(k%(n/4)))),
		      wPrime=Array.from({length:n},(k)=>twid[k].conj().mul(zPrime(k%(n/4)))),
		      ww=Array.from({length:n},(k)=>w[k].add(wPrime[k]));
		return Array.from({length:n},(k)=>u(k%(n/2)).add(ww[k]));
	}
	function newsplitS(sign,n,input){
		if(n===1)return dft1(input);
		if(n===2)return dft2(input);
		const u=newsplitS2(sign,n/2,(i)=>input(i*2)),
		      z=newsplitS(sign,n/4,(i)=>input(i*4+1)),
		      zPrime=newsplitS(sign,n/4,(i)=>input((n+i*4-1)%n)),
		      w=Array.from({length:n},(k)=>t(n,sign*k).mul(z(k%(n/4)))),
		      wPrime=Array.from({length:n},(k)=>t(n,sign*k).conj().mul(zPrime(k%(n/4)))),
		      ww=Array.from({length:n},(k)=>w[k].add(wPrime[k]));
		return Array.from({length:n},(k)=>u(k%(n/2)).add(ww[k]));
	}
	function newsplitS2(sign,n,input){
		if(n===1)return dft1(input);
		if(n===2)return dft2(input);
		const u=newsplitS4(sign,n/2,(i)=>input(i*2)),
		      z=newsplitS(sign,n/4,(i)=>input(i*4+1)),
		      zPrime=newsplitS(sign,n/4,(i)=>input((n+i*4-1)%n)),
		      w=Array.from({length:n},(k)=>t(n,sign*k).mul(z(k%(n/4)))),
		      wPrime=Array.from({length:n},(k)=>t(n,sign*k).conj().mul(zPrime(k%(n/4)))),
		      ww=Array.from({length:n},(k)=>w[k].add(wPrime[k]).mul(sdiv2(n,k)));
		return Array.from({length:n},(k)=>u(k%(n/2)).add(ww[k]));
	}
	function newsplitS4(sign,n,input){
		if(n===1)return dft1(input);
		if(n===2){
			const f=dft2(input);
			return Array.from({length:2},(k)=>f(k).mul(sinv(8,k)));
		}
		const u=newsplitS2(sign,n/2,(i)=>input(i*2)),
		      z=newsplitS(sign,n/4,(i)=>input(i*4+1)),
		      zPrime=newsplitS(sign,n/4,(i)=>input((n+i*4-1)%n)),
		      w=Array.from({length:n},(k)=>t(n,sign*k).mul(z(k%(n/4)))),
		      wPrime=Array.from({length:n},(k)=>t(n,sign*k).conj().mul(zPrime(k%(n/4)))),
		      ww=Array.from({length:n},(k)=>w[k].add(wPrime[k]));
		return Array.from({length:n},(k)=>u(k%(n/2)).add(ww[k]).mul(sdiv4(n,k)));
	}
	return newsplit0(sign,n,input);
}
function dft(sign,n,input){
	function cooleyTukey(sign,n1,n2,input){
		const tmp1=Array.from({length:n2},(i2)=>dft(sign,n1,(i1)=>input(i1*n2+i2))),
		      tmp2=Array.from({length:n1},(i1)=>
			Array.from({length:n2},(i2)=>exp(n,sign*i1*i2).mul(tmp1[i2][i1]))
		),
		      tmp3=Array.from({length:n1},(i1)=>dft(sign,n2,tmp2[i1]));
		return(i)=>tmp3[i%n1][Math.floor(i/n1)];
	}
	function splitRadixDit(sign,n,input){
		const f0=dft(sign,n/2,(i)=>input(i*2)),
		      f10=dft(sign,n/4,(i)=>input(i*4+1)),
		      f11=dft(sign,n/4,(i)=>input((n+i*4-1)%n)),
		      g10=Array.from({length:n},(k)=>exp(n,sign*k).mul(f10(k%(n/4)))),
		      g11=Array.from({length:n},(k)=>exp(n,-sign*k).mul(f11(k%(n/4)))),
		      g1=Array.from({length:n},(k)=>g10[k].add(g11[k]));
		return Array.from({length:n},(k)=>f0(k%(n/2)).add(g1[k]));
	}
	function splitRadixDif(sign,n,input){
		const n2=n/2,
		      n4=n/4,
		      x0=Array.from({length:n2},(i)=>input(i).add(input(i+n2))),
		      x10=Array.from({length:n4},(i)=>input(i).sub(input(i+n2))),
		      x11=Array.from({length:n4},(i)=>input(i+n4).sub(input(i+n2+n4))),
		      x1=(k,i)=>exp(n,k*i*sign).mul(x10[i].add(exp(4,k*sign).mul(x11[i]))),
		      f0=dft(sign,n2,x0),
		      f1=Array.from({length:4},(k)=>dft(sign,n4,(i)=>x1(k,i)));
		return Array.from({length:n},(k)=>{
			if(k%2===0){
				return f0(k/2);
			}else{
				const kPrime=k%4;
				return f1[kPrime][(k-kPrime)/4];
			}
		});
	}
	function primeFactor(sign,n1,n2,input){
		const tmp1=Array.from({length:n2},(i2)=>
			dft(sign,n1,(i1)=>input((i1*n2+i2*n1)%n))
		);
		const tmp2=Array.from({length:n1},(i1)=>dft(sign,n2,(k2)=>tmp1[k2][i1]));
		return(i)=>tmp2[i%n1][i%n2];
	}
	function algorithm(sign,n){
		const r=chooseFactor(n);
		if(Magic.raderList.includes(n)){
			return dftRader(sign,n);
		}else if(r===1){
			return dftPrime(sign,n);
		}else if(gcd(r,n/r)===1){
			return primeFactor(sign,r,n/r);
		}else if(n%4===0&&n>4){
			if(Magic.newsplit&&isPowerOfTwo(n)){
				return newsplit(sign,n);
			}else if(Magic.difSplitRadix){
				return splitRadixDif(sign,n);
			}else{
				return splitRadixDit(sign,n);
			}
		}else{
			return cooleyTukey(sign,r,n/r);
		}
	}
	return Array.from({length:n},(i)=>algorithm(sign,n)(i));
}

//■■■■ [GitHub] Photosounder / rouziclib / rouziclib / math / dsp.c ■■■■ Converted to JavaScript ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
//Helper functions…
function sq(x){return x*x;}
function rangelimit(x,min,max){return Math.min(Math.max(x,min),max);}
function erf(x){ //Approximation of the error function
	const a1= 0.254829592,
	      a2=-0.284496736,
	      a3= 1.421413741,
	      a4=-1.453152027,
	      a5= 1.061405429,
	      p = 0.3275911,
	      sign=x<0?-1:1;
	x=Math.abs(x);
	const t=1/(1+p*x),
				y=1-(((((a5*t+a4)*t)+a3)*t+a2)*t+a1)*t*Math.exp(-x*x);
	return sign*y;
}
//Main functions…
function is_composite(x){
	const p=[2,3,5];
	for(let i=0;i<p.length;i++){
		while(x%p[i]===0)x/=p[i];
	}
	return x;
}
function next_composite(x){
	while(is_composite(x)!==1)
		x++;
	return x;
}
function next_fast_fft_size(x){
	const fft_size=[
		60,64,72,75,80,100,120,125,128,135,144,150,160,180,200,225,240,256,300,320,375,400,405,432,500,
		512,540,576,625,640,648,720,768,800,810,864,900,960,1000,1080,1152,1200,1280,1296,1350,1440,1500,
		1600,1620,1728,1800,1920,2000,2025,2160,2304,2400,2560,2700,2880,3000,3200,3375,3600,3840,4000,4050,
		4320,4500,4800,5120,5184,5400,5625,6000,6400,6480,6912,7200,7500,8000,8100,8640,8748,9000,9375,9600,
		10000,10368,10800,10935,11520,12500,12960,13500,13824,14400,14580,15000,15625,16000,16200,17280,18000,
		19200,19440,20000,20736,21600,23040,24000,24300,25920,27000,28800,30000,30375,32000,32400,34560,36000,
		38400,40000,40500,43200,45000,46080,48000,50000,50625,51840,54000,57600,60000,62500,64000,67500,72000,
		72900,75000,80000,84375,90000,96000,97200,100000,103680,104976,112500,120000,125000,129600,131220,140625,
		150000,160000,162000,164025,174960,187500,200000,202500,207360,216000,218700,234375,250000,259200,270000,
		273375,291600,312500,324000,337500,345600,360000,364500,390625,405000,421875,432000,450000,455625,486000,
		518400,540000,562500,576000,607500,648000,675000,703125,720000,759375,781250,800000,810000,843750,864000,
		900000,937500,1000000,1012500,1080000,1125000,1171875,1200000,1250000,1265625,1280000,1350000,1406250,
		1440000,1500000,1562500,1600000,1687500,1875000,2000000,2109375,2250000,2343750,2400000,2500000,2531250,
		2812500,3125000,3200000,3515625,3750000,3906250,4000000,4100625,4687500,5000000,5062500,5184000,5400000,
		5467500,5859375,6250000,6328125,6480000,6561000,6635520],
	     count=fft_size.length;
	for(let i=0;i<count;i++){
		if(fft_size[i]>=x)return fft_size[i];
	}
	return next_composite(x);
}
function calc_2D_fast_convolution_dim(a,b){
	return func1_xyi(sub_xyi(add_xyi(a,b),set_xyi(1)),next_fast_fft_size);
}
function array_sum(s){
	let sum=0;const n=s.length;
	for(let i=0;i<n;i++)sum+=s[i];
	return sum;
}
function array_sum_sq(s){
	let sum=0;const n=s.length;
	for(let i=0;i<n;i++)sum+=s[i]*s[i];
	return sum;
}
function root_mean_square(s){
	return Math.sqrt(array_sum_sq(s)/s.length);
}
function root_mean_squaref_chan(s,n,ic,chan_count){
	let sum=0;
	for(let i=ic;i<n*chan_count;i+=chan_count)
		sum+=s[i]*s[i];
	return Math.sqrt(sum/n);
}
function db_to_vol(db){return Math.pow(10,db/20);}
function vol_to_db(vol){return 20*Math.log10(vol);}
function sinc(x,fc){
	if(x===0)return 1;
	const a=2*Math.PI*x*fc;
	return Math.sin(a)/a;
}
function blackman(x,range){
	x=Math.abs(x/range);
	if(x>=1)return 0;
	x*=Math.PI;
	return.42+.5*Math.cos(x)+.08*Math.cos(2*x);
}
function gaussian(x){return Math.exp(-x*x);}
function short_gaussian_window(x,range,w){
	x=Math.abs(x/range);
	if(x>=1)return 0;
	w*=.7071067811865475244;
	if(w<1e-4)return sq(1-sq(x));
	return sq(gaussian(x*w)-gaussian(w))/sq(1-gaussian(w));
}
function short_erf(x,w){
	x=rangelimit(x,-w,w)*Math.sqrt(.5);
	w*=Math.sqrt(.5);
	return(Math.sqrt(Math.PI)*Math.exp(w*w)*(Math.sqrt(2)*Math.exp(w*w)*erf(Math.sqrt(2)*x)-4*erf(x))+4*x)/
				(Math.sqrt(Math.PI)*Math.exp(w*w)*(Math.sqrt(2)*Math.exp(w*w)*erf(Math.sqrt(2)*w)-4*erf(w))+4*w);
}
function cumulative_squared_parabola(x){
	x=rangelimit(x,-1,1);
	const x2=x*x;
	return((.375*x2-1.25)*x2+1.875)*x;
}
function cumulative_cubed_parabola(x){
	x=rangelimit(x,-1,1);
	const x2=x*x;
	return(((-.3125*x2+1.3125)*x2-2.1875)*x2+2.1875)*x;
}
function ramp_kernel(x){
	if(x===0)return.25;
	return.5*Math.sin(Math.PI*x)/(Math.PI*x)-.25*sq(Math.sin(Math.PI*x*.5)/(Math.PI*x*.5));
}
function ramp_kernel_discrete(x){
	if(x===0)return.25;
	if(x&1)return-1/sq(x*Math.PI);
	return 0;
}
//Helper functions for 2D convolution (assuming xyi_t is a 2D point object)…
function set_xyi(value){return{x:value,y:value};}
function add_xyi(a,b){return{x:a.x+b.x,y:a.y+b.y};}
function sub_xyi(a,b){return{x:a.x-b.x,y:a.y-b.y};}
function func1_xyi(point,func){return{x:func(point.x),y:func(point.y)};}

//■■■■ The ARSS (Analysis & Resynthesis Sound Spectrograph) by Michel Rouzic ■■■■ arss.sourceforge.net/code ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
//==== util.c ===========================================================================================================================================================================
function roundoff(double x){ //nearbyint() replacement, with the exception that the result contains a non-zero fractional part
	return x>0?x+.5:x-.5;
}
const fmod=(n,m)=>(n%m+m)%m,
      roundup=Math.ceil;
function smallprimes(x){ //Returns 1 if x is only made of small primes
	const p=[2,3];
	for(let i=0;i<2;i++){
		while(x%p[i]===0)x/=p[i];
	}
	return x===1?1:x;
}
function nextsprime(x){ //Returns the next integer only made of small primes
	while(smallprimes(x)!==1)x++;
	return x;
}
let LOGBASE=2; //Example base, adjust as needed…
function log_b(x){ //Returns the logarithm of x with a given base (LOGBASE)
	if(LOGBASE===1){
		return x;
	}else{
		if(x===0)
			console.warn("Warning: log(0) returning -Infinity");
		return Math.log(x)/Math.log(LOGBASE);
	}
}
let RAND_MAX=2147483647; //Example RAND_MAX, adjust as needed…
function rand_u32(){ //Returns a random 32-bit unsigned integer
	if(RAND_MAX===2147483647){
		return Math.floor(Math.random()*(RAND_MAX+1));
	}else if(RAND_MAX===32767){
		return((Math.floor(Math.random()*256)<<24)|
		       (Math.floor(Math.random()*256)<<16)|
		       (Math.floor(Math.random()*256)<<8)|
		       (Math.floor(Math.random()*256)));
	}else{
		console.error(`Unhandled RAND_MAX value: ${RAND_MAX}\nPlease signal this error to the developer.`);
		return Math.floor(Math.random()*(RAND_MAX+1));
	}
}
function dblrand(){ //Range is +/- 1.0
	return(rand_u32()*(1/2147483648))-1;
}
//==== dsp.c ============================================================================================================================================================================
const PI_D=3.1415926535897932,
      LOGBASE_D=2,
      LOOP_SIZE_SEC_D=10,
      BMSQ_LUT_SIZE_D=16000,
      TRANSITION_BW_SYNT=16; //Defines the transition bandwidth for the low-pass filter on the envelopes during synthesisation

function fft(input,out,N,method){
	/* method :
	 * 0 = DFT
	 * 1 = IDFT
	 * 2 = DHT
	 */
	const p=fftw_plan_r2r_1d(N,input,out,method,FFTW_ESTIMATE); //fftw_plan_r2r_1d is a basic FFTW interface call for planning a 1-d real-to-real FFT operation on a single, contiguous data sequence.
	fftw_execute(p);
	fftw_destroy_plan(p);
}
function normi(s,xs,ys,ratio){
	let max=0,
	    maxx=0,maxy=0;
	for(let iy=0;iy<ys;iy++){
		for(let ix=0;ix<xs;ix++){
			if(Math.abs(s[iy][ix])>max){
				max=Math.abs(s[iy][ix]);
				maxx=ix;
				maxy=iy;
			}
		}
	}
	if(DEBUG)
		console.log(`norm:${max.toFixed(3)}(Y:${maxy}X:${maxx})`);
	if(max!==0){
		max/=ratio;
		max=1/max;
	}else{
		max=0;
	}
	for(let iy=0;iy<ys;iy++){
		for(let ix=0;ix<xs;ix++){
			s[iy][ix]*=max;
		}
	}
	if(DEBUG){
		console.log(`ex:${s[0][0].toFixed(3)}`);
	}
}
function log_pos(x,min,max){
	if(LOGBASE===1){
		return x*(max-min)+min;
	}else{
		return(max-min)*(min*Math.pow(LOGBASE,x*(Math.log(max)-Math.log(min))/Math.log(2.0))-min)/(min*Math.pow(LOGBASE,(Math.log(max)-Math.log(min))/Math.log(2.0))-min)+min;
	}
}
function log_pos_inv(x,min,max){
	if(LOGBASE===1.0){
		return(x-min)/(max-min);
	}else{
		return Math.log(((x-min)*(min*Math.pow(LOGBASE,(Math.log(max)-Math.log(min))/Math.log(2.0))-min)/(max-min)+min)/Math.log(LOGBASE))*Math.log(2.0)/(Math.log(max)-Math.log(min));
	}
}
function freqarray(basefreq,bands,bandsperoctave){
	let freq=new Array(bands),
	    maxfreq;
	if(LOGBASE===1.0){
		maxfreq=bandsperoctave;
	}else{
		maxfreq=basefreq*Math.pow(LOGBASE,(bands-1)/bandsperoctave);
	}
	for(let i=0;i<bands;i++){
		freq[i]=log_pos(i/(bands-1),basefreq,maxfreq);
	}
	if(log_pos(bands/(bands-1),basefreq,maxfreq)>0.5){
		console.log("Warning:Upper frequency limit above Nyquist frequency");
	}
	return freq;
}
function blackman_downsampling(inArray,Mi,Mo){
	let out=new Array(Mo).fill(0),
	    ratio=Mi/Mo,
	    ratio_i=1.0/ratio;
	for(let i=0;i<Mo;i++){
		let pos_in=i*ratio,
		    coef_sum=0;
		for(let j=Math.round(pos_in-ratio);j<=pos_in+ratio;j++){
			if(j>=0&&j<Mi){
				let x=j-pos_in+ratio;
				let coef=0.42-0.5*Math.cos(pi*x*ratio_i)+0.08*Math.cos(2*pi*x*ratio_i);
				coef_sum+=coef;
				out[i]+=inArray[j]*coef;
			}
		}
		out[i]/=coef_sum;
	}
	return out;
}
function bmsq_lut(size){
	let i, //general purpose iterator
	    coef, //Blackman square final coefficient
	    foo;
	const bar=Math.PI*(3/size)*(1/1.5);
	const f1=-0.6595044010905501, //Blackman Square coefficients
	      f2= 0.1601741366715479,
	      f4=-0.0010709178680006,
	      f5= 0.0001450093579222,
	      f7= 0.0001008528049040,
	      f8= 0.0000653092892874,
	      f10=0.0000293385615146,
	      f11=0.0000205351559060,
	      f13=0.0000108567682890,
	      f14=0.0000081549460136,
	      f16=0.0000048519309366,
	      f17=0.0000038284344102,
	      f19=0.0000024753630724;
	size++; //allows to read value 3.0
	const lut=new Array(size).fill(0);
	for(i=0;i<size;i++){
		foo=i*bar;
		coef=0;
		coef+=Math.cos(foo)*f1-f1;
		coef+=Math.cos(2*foo)*f2-f2;
		coef+=Math.cos(4*foo)*f4-f4;
		coef+=Math.cos(5*foo)*f5-f5;
		coef+=Math.cos(7*foo)*f7-f7;
		coef+=Math.cos(8*foo)*f8-f8;
		coef+=Math.cos(10*foo)*f10-f10;
		coef+=Math.cos(11*foo)*f11-f11;
		coef+=Math.cos(13*foo)*f13-f13;
		coef+=Math.cos(14*foo)*f14-f14;
		coef+=Math.cos(16*foo)*f16-f16;
		coef+=Math.cos(17*foo)*f17-f17;
		coef+=Math.cos(19*foo)*f19-f19;
		lut[i]=coef;
	}
	return lut;
}
function blackman_square_interpolation(inArray,outArray,Mi,Mo,lut,lutSize){
	let i,j, //general purpose iterators
	    posIn, //position in the original signal
	    x, //position of the iterator in the blackman_square(x)formula
	    ratio, //scaling ratio(>1.0)
	    ratioI, //ratio^-1
	    coef, //Blackman square final coefficient
	    posLut, //Index on the look-up table
	    posLuti, //Integer index on the look-up table
	    modPos, //modulo of the position on the look-up table
	    y0,y1; //values of the two closest values on the LUT
	const foo=lutSize/3;
	let jStart,jStop; //boundary values for the j loop
	/*
	 * Mi is the original signal's length
	 * Mo is the output signal's length
	 */
	ratio=Mi/Mo;
	ratioI=1/ratio;
	for(i=0;i<Mo;i++){
		posIn=i*ratio;
		jStop=Math.floor(posIn+1.5);
		jStart=jStop-2;
		if(jStart<0)jStart=0;
		if(jStop>=Mi) //The boundary check is done after jStart is calculated to avoid miscalculating it
			jStop=Mi-1;
		for(j=jStart;j<=jStop;j++){
			x=j-posIn+1.5; //calculate position within the Blackman square function in the[0.0;3.0]range
			posLut=x*foo;
			posLuti=Math.floor(posLut);
			modPos=posLut%1; //modulo of the index
			y0=lut[posLuti]; //interpolate linearly between the two closest values
			y1=lut[posLuti+1];
			coef=y0+modPos*(y1-y0); //linear interpolation
			outArray[i]+=inArray[j]*coef; //convolve
		}
	}
}
function anal(s,sampleCount,sampleRate,Xsize,bands,bpo,pixpersec,basefreq){
	let i,ib,Mb,Mc,Md,Fa,Fd,
	    out,h,freq,t,coef,La,Ld,Li,maxfreq;
	freq=freqarray(basefreq,bands,bpo);
	if(LOGBASE===1)
		maxfreq=bpo;
	else
		maxfreq=basefreq*Math.pow(LOGBASE,(bands-1)/bpo);
	Xsize=sampleCount*pixpersec;
	if(fmod(sampleCount*pixpersec,1)!==0)
		Xsize++;
	console.log(`Image size : ${Xsize}×${bands}`);
	out=new Array(bands);
	clocka=new Date().getTime();
	if(LOGBASE===1)
		Mb=sampleCount-1+Math.round(5/(freq[1]-freq[0]));
	else
		Mb=sampleCount-1+Math.round(2*5/((freq[0]*Math.pow(LOGBASE,-1/bpo))*(1-Math.pow(LOGBASE,-1/bpo))));
	if(Mb%2===1)
		Mb++;
	Mb=Math.round(nextsprime(Math.round(Mb*pixpersec))/pixpersec);
	Md=Math.round(Mb*pixpersec);
	s=new Float64Array(Mb);
	s.fill(0,sampleCount);
	fft(s,s,Mb,0);
	for(ib=0;ib<bands;ib++){
		Fa=Math.round(log_pos((ib-1)/(bands-1),basefreq,maxfreq)*Mb);
		Fd=Math.round(log_pos((ib+1)/(bands-1),basefreq,maxfreq)*Mb);
		La=log_pos_inv(Fa/Mb,basefreq,maxfreq);
		Ld=log_pos_inv(Fd/Mb,basefreq,maxfreq);
		if(Fd>Mb/2)
			Fd=Mb/2;
		if(Fa<1)
			Fa=1;
		Mc=(Fd-Fa)*2+1;
		if(Md>Mc)
			Mc=Md;
		if(Md<Mc)
			Mc=nextsprime(Mc);
		console.log(`${ib+1}/${bands} (FFT size:${Mc}) ${Fa*sampleRate/Mb} Hz - ${Fd*sampleRate/Mb} Hz\r`);
		out[bands-ib-1]=new Float64Array(Mc).fill(0);
		for(i=0;i<Fd-Fa;i++){
			Li=log_pos_inv((i+Fa)/Mb,basefreq,maxfreq);
			Li=(Li-La)/(Ld-La);
			coef=.5-.5*Math.cos(2*Math.PI*Li);
			out[bands-ib-1][i+1]=s[i+1+Fa]*coef;
			out[bands-ib-1][Mc-1-i]=s[Mb-Fa-1-i]*coef;
		}
		h=new Float64Array(Mc).fill(0);
		for(i=0;i<Fd-Fa;i++){
			h[i+1]=out[bands-ib-1][Mc-1-i];
			h[Mc-1-i]=-out[bands-ib-1][i+1];
		}
		fft(out[bands-ib-1],out[bands-ib-1],Mc,1);
		fft(h,h,Mc,1);
		for(i=0;i<Mc;i++)
			out[bands-ib-1][i]=Math.sqrt(out[bands-ib-1][i]*out[bands-ib-1][i]+h[i]*h[i]);
		if(Mc<Md)
			out[bands-ib-1]=new Float64Array(Md).set(out[bands-ib-1]);
		if(Mc>Md){
			t=out[bands-ib-1];
			out[bands-ib-1]=blackman_downsampling(out[bands-ib-1],Mc,Md);
		}
		out[bands-ib-1]=new Float64Array(Xsize).set(out[bands-ib-1]);
	}
	console.log("\n");
	normi(out,Xsize,bands,1);
	return out;
}
function wsinc_max(length,bw){
	let i,
	    bwl, //integer transition bandwidth
	    tbw, //double transition bandwidth
	    h, //kernel
	    x, //position in the antiderivative of the Blackman function
	    coef; //coefficient obtained from the function
	tbw=bw*(length-1);
	bwl=Math.ceil(tbw); //roundup
	h=new Float64Array(length).fill(1);
	for(i=0;i<bwl;i++){
		x=i/tbw; //position calculation between 0.0 and 1.0
		coef=.42*x-(.5/(2*Math.PI))*Math.sin(2*Math.PI*x)+(.08/(4*Math.PI))*Math.sin(4*Math.PI*x); //antiderivative of the Blackman function
		coef*=1/.42;
		h[i+1]=coef;
		h[length-1-i]=coef;
	}
	return h;
}
function synt_sine(d,Xsize,bands,samplecount,samplerate,basefreq,pixpersec,bpo){
	let s,freq,filter,sband,sine=new Array(4),rphase,
	    i,ib,
	    Fc,Bc,Mh,Mn,sbsize;
	freq=freqarray(basefreq,bands,bpo);
	clocka=new Date().getTime();
	sbsize=nextsprime(Xsize*2); //In Circular mode, keep it to sbsize=Xsize*2;
	samplecount=Math.round(Xsize/pixpersec);
	console.log(`Sound duration : ${(samplecount/samplerate).toFixed(3)} s`);
	samplecount=Math.round(.5*sbsize/pixpersec); //Do not change this value as it would stretch envelopes
	s=new Float64Array(samplecount).fill(0); //allocation of the sound signal
	sband=new Float64Array(sbsize); //allocation of the shifted band
	Bc=Math.round(.25*sbsize);
	Mh=(sbsize+1)>>1;
	Mn=(samplecount+1)>>1;
	filter=wsinc_max(Mh,1/TRANSITION_BW_SYNT); //generation of the frequency-domain filter
	for(ib=0;ib<bands;ib++){
		sband.fill(0); //reset sband
		//********Frequency shifting********
		rphase=dblrand()*Math.PI; //random phase between-pi and+pi
		for(i=0;i<4;i++)// generating the random sine LUT
			sine[i]=Math.cos(i*2*Math.PI*.25+rphase);
		for(i=0;i<Xsize;i++){ //envelope sampling rate*2 and frequency shifting by 0.25
			if((i&1)===0){
				sband[i<<1]=d[bands-ib-1][i]*sine[0];
				sband[(i<<1)+1]=d[bands-ib-1][i]*sine[1];
			}else{
				sband[i<<1]=d[bands-ib-1][i]*sine[2];
				sband[(i<<1)+1]=d[bands-ib-1][i]*sine[3];
			}
		}
		//--------Frequency shifting--------
		fft(sband,sband,sbsize,0); //FFT of the envelope
		Fc=Math.round(freq[ib]*samplecount); //band's centre index(envelope's DC element)
		console.log(`${ib+1}/${bands}   ${(Fc*samplerate/samplecount).toFixed(2)} Hz\r`);
		//********Write FFT********
		for(i=1;i<Mh;i++){
			if(Fc-Bc+i>0&&Fc-Bc+i<Mn){ //if we're between frequencies 0 and 0.5 of the new signal and not at Fc
				s[i+Fc-Bc]+=sband[i]*filter[i]; //Real part
				s[samplecount-(i+Fc-Bc)]+=sband[sbsize-i]*filter[i]; //Imaginary part
			}
		}
		//--------Write FFT--------
	}
	console.log("\n");
	fft(s,s,samplecount,1); //IFFT of the final sound
	samplecount=Math.round(Xsize/pixpersec); //chopping tails by ignoring them
	normi(s,samplecount,1,1);
	return s;
}
function synt_noise(d,Xsize,bands,samplecount,samplerate,basefreq,pixpersec,bpo){
	let i,ib,il,
	    s,coef,noise,loop_size_sec=LOOP_SIZE_SEC,loop_size,loop_size_min,pink_noise,mag,phase,envelope,lut,
	    freq,maxfreq,Fa,Fd,La,Ld,Li;
	freq=freqarray(basefreq,bands,bpo);
	if(LOGBASE===1)
		maxfreq=bpo; //in linear mode, we use bpo to store the maxfreq
	else
		maxfreq=basefreq*Math.pow(LOGBASE,(bands-1)/bpo);
	clocka=new Date().getTime();
	samplecount=Math.round(Xsize/pixpersec); //calculation of the length of the final signal
	console.log(`Sound duration : ${(samplecount/samplerate).toFixed(3)} s`);
	s=new Float64Array(samplecount).fill(0); //allocation of the final signal
	envelope=new Float64Array(samplecount).fill(0); //allocation of the interpolated envelope
	//********Loop size calculation********
	loop_size=loop_size_sec*samplerate;
	if(LOGBASE===1)
		loop_size_min=Math.round(4*5/(freq[1]-freq[0])); //linear mode
	else
		loop_size_min=Math.round(2*5/((freq[0]*Math.pow(2,-1/bpo))*(1-Math.pow(2,-1/bpo)))); //estimate of the longest FIR's length
	if(loop_size_min>loop_size)
		loop_size=loop_size_min;
	loop_size=nextsprime(loop_size); //enlarge the loop_size to the next multiple of short primes
	//--------Loop size calculation--------
	//********Pink noise generation********
	pink_noise=new Float64Array(loop_size).fill(0);
	for(i=1;i<(loop_size+1)>>1;i++){
		mag=Math.pow(i,.5-.5*LOGBASE); //FIXME: something's not necessarily right with that formula
		phase=dblrand()*Math.PI; //random phase between-pi and+pi
		pink_noise[i]=mag*Math.cos(phase); //real part
		pink_noise[loop_size-i]=mag*Math.sin(phase); //imaginary part
	}
	//--------Pink noise generation--------
	noise=new Float64Array(loop_size); //allocate noise
	lut=bmsq_lut(BMSQ_LUT_SIZE); //Blackman Square look-up table initialization
	for(ib=0;ib<bands;ib++){
		console.log(`${ib+1}/${bands}\r`);
		noise.fill(0); //reset filtered noise
		//********Filtering********
		Fa=Math.round(log_pos((ib-1)/(bands-1),basefreq,maxfreq)*loop_size);
		Fd=Math.round(log_pos((ib+1)/(bands-1),basefreq,maxfreq)*loop_size);
		La=log_pos_inv(Fa/loop_size,basefreq,maxfreq);
		Ld=log_pos_inv(Fd/loop_size,basefreq,maxfreq);
		if(Fd>loop_size/2)
			Fd=loop_size/2; //stop reading if reaching the Nyquist frequency
		if(Fa<1)
			Fa=1;
		console.log(`${ib+1}/${bands}   ${(Fa*samplerate/loop_size).toFixed(2)} Hz - ${(Fd*samplerate/loop_size).toFixed(2)} Hz\r`);
		for(i=Fa;i<Fd;i++){
			Li=log_pos_inv(i/loop_size,basefreq,maxfreq); //calculation of the logarithmic position
			Li=(Li-La)/(Ld-La);
			coef=.5-.5*Math.cos(2*Math.PI*Li); //Hann function
			noise[i+1]=pink_noise[i+1]*coef;
			noise[loop_size-1-i]=pink_noise[loop_size-1-i]*coef;
		}
		//--------Filtering--------
		fft(noise,noise,loop_size,1); //IFFT of the filtered noise
		envelope.fill(0); //blank the envelope
		blackman_square_interpolation(d[bands-ib-1],envelope,Xsize,samplecount,lut,BMSQ_LUT_SIZE); //interpolation of the envelope
		il=0;
		for(i=0;i<samplecount;i++){
			s[i]+=envelope[i]*noise[il]; //modulation
			il++; //increment loop iterator
			if(il===loop_size) //if the array iterator has reached the end of the array, reset it
				il=0;
		}
	}
	console.log("\n");
	normi(s,samplecount,1,1);
	return s;
}
function brightness_control(image,width,height,ratio){
	//Adjusts the brightness of the image using a logarithmic-like transformation
	for(let iy=0;iy<width;iy++){
		for(let ix=0;ix<height;ix++){
			image[iy][ix]=Math.pow(image[iy][ix],ratio);
		}
	}
}

//■■■■ Virtual ANS (Alexander Nikolayevich Scriabin) by Alexander Zolotov → pixilang\lib_dsp\dsp.h & dsp_functions.cpp ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
const M_PI=Math.PI;
//Half of sine table (assuming it's predefined)
const g_hsin_tab=new Uint8Array([
	0,3,6,9,12,15,18,21,25,28,31,34,37,40,43,46,
	49,52,56,59,62,65,68,71,74,77,80,83,86,89,92,95,
	97,100,103,106,109,112,115,117,120,123,126,128,131,134,136,139,
	142,144,147,149,152,154,157,159,162,164,167,169,171,174,176,178,
	180,183,185,187,189,191,193,195,197,199,201,203,205,207,209,211,
	212,214,216,217,219,221,222,224,225,227,228,229,231,232,233,235,
	236,237,238,239,240,242,243,243,244,245,246,247,248,249,249,250,
	251,251,252,252,253,253,254,254,254,255,255,255,255,255,255,255,
	255,255,255,255,255,255,255,255,254,254,254,253,253,252,252,251,
	251,250,249,249,248,247,246,245,245,244,243,242,241,240,238,237,
	236,235,234,232,231,230,228,227,225,224,222,221,219,218,216,214,
	213,211,209,207,205,203,201,200,198,196,194,191,189,187,185,183,
	181,179,176,174,172,169,167,165,162,160,157,155,152,150,147,145,
	142,139,137,134,131,129,126,123,120,118,115,112,109,106,104,101,
	98,95,92,89,86,83,80,77,74,71,68,65,62,59,56,53,
	50,47,44,41,37,34,31,28,25,22,19,16,12,9,6,3
]);
const FFT_FLAG_INVERSE=1;
//FFT function for both float and double arrays
function fft(flags,fi,fr,size){do_fft(flags,fi,fr,size);}
//Template function for FFT
function do_fft(flags,fi,fr,size){
	let size2=size/2,
	    isign=-1;
	if(flags&FFT_FLAG_INVERSE)
		isign=1;
	//Bit-reversal permutation
	for(let i=1,j=size2;i<size-1;i++){
		if(i<j){
			[fr[j],fr[i]]=[fr[i],fr[j]];
			[fi[j],fi[i]]=[fi[i],fi[j]];
		}
		let k=size2;
		while(k<=j){
			j-=k;
			k>>=1;
		}
		j+=k;
	}
	//Danielson-Lanczos section
	let mmax=1;
	while(mmax<size){
		let istep=mmax<<1,
		    theta=isign*M_PI/mmax,
		    wtemp=Math.sin(.5*theta),
		    wpr=-2*wtemp*wtemp,
		    wpi=Math.sin(theta),
		    wr=1,
		    wi=0;
		for(let m=0;m<mmax;m++){
			for(let i=m;i<size;i+=istep){
				let j=i+mmax,
				    tr=wr*fr[j]-wi*fi[j],
				    ti=wr*fi[j]+wi*fr[j];
				fr[j]=fr[i]-tr;
				fi[j]=fi[i]-ti;
				fr[i]+=tr;
				fi[i]+=ti;
			}
			wtemp=wr;
			wr+=wtemp*wpr-wi*wpi;
			wi+=wtemp*wpi+wi*wpr;
		}
		mmax=istep;
	}
	if(flags&FFT_FLAG_INVERSE){
		for(let i=0;i<size;i++){
			fr[i]/=size;
			fi[i]=-fi[i]/size;
		}
	}
}
//DSP curve function
function dsp_curve(x,type){
	switch(type){
		case'exponential1':let y2=32768-x;return 32768-((y2*y2)>>15);break;
		case'exponential2':return(x*x)>>15;break;
		case'spline':return Math.round(((Math.sin((x/32768)*M_PI-M_PI/2)+1)/2)*32768);break;
		case'rectangular':return(x<16384)?0:32768;break;
		default:break;
	}
	return x;
}
//DSP window function
function dsp_window_function(dest,size,type){
	let size2=size-1;
	switch(type){
		case'triangular':
			for(let i=0;i<size;i++)
				dest[i]=(i<size/2)?(i/size*2):((1.0-i/size)*2);
			break;
		case'sine':
			for(let i=0;i<size;i++)
				dest[i]=Math.sin(M_PI*i/size2);
			break;
		case'hann':
			for(let i=0;i<size;i++){
				let v=Math.sin(M_PI*i/size2);
				dest[i]=v*v;
			}
			break;
		case'blackman':
			for(let i=0;i<size;i++)
				dest[i]=0.42-0.5*Math.cos(2*M_PI*i/size2)+0.08*Math.cos(4*M_PI*i/size2);
			break;
		case'nuttall':
			for(let i=0;i<size;i++){
				let v=(2.0*M_PI*i)/size2;
				dest[i]=0.355768-0.487396*Math.cos(v)+0.144232*Math.cos(2*v)-0.012604*Math.cos(3*v);
			}
			break;
		default:break;
	}
	let sum=dest.reduce((a,b)=>a+b,0);
	return size/sum;//Amplitude correction
}
//Catmull-Rom cubic spline interpolation
function catmull_rom_spline_interp_int16(y0,y1,y2,y3,x){
	x>>=2;
	let a=(3*(y1-y2)-y0+y3)/2,
	    b=2*y2+y0-(5*y1+y3)/2,
	    c2=(y2-y0)/2;
	return(((((((a*x)>>13)+b)*x)>>13)+c2)*x)>>13+y1;
}
function catmull_rom_spline_interp_float32(y0,y1,y2,y3,x){
	let a=(3*(y1-y2)-y0+y3)/2,
	    b=2*y2+y0-(5*y1+y3)/2,
	    c2=(y2-y0)/2;
	return((a*x+b)*x+c2)*x+y1;
}
//Linear interpolation
function DSP_LINEAR_INTERP(Y0,Y1,X,X_ONE){return((Y0*(X_ONE-X))+(Y1*X))/X_ONE;}

/*
There's also pixilang\lib_sunvox\psynth\psynths_fft.cpp…
also check DtBlkFx (lonni.neocities.org/DtBlkFx_Revision, github.com/dozius/DtBlkFx)
*/
Leave a Comment