//------- FFT functions

// njr: found fft in javascript, saving me translating:
/*=========================================================
  Filename: dspUtils-08.js
  Rev: 8
  By: A.R.Collins
  Description: Digital Signal Processing utilities.

  License: Released into the public domain
  latest version at
  <http://www.arc.id.au>
  =========================================================*/
function fft(Ind, Npair, Ar, Ai) {
/*=========================================
 * Calculate the floating point complex FFT
 * Ind = +1 => FORWARD FFT
 * Ind = -l => INVERSE FFT
 * Data is passed in Npair Complex pairs
 * where Npair is power of 2 (2^N)
 * data is indexed from 0 to Npair-1
 * Real data in Ar
 * Imag data in Ai.
 *
 * Output data is returned in the same arrays,
 * DC in bin 0, +ve freqs in bins 1..Npair/2
 * -ve freqs in Npair/2+1 .. Npair-1.
 *
 * ref: Rabiner & Gold
 * "THEORY AND APPLICATION OF DIGITAL
 *  SIGNAL PROCESSING" p367
 *
 * Translated to JavaScript by A.R.Collins
 * <http://www.arc.id.au>
 *========================================*/

	var Pi = Math.PI;
	var Num1, Num2, I, J, K, L, M, Le, Le1;
	var Tr, Ti, Ur, Ui, Xr, Xi;

	M = isPwrOf2(Npair);
	if (M < 0) {
		alert("Npair must be power of 2 from 4 to 4096");
	return;
	}

	Num1 = Npair-1;
	Num2 = Npair/2;
	// if IFT conjugate prior to transforming:
	if (Ind < 0) {
		for (I = 0; I < Npair; I++)
			Ai[I] *= -1;
	}

	J = 0;    // In place bit reversal of input data
	for (I = 0; I < Num1; I++) {
		if (I < J) {
			Tr = Ar[J];
			Ti = Ai[J];
			Ar[J] = Ar[I];
			Ai[J] = Ai[I];
			Ar[I] = Tr;
			Ai[I] = Ti;
		}
		K = Num2;
		while (K < J+1) {
			J = J-K;
			K = K/2;
		}
		J = J+K;
	}

	Le = 1;
	for (L = 1; L <= M; L++) {
		Le1 = Le;
		Le += Le;
		Ur = 1;
		Ui = 0;
		Wr = Math.cos(Pi/Le1);
		Wi = -Math.sin(Pi/Le1);
		for (J = 1; J <= Le1; J++) {
			for (I = J-1; I <= Num1; I += Le) {
				Ip = I+Le1;
				Tr = Ar[Ip]*Ur-Ai[Ip]*Ui;
				Ti = Ar[Ip]*Ui+Ai[Ip]*Ur;
				Ar[Ip] = Ar[I]-Tr;
				Ai[Ip] = Ai[I]-Ti;
				Ar[I] = Ar[I]+Tr;
				Ai[I] = Ai[I]+Ti;
			}
			Xr = Ur*Wr-Ui*Wi;
			Xi = Ur*Wi+Ui*Wr;
			Ur = Xr;
			Ui = Xi;
		}
	}
	// conjugate and normalise
	if (Ind < 0) {
		for (I = 0; I < Npair; I++)
			Ai[I] *= -1;
	}
	else {
		for (I = 0; I < Npair; I++) {
			Ar[I] /= Npair;
			Ai[I] /= Npair;
		}
	}
}

function isPwrOf2(n) {
	var m = -1;
	for (m = 2; m < 13; m++)
		if (Math.pow(2,m) == n)
			return m;
	return -1;
}


//
// getMagDB: perform fft and return magnitude of entire spectrum
// note: length of returned array is one greater than fft size, to include nyquist
//
// njr Nov 2010
//
function getMagDB(signal, fftSize) {
	var halfSignalSpec = doRealFFT(signal, fftSize);

	// impulse for normalization
	var impulse = [1];
	impulse[0] = 1.0;
	halfImpulseSpec = doRealFFT(impulse, fftSize);
	
	var halfSigMag = getMag(halfSignalSpec);
	var halfImpMag = getMag(halfImpulseSpec);

	// magnitude in dB
	var dbMult = 20 / Math.log(10);
	var mag = [fftSize + 1];
	for (var idx = 0, len = halfSigMag.length; idx < len; idx++) {
		var dbVal = dbMult * Math.log(halfSigMag[idx] / halfImpMag[idx]);
		if (dbVal == -Infinity)	// log(0) returns -Infinity
			dbVal = -200;
		mag[idx] = mag[ fftSize - idx ] = dbVal;
	}

	return mag;
}


//
// getMag: get magnitude of half spectrum
// returns array one greater than the size of the complex input, to include nyquist
//
// njr Nov 2010
//
function getMag(halfSpectrum){
	var rPart = halfSpectrum[0];
	var iPart = halfSpectrum[1];
	var len = rPart.length;
	var data = [len + 1];
	for (idx = 1; idx < len; idx++)
		data[idx] = Math.sqrt(rPart[idx] * rPart[idx] + iPart[idx] * iPart[idx]);
	data[0] = Math.abs(rPart[0]) * 2;
	data[len] = Math.abs(iPart[0]) * 2;
	return data;
}


//
// doRealFFT: do real fft returns half spectrum of fft of real data (nyauist in iPart[0])
// input may be smaller than len, and will be padded; len must be power of 2
//
// njr Nov 2010
//
function doRealFFT(input, len) {
	// copy and zero pad to len
	var inputLen = input.length;
	var data = [len];
	for (idx = 0; idx < inputLen; idx++)
		data[idx] = input[idx];
	for (idx = inputLen; idx < len; idx++)
		data[idx] = 0;

	// perform fft on real data for half spectrum
	var halfLen = len >> 1;
	var rPart = [halfLen];
	var iPart = [halfLen];
	for (idx = 0, jdx = 0; idx < halfLen; idx++) {
		iPart[idx] = data[jdx++];
		rPart[idx] = data[jdx++];
	}
	fft(1, len >> 1, rPart, iPart);
	complexToReal(rPart, iPart);

	return [ rPart, iPart ];
}


//
// realToComplex: convert real data in complex array to real, in place
//
// njt Nov 2010 - based on Chamberlin
//
function realToComplex(rD, iD) {
	var n = rD.length;
	var piOverN = Math.PI / n;
	var t1 = rD[0];
	var t2 = iD[0];
	rD[0] = t1 + t2;
	iD[0] = t1 - t2;
	for (n1 = 1, n1max = n * 0.5, n2 = n - 1; n1 <= n1max; n1++, n2--) {
		var temp = -n1 * piOverN;
		var c = Math.cos(temp);
		var s = Math.sin(temp);
		var t1 = (rD[n1] + rD[n2]) * 0.5;
		var t4 = (iD[n1] - iD[n2]) * 0.5;
		var t5 = (rD[n2] - rD[n1]) * 0.5;
		var t6 = (-iD[n1] - iD[n2]) * 0.5;
		var t2 = t5 * s + t6 * c;
		var t3 = t6 * s - t5 * c;
		rD[n1] = t1 + t2;
		rD[n2] = t1 - t2;
		iD[n1] = t3 + t4;
		iD[n2] = t3 - t4;
	}
}


//
// complexToReal: convert complex array to real data, in place
//
// njt Nove 2010 - based on Chamberlin
//
function complexToReal(rD, iD) {
	var n = rD.length;
	var oneOverN = 1 / n;
	var piOverN = Math.PI * oneOverN;
	t1 = rD[0];
	t2 = iD[0];
	rD[0] = (t1 + t2) * oneOverN * 0.5;
	iD[0] = (t1 - t2) * oneOverN * 0.5;
	for (n1 = 1, n1max = n * 0.5, n2 = n - 1; n1 <= n1max; n1++, n2--) {
		var temp = -n1 * piOverN;
		var c = Math.cos(temp);
		var s = Math.sin(temp);
		var t1 = (rD[n1] + rD[n2]) * 0.5;
		var t2 = (rD[n1] - rD[n2]) * 0.5;
		var t3 = (iD[n1] + iD[n2]) * 0.5;
		var t4 = (iD[n1] - iD[n2]) * 0.5;
		var t5 = t2 * s - t3 * c;
		var t6 = t2 * c + t3 * s;
		rD[n1] = (t1 - t5) * oneOverN;
		rD[n2] = (t1 + t5) * oneOverN;
		iD[n1] = (t4 - t6) * oneOverN;
		iD[n2] = (-t4 - t6) * oneOverN;
	}
}

//------- Kaiser window

//
// kaiser: fills a table with the right half of a kaiser window
// win points to window array (length must be at least numTaps + 1 >> 1)
//
// njr Nov 2010
//
function kaiserHalf(beta, numTaps)
{
	var rm;
	var evenOffset = 0.5 * (!(numTaps & 1));
	var bd = bessel(beta);
	var halfTaps = (numTaps + 1) >> 1;
	var win = [halfTaps];
	for (var j = 0; j < halfTaps; ++j) {
		rm = 2 * (j + evenOffset) / (numTaps - 1);
		win[j] = bessel(beta * Math.sqrt(1 - rm * rm)) / bd;
	}
	return win;
}


//
// kaiser: returns a Kaiser window of length numTaps
//
// njr Nov 2010 - based on Rabiner & Gold
//
function kaiser(beta, numTaps)
{
	var rm;
	var evenOffset = 0.5 * (!(numTaps & 1));
	var bd = bessel(beta);
	var halfTaps = (numTaps + 1) >> 1;
	var win = [numTaps];
	var jdx = halfTaps - 1;
	var kdx = numTaps - halfTaps;
	for (var idx = 0; idx < halfTaps; ++idx) {
		rm = 2 * (idx + evenOffset) / (numTaps - 1);
		win[kdx++] = win[jdx--] = bessel(beta * Math.sqrt(1 - rm * rm)) / bd;
	}
	return win;
}


//
// bessel: evaluate the bessel function
//
// njr Nov 2010 - based on Rabiner & Gold
//
function bessel(x)
{
	var y, t, e, de, sde;
	var	i;

	y = x / 2;
	t = 1.E-08;
	e=1;
	de=1;
	for (i = 1; i < 25; ++i) {
		de = de * y / i;
		sde = de * de;
		e = e + sde;
		if ((e * t - sde) > 0.0)
			break;
	}
	return e;
}


//
// estimateBeta: given a stopband attenuation value (in +dB), estimate beta for kaiser window
//
// njr Nov 2010 - based on Rabiner & Gold
//
function estimateBeta(stopbandAtten)
{
	var beta;

	if (stopbandAtten >= 50)
		beta = 0.1102 * (stopbandAtten - 8.7);
	else if (stopbandAtten > 21)
		beta = 0.5842 * Math.pow(stopbandAtten - 21, 0.4) + 0.07886 * (stopbandAtten - 21);
	else
		beta = 0;

	return beta;
}