Mon 21 Jul 22:43:21 CEST 2025
This commit is contained in:
parent
9375c14b1f
commit
cf918942f3
447
js/numerics/fft.js
Normal file
447
js/numerics/fft.js
Normal file
|
|
@ -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<signal.length; i++){
|
||||
csignal[i]=[signal[i][1], signal[i][0]];
|
||||
}
|
||||
|
||||
//Apply fft
|
||||
var ps=module.exports.fft(csignal);
|
||||
|
||||
//Interchange real and imaginary parts and normalize
|
||||
var res=[];
|
||||
for(var j=0; j<ps.length; j++){
|
||||
res[j]=[ps[j][1]/ps.length, ps[j][0]/ps.length];
|
||||
}
|
||||
return res;
|
||||
},
|
||||
fftMag: fftUtil.fftMag,
|
||||
fftPha: fftUtil.fftPha,
|
||||
fftFreq: fftUtil.fftFreq,
|
||||
};
|
||||
Loading…
Reference in New Issue
Block a user