diff --git a/src/algorithms/math/fast-fourier-transform/README.md b/src/algorithms/math/fast-fourier-transform/README.md new file mode 100644 index 00000000..62d42c12 --- /dev/null +++ b/src/algorithms/math/fast-fourier-transform/README.md @@ -0,0 +1,11 @@ +# Discrete Fourier transform +The Discrete Fourier transform transforms a sequence of `N` complex numbers +**{xn}** := **x0, x1, x2 ..., xN-1**   into another sequence of complex numbers
**{Xk}** := **X0, X1, X2 ..., XN-1**    which is defined by + +![alt text](https://wikimedia.org/api/rest_v1/media/math/render/svg/1af0a78dc50bbf118ab6bd4c4dcc3c4ff8502223) + + +## References + +- [Wikipedia, DFT](https://www.wikiwand.com/en/Discrete_Fourier_transform) +- [Wikipedia, FFT](https://www.wikiwand.com/en/Fast_Fourier_transform) diff --git a/src/algorithms/math/fast-fourier-transform/__test__/fastFourierTransform.test.js b/src/algorithms/math/fast-fourier-transform/__test__/fastFourierTransform.test.js new file mode 100644 index 00000000..b5b89e27 --- /dev/null +++ b/src/algorithms/math/fast-fourier-transform/__test__/fastFourierTransform.test.js @@ -0,0 +1,70 @@ +import ComplexNumber from '../complex'; +import fastFourierTransform from '../fastFourierTransform'; +/** + * @param {ComplexNumber[]} [seq1] + * @param {ComplexNumber[]} [seq2] + * @param {Number} [eps] + * @return {boolean} + */ +function approximatelyEqual(seq1, seq2, eps) { + if (seq1.length !== seq2.length) { return false; } + + for (let i = 0; i < seq1.length; i += 1) { + if (Math.abs(seq1[i].real - seq2[i].real) > eps) { return false; } + if (Math.abs(seq1[i].complex - seq2[i].complex) > eps) { return false; } + } + + return true; +} + +describe('fastFourierTransform', () => { + it('should calculate the radix-2 discrete fourier transform after zero padding', () => { + const eps = 1e-6; + const in1 = [new ComplexNumber(0, 0)]; + const expOut1 = [new ComplexNumber(0, 0)]; + const out1 = fastFourierTransform(in1); + const invOut1 = fastFourierTransform(out1, true); + expect(approximatelyEqual(expOut1, out1, eps)).toBe(true); + expect(approximatelyEqual(in1, invOut1, eps)).toBe(true); + + const in2 = [new ComplexNumber(1, 2), new ComplexNumber(2, 3), + new ComplexNumber(8, 4)]; + const expOut2 = [new ComplexNumber(11, 9), new ComplexNumber(-10, 0), + new ComplexNumber(7, 3), new ComplexNumber(-4, -4)]; + const out2 = fastFourierTransform(in2); + const invOut2 = fastFourierTransform(out2, true); + expect(approximatelyEqual(expOut2, out2, eps)).toBe(true); + expect(approximatelyEqual(in2, invOut2, eps)).toBe(true); + + const in3 = [new ComplexNumber(-83656.9359385182, 98724.08038374918), + new ComplexNumber(-47537.415125808424, 88441.58381765135), + new ComplexNumber(-24849.657029355192, -72621.79007878687), + new ComplexNumber(31451.27290052717, -21113.301128347346), + new ComplexNumber(13973.90836288876, -73378.36721594246), + new ComplexNumber(14981.520420492234, 63279.524958963884), + new ComplexNumber(-9892.575367044381, -81748.44671677813), + new ComplexNumber(-35933.00356823792, -46153.47157161784), + new ComplexNumber(-22425.008561855735, -86284.24507370662), + new ComplexNumber(-39327.43830818355, 30611.949874562706)]; + const expOut3 = [new ComplexNumber(-203215.3322151, -100242.4827503), + new ComplexNumber(99217.0805705, 270646.9331932), + new ComplexNumber(-305990.9040412, 68224.8435751), + new ComplexNumber(-14135.7758282, 199223.9878095), + new ComplexNumber(-306965.6350922, 26030.1025439), + new ComplexNumber(-76477.6755206, 40781.9078990), + new ComplexNumber(-48409.3099088, 54674.7959662), + new ComplexNumber(-329683.0131713, 164287.7995937), + new ComplexNumber(-50485.2048527, -330375.0546527), + new ComplexNumber(122235.7738708, 91091.6398019), + new ComplexNumber(47625.8850387, 73497.3981523), + new ComplexNumber(-15619.8231136, 80804.8685410), + new ComplexNumber(192234.0276101, 160833.3072355), + new ComplexNumber(-96389.4195635, 393408.4543872), + new ComplexNumber(-173449.0825417, 146875.7724104), + new ComplexNumber(-179002.5662573, 239821.0124341)]; + const out3 = fastFourierTransform(in3); + const invOut3 = fastFourierTransform(out3, true); + expect(approximatelyEqual(expOut3, out3, eps)).toBe(true); + expect(approximatelyEqual(in3, invOut3, eps)).toBe(true); + }); +}); diff --git a/src/algorithms/math/fast-fourier-transform/complex.js b/src/algorithms/math/fast-fourier-transform/complex.js new file mode 100644 index 00000000..7a6279e6 --- /dev/null +++ b/src/algorithms/math/fast-fourier-transform/complex.js @@ -0,0 +1,36 @@ +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); + } +} diff --git a/src/algorithms/math/fast-fourier-transform/fastFourierTransform.js b/src/algorithms/math/fast-fourier-transform/fastFourierTransform.js new file mode 100644 index 00000000..3869d0d9 --- /dev/null +++ b/src/algorithms/math/fast-fourier-transform/fastFourierTransform.js @@ -0,0 +1,76 @@ +import ComplexNumber from './complex'; + +/** + * Return the no of bits used in the binary representation of input + * @param {Number} [input] + * @return {Number} + */ +function bitLength(input) { + let bitlen = 0; + while ((1 << bitlen) <= input) { + bitlen += 1; + } + return bitlen; +} + +/** + * Returns the number which is the flipped binary representation of input + * @param {Number} [input] + * @param {Number} [bitlen] + * @return {Number} + */ +function reverseBits(input, bitlen) { + let reversedBits = 0; + for (let i = 0; i < bitlen; i += 1) { + reversedBits *= 2; + if (Math.floor(input / (1 << i)) % 2 === 1) { reversedBits += 1; } + } + return reversedBits; +} + +/** + * Returns the radix-2 fast fourier transform of the given array + * Optionally computes the radix-2 inverse fast fourier transform + * @param {ComplexNumber[]} [inputData] + * @param {Boolean} [inverse] + * @return {ComplexNumber[]} + */ +export default function fastFourierTransform(inputData, inverse = false) { + const bitlen = bitLength(inputData.length - 1); + const N = 1 << bitlen; + + while (inputData.length < N) { inputData.push(new ComplexNumber(0, 0)); } + + + const output = []; + for (let i = 0; i < N; i += 1) { output[i] = inputData[reverseBits(i, bitlen)]; } + + for (let blockLength = 2; blockLength <= N; blockLength *= 2) { + let phaseStep; + if (inverse) { + phaseStep = new ComplexNumber(Math.cos(2 * Math.PI / blockLength), + -1 * Math.sin(2 * Math.PI / blockLength)); + } else { + phaseStep = new ComplexNumber(Math.cos(2 * Math.PI / blockLength), + Math.sin(2 * Math.PI / blockLength)); + } + + for (let blockStart = 0; blockStart < N; blockStart += blockLength) { + let phase = new ComplexNumber(1, 0); + + for (let idx = blockStart; idx < blockStart + blockLength / 2; idx += 1) { + const upd1 = output[idx].add(output[idx + blockLength / 2].multiply(phase)); + const upd2 = output[idx].subtract(output[idx + blockLength / 2].multiply(phase)); + output[idx] = upd1; + output[idx + blockLength / 2] = upd2; + phase = phase.multiply(phaseStep); + } + } + } + if (inverse) { + for (let idx = 0; idx < N; idx += 1) { + output[idx] /= N; + } + } + return output; +}