Incorrect output for realTransform (size >4)

[xenova]

Hi there 👋 While experimenting with your library, I encountered an issue where the output would be different (to numpy's implementation) by a single number. I would greatly appreciate any advice on the matter, thanks!

Reproduction:

import FFT from 'https://cdn.jsdelivr.net/npm/[email protected]/+esm'

const input = new Float32Array(8);
input.fill(1);
const f = new FFT(input.length);

const out = f.createComplexArray();
const a = performance.now();
f.realTransform(out, input);
const b = performance.now();
console.log(out)
// [8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0]
//                                      ^
//                                  incorrect

JSFiddle link (including a comparison with another library, which matches the numpy implementation).
https://jsfiddle.net/9ex3bkzt/

[xenova]

If it helps debug, it does work for input size = 4:

import FFT from 'https://cdn.jsdelivr.net/npm/[email protected]/+esm'

const input = new Float32Array(4);
input.fill(1);
const f = new FFT(input.length);

const out = f.createComplexArray();
const a = performance.now();
f.realTransform(out, input);
const b = performance.now();
console.log(b-a, out)
// [4, 0, 0, 0, 0, 0, 0, 0]
// correct

[xenova]

Further investigation shows this only occurs for realTransform:

import FFT from 'https://cdn.jsdelivr.net/npm/[email protected]/+esm'

const input = new Float32Array([1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0]);
const f = new FFT(input.length/2);

const out = f.createComplexArray();
const a = performance.now();
f.transform(out, input);
const b = performance.now();
console.log(b-a, Array.from(out))
// [8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

[xenova]

Also, if you replace this block of code:

fft.js/lib/fft.js

Lines 339 to 439 in 4a18cf8

	// Loop through offsets in the data
	for (outOff = 0; outOff < size; outOff += len) {
	for (var i = 0, k = 0; i <= hquarterLen; i += 2, k += step) {
	var A = outOff + i;
	var B = A + quarterLen;
	var C = B + quarterLen;
	var D = C + quarterLen;
	// Original values
	var Ar = out[A];
	var Ai = out[A + 1];
	var Br = out[B];
	var Bi = out[B + 1];
	var Cr = out[C];
	var Ci = out[C + 1];
	var Dr = out[D];
	var Di = out[D + 1];
	// Middle values
	var MAr = Ar;
	var MAi = Ai;
	var tableBr = table[k];
	var tableBi = inv * table[k + 1];
	var MBr = Br * tableBr - Bi * tableBi;
	var MBi = Br * tableBi + Bi * tableBr;
	var tableCr = table[2 * k];
	var tableCi = inv * table[2 * k + 1];
	var MCr = Cr * tableCr - Ci * tableCi;
	var MCi = Cr * tableCi + Ci * tableCr;
	var tableDr = table[3 * k];
	var tableDi = inv * table[3 * k + 1];
	var MDr = Dr * tableDr - Di * tableDi;
	var MDi = Dr * tableDi + Di * tableDr;
	// Pre-Final values
	var T0r = MAr + MCr;
	var T0i = MAi + MCi;
	var T1r = MAr - MCr;
	var T1i = MAi - MCi;
	var T2r = MBr + MDr;
	var T2i = MBi + MDi;
	var T3r = inv * (MBr - MDr);
	var T3i = inv * (MBi - MDi);
	// Final values
	var FAr = T0r + T2r;
	var FAi = T0i + T2i;
	var FBr = T1r + T3i;
	var FBi = T1i - T3r;
	out[A] = FAr;
	out[A + 1] = FAi;
	out[B] = FBr;
	out[B + 1] = FBi;
	// Output final middle point
	if (i === 0) {
	var FCr = T0r - T2r;
	var FCi = T0i - T2i;
	out[C] = FCr;
	out[C + 1] = FCi;
	continue;
	}
	// Do not overwrite ourselves
	if (i === hquarterLen)
	continue;
	// In the flipped case:
	// MAi = -MAi
	// MBr=-MBi, MBi=-MBr
	// MCr=-MCr
	// MDr=MDi, MDi=MDr
	var ST0r = T1r;
	var ST0i = -T1i;
	var ST1r = T0r;
	var ST1i = -T0i;
	var ST2r = -inv * T3i;
	var ST2i = -inv * T3r;
	var ST3r = -inv * T2i;
	var ST3i = -inv * T2r;
	var SFAr = ST0r + ST2r;
	var SFAi = ST0i + ST2i;
	var SFBr = ST1r + ST3i;
	var SFBi = ST1i - ST3r;
	var SA = outOff + quarterLen - i;
	var SB = outOff + halfLen - i;
	out[SA] = SFAr;
	out[SA + 1] = SFAi;
	out[SB] = SFBr;
	out[SB + 1] = SFBi;
	}
	}

with this block of code:

fft.js/lib/fft.js

Lines 151 to 222 in 4a18cf8

	// Loop through offsets in the data
	for (outOff = 0; outOff < size; outOff += len) {
	// Full case
	var limit = outOff + quarterLen;
	for (var i = outOff, k = 0; i < limit; i += 2, k += step) {
	const A = i;
	const B = A + quarterLen;
	const C = B + quarterLen;
	const D = C + quarterLen;
	// Original values
	const Ar = out[A];
	const Ai = out[A + 1];
	const Br = out[B];
	const Bi = out[B + 1];
	const Cr = out[C];
	const Ci = out[C + 1];
	const Dr = out[D];
	const Di = out[D + 1];
	// Middle values
	const MAr = Ar;
	const MAi = Ai;
	const tableBr = table[k];
	const tableBi = inv * table[k + 1];
	const MBr = Br * tableBr - Bi * tableBi;
	const MBi = Br * tableBi + Bi * tableBr;
	const tableCr = table[2 * k];
	const tableCi = inv * table[2 * k + 1];
	const MCr = Cr * tableCr - Ci * tableCi;
	const MCi = Cr * tableCi + Ci * tableCr;
	const tableDr = table[3 * k];
	const tableDi = inv * table[3 * k + 1];
	const MDr = Dr * tableDr - Di * tableDi;
	const MDi = Dr * tableDi + Di * tableDr;
	// Pre-Final values
	const T0r = MAr + MCr;
	const T0i = MAi + MCi;
	const T1r = MAr - MCr;
	const T1i = MAi - MCi;
	const T2r = MBr + MDr;
	const T2i = MBi + MDi;
	const T3r = inv * (MBr - MDr);
	const T3i = inv * (MBi - MDi);
	// Final values
	const FAr = T0r + T2r;
	const FAi = T0i + T2i;
	const FCr = T0r - T2r;
	const FCi = T0i - T2i;
	const FBr = T1r + T3i;
	const FBi = T1i - T3r;
	const FDr = T1r - T3i;
	const FDi = T1i + T3r;
	out[A] = FAr;
	out[A + 1] = FAi;
	out[B] = FBr;
	out[B + 1] = FBi;
	out[C] = FCr;
	out[C + 1] = FCi;
	out[D] = FDr;
	out[D + 1] = FDi;
	}
	}

The output is as expected, so it seems to be an error in the optimization.

[daniele-roncaglioni]

I think you have to call f.completeSpectrum(out); as mentioned in the readme, then it correctly sets the values in the right part of the array.

[xenova]

Oh wow, that was it... thanks so much @daniele-roncaglioni!