Use existing complex numbers implementation for Fourier.

This commit is contained in:
Oleksii Trekhleb 2018-08-13 17:57:17 +03:00
parent 6f10b0e10f
commit 7d3115edaa
3 changed files with 74 additions and 81 deletions

View File

@ -1,5 +1,6 @@
import ComplexNumber from '../complex';
import fastFourierTransform from '../fastFourierTransform'; import fastFourierTransform from '../fastFourierTransform';
import ComplexNumber from '../../complex-number/ComplexNumber';
/** /**
* @param {ComplexNumber[]} [seq1] * @param {ComplexNumber[]} [seq1]
* @param {ComplexNumber[]} [seq2] * @param {ComplexNumber[]} [seq2]
@ -20,48 +21,62 @@ function approximatelyEqual(seq1, seq2, eps) {
describe('fastFourierTransform', () => { describe('fastFourierTransform', () => {
it('should calculate the radix-2 discrete fourier transform after zero padding', () => { it('should calculate the radix-2 discrete fourier transform after zero padding', () => {
const eps = 1e-6; const eps = 1e-6;
const in1 = [new ComplexNumber(0, 0)]; const in1 = [new ComplexNumber({ real: 0, imaginary: 0 })];
const expOut1 = [new ComplexNumber(0, 0)]; const expOut1 = [new ComplexNumber({ real: 0, imaginary: 0 })];
const out1 = fastFourierTransform(in1); const out1 = fastFourierTransform(in1);
const invOut1 = fastFourierTransform(out1, true); const invOut1 = fastFourierTransform(out1, true);
expect(approximatelyEqual(expOut1, out1, eps)).toBe(true); expect(approximatelyEqual(expOut1, out1, eps)).toBe(true);
expect(approximatelyEqual(in1, invOut1, eps)).toBe(true); expect(approximatelyEqual(in1, invOut1, eps)).toBe(true);
const in2 = [new ComplexNumber(1, 2), new ComplexNumber(2, 3), const in2 = [
new ComplexNumber(8, 4)]; new ComplexNumber({ real: 1, imaginary: 2 }),
const expOut2 = [new ComplexNumber(11, 9), new ComplexNumber(-10, 0), new ComplexNumber({ real: 2, imaginary: 3 }),
new ComplexNumber(7, 3), new ComplexNumber(-4, -4)]; new ComplexNumber({ real: 8, imaginary: 4 }),
];
const expOut2 = [
new ComplexNumber({ real: 11, imaginary: 9 }),
new ComplexNumber({ real: -10, imaginary: 0 }),
new ComplexNumber({ real: 7, imaginary: 3 }),
new ComplexNumber({ real: -4, imaginary: -4 }),
];
const out2 = fastFourierTransform(in2); const out2 = fastFourierTransform(in2);
const invOut2 = fastFourierTransform(out2, true); const invOut2 = fastFourierTransform(out2, true);
expect(approximatelyEqual(expOut2, out2, eps)).toBe(true); expect(approximatelyEqual(expOut2, out2, eps)).toBe(true);
expect(approximatelyEqual(in2, invOut2, eps)).toBe(true); expect(approximatelyEqual(in2, invOut2, eps)).toBe(true);
const in3 = [new ComplexNumber(-83656.9359385182, 98724.08038374918), const in3 = [
new ComplexNumber(-47537.415125808424, 88441.58381765135), new ComplexNumber({ real: -83656.9359385182, imaginary: 98724.08038374918 }),
new ComplexNumber(-24849.657029355192, -72621.79007878687), new ComplexNumber({ real: -47537.415125808424, imaginary: 88441.58381765135 }),
new ComplexNumber(31451.27290052717, -21113.301128347346), new ComplexNumber({ real: -24849.657029355192, imaginary: -72621.79007878687 }),
new ComplexNumber(13973.90836288876, -73378.36721594246), new ComplexNumber({ real: 31451.27290052717, imaginary: -21113.301128347346 }),
new ComplexNumber(14981.520420492234, 63279.524958963884), new ComplexNumber({ real: 13973.90836288876, imaginary: -73378.36721594246 }),
new ComplexNumber(-9892.575367044381, -81748.44671677813), new ComplexNumber({ real: 14981.520420492234, imaginary: 63279.524958963884 }),
new ComplexNumber(-35933.00356823792, -46153.47157161784), new ComplexNumber({ real: -9892.575367044381, imaginary: -81748.44671677813 }),
new ComplexNumber(-22425.008561855735, -86284.24507370662), new ComplexNumber({ real: -35933.00356823792, imaginary: -46153.47157161784 }),
new ComplexNumber(-39327.43830818355, 30611.949874562706)]; new ComplexNumber({ real: -22425.008561855735, imaginary: -86284.24507370662 }),
const expOut3 = [new ComplexNumber(-203215.3322151, -100242.4827503), new ComplexNumber({ real: -39327.43830818355, imaginary: 30611.949874562706 }),
new ComplexNumber(99217.0805705, 270646.9331932), ];
new ComplexNumber(-305990.9040412, 68224.8435751),
new ComplexNumber(-14135.7758282, 199223.9878095), const expOut3 = [
new ComplexNumber(-306965.6350922, 26030.1025439), new ComplexNumber({ real: -203215.3322151, imaginary: -100242.4827503 }),
new ComplexNumber(-76477.6755206, 40781.9078990), new ComplexNumber({ real: 99217.0805705, imaginary: 270646.9331932 }),
new ComplexNumber(-48409.3099088, 54674.7959662), new ComplexNumber({ real: -305990.9040412, imaginary: 68224.8435751 }),
new ComplexNumber(-329683.0131713, 164287.7995937), new ComplexNumber({ real: -14135.7758282, imaginary: 199223.9878095 }),
new ComplexNumber(-50485.2048527, -330375.0546527), new ComplexNumber({ real: -306965.6350922, imaginary: 26030.1025439 }),
new ComplexNumber(122235.7738708, 91091.6398019), new ComplexNumber({ real: -76477.6755206, imaginary: 40781.9078990 }),
new ComplexNumber(47625.8850387, 73497.3981523), new ComplexNumber({ real: -48409.3099088, imaginary: 54674.7959662 }),
new ComplexNumber(-15619.8231136, 80804.8685410), new ComplexNumber({ real: -329683.0131713, imaginary: 164287.7995937 }),
new ComplexNumber(192234.0276101, 160833.3072355), new ComplexNumber({ real: -50485.2048527, imaginary: -330375.0546527 }),
new ComplexNumber(-96389.4195635, 393408.4543872), new ComplexNumber({ real: 122235.7738708, imaginary: 91091.6398019 }),
new ComplexNumber(-173449.0825417, 146875.7724104), new ComplexNumber({ real: 47625.8850387, imaginary: 73497.3981523 }),
new ComplexNumber(-179002.5662573, 239821.0124341)]; new ComplexNumber({ real: -15619.8231136, imaginary: 80804.8685410 }),
new ComplexNumber({ real: 192234.0276101, imaginary: 160833.3072355 }),
new ComplexNumber({ real: -96389.4195635, imaginary: 393408.4543872 }),
new ComplexNumber({ real: -173449.0825417, imaginary: 146875.7724104 }),
new ComplexNumber({ real: -179002.5662573, imaginary: 239821.0124341 }),
];
const out3 = fastFourierTransform(in3); const out3 = fastFourierTransform(in3);
const invOut3 = fastFourierTransform(out3, true); const invOut3 = fastFourierTransform(out3, true);
expect(approximatelyEqual(expOut3, out3, eps)).toBe(true); expect(approximatelyEqual(expOut3, out3, eps)).toBe(true);

View File

@ -1,36 +0,0 @@
export default class ComplexNumber {
/**
* @param {Number} [real]
* @param {Number} [imaginary]
*/
constructor(real, imaginary) {
this.real = real;
this.imaginary = imaginary;
}
/**
* @param {ComplexNumber} [addend]
* @return {ComplexNumber}
*/
add(addend) {
return new ComplexNumber(this.real + addend.real, this.imaginary + addend.imaginary);
}
/**
* @param {ComplexNumber} [subtrahend]
* @return {ComplexNumber}
*/
subtract(subtrahend) {
return new ComplexNumber(this.real - subtrahend.real, this.imaginary - subtrahend.imaginary);
}
/**
* @param {ComplexNumber} [multiplicand]
* @return {ComplexNumber}
*/
multiply(multiplicand) {
const real = this.real * multiplicand.real - this.imaginary * multiplicand.imaginary;
const imaginary = this.real * multiplicand.imaginary + this.imaginary * multiplicand.real;
return new ComplexNumber(real, imaginary);
}
}

View File

@ -1,7 +1,8 @@
import ComplexNumber from './complex'; import ComplexNumber from '../complex-number/ComplexNumber';
/** /**
* Return the no of bits used in the binary representation of input * Return the no of bits used in the binary representation of input.
*
* @param {Number} [input] * @param {Number} [input]
* @return {Number} * @return {Number}
*/ */
@ -14,7 +15,8 @@ function bitLength(input) {
} }
/** /**
* Returns the number which is the flipped binary representation of input * Returns the number which is the flipped binary representation of input.
*
* @param {Number} [input] * @param {Number} [input]
* @param {Number} [bitlen] * @param {Number} [bitlen]
* @return {Number} * @return {Number}
@ -29,8 +31,9 @@ function reverseBits(input, bitlen) {
} }
/** /**
* Returns the radix-2 fast fourier transform of the given array * Returns the radix-2 fast fourier transform of the given array.
* Optionally computes the radix-2 inverse fast fourier transform * Optionally computes the radix-2 inverse fast fourier transform.
*
* @param {ComplexNumber[]} [inputData] * @param {ComplexNumber[]} [inputData]
* @param {Boolean} [inverse] * @param {Boolean} [inverse]
* @return {ComplexNumber[]} * @return {ComplexNumber[]}
@ -39,8 +42,12 @@ export default function fastFourierTransform(inputData, inverse = false) {
const bitlen = bitLength(inputData.length - 1); const bitlen = bitLength(inputData.length - 1);
const N = 1 << bitlen; const N = 1 << bitlen;
while (inputData.length < N) { inputData.push(new ComplexNumber(0, 0)); } while (inputData.length < N) {
inputData.push(new ComplexNumber({
real: 0,
imaginary: 0,
}));
}
const output = []; const output = [];
for (let i = 0; i < N; i += 1) { output[i] = inputData[reverseBits(i, bitlen)]; } for (let i = 0; i < N; i += 1) { output[i] = inputData[reverseBits(i, bitlen)]; }
@ -48,15 +55,22 @@ export default function fastFourierTransform(inputData, inverse = false) {
for (let blockLength = 2; blockLength <= N; blockLength *= 2) { for (let blockLength = 2; blockLength <= N; blockLength *= 2) {
let phaseStep; let phaseStep;
if (inverse) { if (inverse) {
phaseStep = new ComplexNumber(Math.cos(2 * Math.PI / blockLength), phaseStep = new ComplexNumber({
-1 * Math.sin(2 * Math.PI / blockLength)); real: Math.cos(2 * Math.PI / blockLength),
imaginary: -1 * Math.sin(2 * Math.PI / blockLength),
});
} else { } else {
phaseStep = new ComplexNumber(Math.cos(2 * Math.PI / blockLength), phaseStep = new ComplexNumber({
Math.sin(2 * Math.PI / blockLength)); real: Math.cos(2 * Math.PI / blockLength),
imaginary: Math.sin(2 * Math.PI / blockLength),
});
} }
for (let blockStart = 0; blockStart < N; blockStart += blockLength) { for (let blockStart = 0; blockStart < N; blockStart += blockLength) {
let phase = new ComplexNumber(1, 0); let phase = new ComplexNumber({
real: 1,
imaginary: 0,
});
for (let idx = blockStart; idx < blockStart + blockLength / 2; idx += 1) { for (let idx = blockStart; idx < blockStart + blockLength / 2; idx += 1) {
const upd1 = output[idx].add(output[idx + blockLength / 2].multiply(phase)); const upd1 = output[idx].add(output[idx + blockLength / 2].multiply(phase));