Various collection of (transpiled) FFT algorithms
/* (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