From cf918942f3ed58af40edb54c0ae2115356512c2f Mon Sep 17 00:00:00 2001 From: sbosse Date: Mon, 21 Jul 2025 23:09:18 +0200 Subject: [PATCH] Mon 21 Jul 22:43:21 CEST 2025 --- js/numerics/fft.js | 447 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 447 insertions(+) create mode 100644 js/numerics/fft.js diff --git a/js/numerics/fft.js b/js/numerics/fft.js new file mode 100644 index 0000000..9224020 --- /dev/null +++ b/js/numerics/fft.js @@ -0,0 +1,447 @@ +/*===========================================================================*\ + * Fast Fourier Transform (Cooley-Tukey Method) + * + * (c) Vail Systems. Joshua Jung and Ben Bryan. 2015 + * + * This code is not designed to be highly optimized but as an educational + * tool to understand the Fast Fourier Transform. + * + * Can be used with Vectors (typed arrays), Arrays + * + * +\*===========================================================================*/ +var Vector = Require('numerics/vector'); +var Matrix = Require('numerics/matrix'); + +//------------------------------------------------ +// Note: Some of this code is not optimized and is +// primarily designed as an educational and testing +// tool. +// To get high performace would require transforming +// the recursive calls into a loop and then loop +// unrolling. All of this is best accomplished +// in C or assembly. +//------------------------------------------------- + +//------------------------------------------------- +// The following code assumes a complex number is +// an array: [real, imaginary] +//------------------------------------------------- +var complex = { + //------------------------------------------------- + // Add two complex numbers + //------------------------------------------------- + add : function (a, b) + { + return [a[0] + b[0], a[1] + b[1]]; + }, + + //------------------------------------------------- + // Subtract two complex numbers + //------------------------------------------------- + subtract : function (a, b) + { + return [a[0] - b[0], a[1] - b[1]]; + }, + + //------------------------------------------------- + // Multiply two complex numbers + // + // (a + bi) * (c + di) = (ac - bd) + (ad + bc)i + //------------------------------------------------- + multiply : function (a, b) + { + return [(a[0] * b[0] - a[1] * b[1]), + (a[0] * b[1] + a[1] * b[0])]; + }, + + //------------------------------------------------- + // Calculate |a + bi| + // + // sqrt(a*a + b*b) + //------------------------------------------------- + magnitude : function (c) + { + return Math.sqrt(c[0]*c[0] + c[1]*c[1]); + }, + + phase : function (c) + { + return c[0]!=0?Math.atan(c[1]/c[0])*180/Math.PI:(c[1]>0?90:-90); + } + +} +var mapExponent = {}; +var fftUtil = { + //------------------------------------------------- + // By Eulers Formula: + // + // e^(i*x) = cos(x) + i*sin(x) + // + // and in DFT: + // + // x = -2*PI*(k/N) + //------------------------------------------------- + exponent : function (k, N) { + var x = -2 * Math.PI * (k / N); + + mapExponent[N] = mapExponent[N] || {}; + mapExponent[N][k] = mapExponent[N][k] || [Math.cos(x), Math.sin(x)];// [Real, Imaginary] + + return mapExponent[N][k]; + }, + + //------------------------------------------------- + // Calculate FFT Magnitude for complex numbers. + //------------------------------------------------- + fftMag : function (fftBins) { + if (isArray(fftBins)) { + var ret = fftBins.map(complex.magnitude); + return ret.slice(0, ret.length / 2); + } else if (isVector(fftBins) || isMatrix(fftBins)) { + // complex matrix (2 columns) + if (fftBins.columns != 2) throw "fft.fftMag: Complex matrix columns != 2"; + var ret = Vector(fftBins.rows,{dtn:fftBins.dtn}) + return ret.eval(function (v,i) { + return complex.magnitude([fftBins.get(i,0),fftBins.get(i,1)]) }); + } + }, + fftPha : function (fftBins) { + if (isArray(fftBins)) { + var ret = fftBins.map(complex.phase); + return ret.slice(0, ret.length / 2); + } else if (isVector(fftBins) || isMatrix(fftBins)) { + // complex matrix (2 columns) + if (fftBins.columns != 2) throw "fft.fftMag: Complex matrix columns != 2"; + var ret = Vector(fftBins.rows,{dtn:fftBins.dtn}) + return ret.eval(function (v,i) { + return complex.phase([fftBins.get(i,0),fftBins.get(i,1)]) }); + } + }, + + //------------------------------------------------- + // Calculate Frequency Bins + // + // Returns an array of the frequencies (in hertz) of + // each FFT bin provided, assuming the sampleRate is + // samples taken per second. + //------------------------------------------------- + fftFreq : function (fftBins, sampleRate) { + var stepFreq = sampleRate / (fftBins.length); + var ret = fftBins.slice(0, fftBins.length / 2); + + return ret.map(function (__, ix) { + return ix * stepFreq; + }); + } +} + +// Bit-twiddle +var REVERSE_TABLE = new Array(256); +var INT_BITS = 32; //Number of bits in an integer + +(function(tab) { + for(var i=0; i<256; ++i) { + var v = i, r = i, s = 7; + for (v >>>= 1; v; v >>>= 1) { + r <<= 1; + r |= v & 1; + --s; + } + tab[i] = (r << s) & 0xff; + } +})(REVERSE_TABLE); + +var twiddle = { + + //Constants + INT_BITS : INT_BITS, + INT_MAX : 0x7fffffff, + INT_MIN : -1<<(INT_BITS-1), + + //Returns -1, 0, +1 depending on sign of x + sign : function(v) { + return (v > 0) - (v < 0); + }, + + //Computes absolute value of integer + abs : function(v) { + var mask = v >> (INT_BITS-1); + return (v ^ mask) - mask; + }, + + //Computes minimum of integers x and y + min : function(x, y) { + return y ^ ((x ^ y) & -(x < y)); + }, + + //Computes maximum of integers x and y + max : function(x, y) { + return x ^ ((x ^ y) & -(x < y)); + }, + + //Checks if a number is a power of two + isPow2 : function(v) { + return !(v & (v-1)) && (!!v); + }, + + //Computes log base 2 of v + log2 : function(v) { + var r, shift; + r = (v > 0xFFFF) << 4; v >>>= r; + shift = (v > 0xFF ) << 3; v >>>= shift; r |= shift; + shift = (v > 0xF ) << 2; v >>>= shift; r |= shift; + shift = (v > 0x3 ) << 1; v >>>= shift; r |= shift; + return r | (v >> 1); + }, + + //Computes log base 10 of v + log10 : function(v) { + return (v >= 1000000000) ? 9 : (v >= 100000000) ? 8 : (v >= 10000000) ? 7 : + (v >= 1000000) ? 6 : (v >= 100000) ? 5 : (v >= 10000) ? 4 : + (v >= 1000) ? 3 : (v >= 100) ? 2 : (v >= 10) ? 1 : 0; + }, + + //Counts number of bits + popCount : function(v) { + v = v - ((v >>> 1) & 0x55555555); + v = (v & 0x33333333) + ((v >>> 2) & 0x33333333); + return ((v + (v >>> 4) & 0xF0F0F0F) * 0x1010101) >>> 24; + }, + + //Counts number of trailing zeros + countTrailingZeros : function (v) { + var c = 32; + v &= -v; + if (v) c--; + if (v & 0x0000FFFF) c -= 16; + if (v & 0x00FF00FF) c -= 8; + if (v & 0x0F0F0F0F) c -= 4; + if (v & 0x33333333) c -= 2; + if (v & 0x55555555) c -= 1; + return c; + }, + + //Rounds to next power of 2 + nextPow2 : function(v) { + v += v === 0; + --v; + v |= v >>> 1; + v |= v >>> 2; + v |= v >>> 4; + v |= v >>> 8; + v |= v >>> 16; + return v + 1; + }, + + //Rounds down to previous power of 2 + prevPow2 : function(v) { + v |= v >>> 1; + v |= v >>> 2; + v |= v >>> 4; + v |= v >>> 8; + v |= v >>> 16; + return v - (v>>>1); + }, + + //Computes parity of word + parity : function(v) { + v ^= v >>> 16; + v ^= v >>> 8; + v ^= v >>> 4; + v &= 0xf; + return (0x6996 >>> v) & 1; + }, + + + //Reverse bits in a 32 bit word + reverse : function(v) { + return (REVERSE_TABLE[ v & 0xff] << 24) | + (REVERSE_TABLE[(v >>> 8) & 0xff] << 16) | + (REVERSE_TABLE[(v >>> 16) & 0xff] << 8) | + REVERSE_TABLE[(v >>> 24) & 0xff]; + }, + + //Interleave bits of 2 coordinates with 16 bits. Useful for fast quadtree codes + interleave2 : function(x, y) { + x &= 0xFFFF; + x = (x | (x << 8)) & 0x00FF00FF; + x = (x | (x << 4)) & 0x0F0F0F0F; + x = (x | (x << 2)) & 0x33333333; + x = (x | (x << 1)) & 0x55555555; + + y &= 0xFFFF; + y = (y | (y << 8)) & 0x00FF00FF; + y = (y | (y << 4)) & 0x0F0F0F0F; + y = (y | (y << 2)) & 0x33333333; + y = (y | (y << 1)) & 0x55555555; + + return x | (y << 1); + }, + + //Extracts the nth interleaved component + deinterleave2 : function(v, n) { + v = (v >>> n) & 0x55555555; + v = (v | (v >>> 1)) & 0x33333333; + v = (v | (v >>> 2)) & 0x0F0F0F0F; + v = (v | (v >>> 4)) & 0x00FF00FF; + v = (v | (v >>> 16)) & 0x000FFFF; + return (v << 16) >> 16; + }, + + + //Interleave bits of 3 coordinates, each with 10 bits. Useful for fast octree codes + interleave3 : function(x, y, z) { + x &= 0x3FF; + x = (x | (x<<16)) & 4278190335; + x = (x | (x<<8)) & 251719695; + x = (x | (x<<4)) & 3272356035; + x = (x | (x<<2)) & 1227133513; + + y &= 0x3FF; + y = (y | (y<<16)) & 4278190335; + y = (y | (y<<8)) & 251719695; + y = (y | (y<<4)) & 3272356035; + y = (y | (y<<2)) & 1227133513; + x |= (y << 1); + + z &= 0x3FF; + z = (z | (z<<16)) & 4278190335; + z = (z | (z<<8)) & 251719695; + z = (z | (z<<4)) & 3272356035; + z = (z | (z<<2)) & 1227133513; + + return x | (z << 2); + }, + + //Extracts nth interleaved component of a 3-tuple + deinterleave3 : function(v, n) { + v = (v >>> n) & 1227133513; + v = (v | (v>>>2)) & 3272356035; + v = (v | (v>>>4)) & 251719695; + v = (v | (v>>>8)) & 4278190335; + v = (v | (v>>>16)) & 0x3FF; + return (v<<22)>>22; + }, + + //Computes next combination in colexicographic order (this is mistakenly called nextPermutation on the bit twiddling hacks page) + nextCombination : function(v) { + var t = v | (v - 1); + return (t + 1) | (((~t & -~t) - 1) >>> (twiddle.countTrailingZeros(v) + 1)); + } + +} + +function checkpow2(info,vector) { + if (Math.floor(Math.log2(vector.length)) != Math.log2(vector.length)) + throw ('fft.'+info+' error: vector length must be a power of 2') +} + +module.exports = { + //------------------------------------------------- + // Calculate FFT for vector where vector.length + // is assumed to be a power of 2. + //------------------------------------------------- + fft: function fft(vector) { + if (vector.data) vector=vector.data; // Matrix|Vector + + checkpow2('fft',vector); + + var X = [], + N = vector.length; + + // Base case is X = x + 0i since our input is assumed to be real only. + if (N == 1) { + if (Array.isArray(vector[0])) //If input vector contains complex numbers + return [[vector[0][0], vector[0][1]]]; + else + return [[vector[0], 0]]; + } + + // Recurse: all even samples + var X_evens = fft(vector.filter(even)), + + // Recurse: all odd samples + X_odds = fft(vector.filter(odd)); + + // Now, perform N/2 operations! + for (var k = 0; k < N / 2; k++) { + // t is a complex number! + var t = X_evens[k], + e = complex.multiply(fftUtil.exponent(k, N), X_odds[k]); + + X[k] = complex.add(t, e); + X[k + (N / 2)] = complex.subtract(t, e); + } + + function even(__, ix) { + return ix % 2 == 0; + } + + function odd(__, ix) { + return ix % 2 == 1; + } + + return X; + }, + //------------------------------------------------- + // Calculate FFT for vector where vector.length + // is assumed to be a power of 2. This is the in- + // place implementation, to avoid the memory + // footprint used by recursion. + //------------------------------------------------- + fftInPlace: function(vector) { + if (vector.data) vector=vector.data; // Matrix|Vector + checkpow2('fftInPlace',vector); + + var N = vector.length; + + var trailingZeros = twiddle.countTrailingZeros(N); //Once reversed, this will be leading zeros + + // Reverse bits + for (var k = 0; k < N; k++) { + var p = twiddle.reverse(k) >>> (twiddle.INT_BITS - trailingZeros); + if (p > k) { + var complexTemp = [vector[k], 0]; + vector[k] = vector[p]; + vector[p] = complexTemp; + } else { + vector[p] = [vector[p], 0]; + } + } + + //Do the DIT now in-place + for (var len = 2; len <= N; len += len) { + for (var i = 0; i < len / 2; i++) { + var w = fftUtil.exponent(i, len); + for (var j = 0; j < N / len; j++) { + var t = complex.multiply(w, vector[j * len + i + len / 2]); + vector[j * len + i + len / 2] = complex.subtract(vector[j * len + i], t); + vector[j * len + i] = complex.add(vector[j * len + i], t); + } + } + } + }, + ifft: function ifft(signal){ + if (signal.data) signal=signal.data; // Matrix + checkpow2('ifft',signal); + //Interchange real and imaginary parts + var csignal=[]; + for(var i=0; i