From 6e63a0548e02a7dbd54a11150d0c9bcf0bdce7c3 Mon Sep 17 00:00:00 2001 From: Oleksii Trekhleb Date: Thu, 16 Aug 2018 13:14:40 +0300 Subject: [PATCH] Add Inverse Discrete Fourier Transform. --- README.md | 2 +- .../__test__/FourierTester.js | 51 +++++++++++----- .../inverseDiscreteFourierTransform.test.js | 8 +++ .../discreteFourierTransform.js | 15 +++-- .../inverseDiscreteFourierTransform.js | 58 +++++++++++++++++++ 5 files changed, 112 insertions(+), 22 deletions(-) create mode 100644 src/algorithms/math/fourier-transform/__test__/inverseDiscreteFourierTransform.test.js create mode 100644 src/algorithms/math/fourier-transform/inverseDiscreteFourierTransform.js diff --git a/README.md b/README.md index edac6451..1b7df70f 100644 --- a/README.md +++ b/README.md @@ -73,7 +73,7 @@ a set of rules that precisely define a sequence of operations. * `B` [Radian & Degree](src/algorithms/math/radian) - radians to degree and backwards conversion * `A` [Integer Partition](src/algorithms/math/integer-partition) * `A` [Liu Hui π Algorithm](src/algorithms/math/liu-hui) - approximate π calculations based on N-gons - * `A` [Fourier Transform (DFT, FFT)](src/algorithms/math/fourier-transform) - decompose a function of time (a signal) into the frequencies that make it up + * `A` [Fourier Transform (DFT, IDFT, FFT)](src/algorithms/math/fourier-transform) - decompose a function of time (a signal) into the frequencies that make it up * **Sets** * `B` [Cartesian Product](src/algorithms/sets/cartesian-product) - product of multiple sets * `B` [Fisher–Yates Shuffle](src/algorithms/sets/fisher-yates) - random permutation of a finite sequence diff --git a/src/algorithms/math/fourier-transform/__test__/FourierTester.js b/src/algorithms/math/fourier-transform/__test__/FourierTester.js index 681d7668..84a16bf2 100644 --- a/src/algorithms/math/fourier-transform/__test__/FourierTester.js +++ b/src/algorithms/math/fourier-transform/__test__/FourierTester.js @@ -1,6 +1,6 @@ import ComplexNumber from '../../complex-number/ComplexNumber'; -export const fourierDirectTestCases = [ +export const fourierTestCases = [ { input: [ { amplitude: 1 }, @@ -47,13 +47,13 @@ export const fourierDirectTestCases = [ ], output: [ { - frequency: 0, amplitude: 0.3333, phase: 0, re: 0.3333, im: 0, + frequency: 0, amplitude: 0.33333, phase: 0, re: 0.33333, im: 0, }, { - frequency: 1, amplitude: 0.3333, phase: 0, re: 0.3333, im: 0, + frequency: 1, amplitude: 0.33333, phase: 0, re: 0.33333, im: 0, }, { - frequency: 2, amplitude: 0.3333, phase: 0, re: 0.3333, im: 0, + frequency: 2, amplitude: 0.33333, phase: 0, re: 0.33333, im: 0, }, ], }, @@ -179,13 +179,13 @@ export const fourierDirectTestCases = [ frequency: 0, amplitude: 1.75, phase: 0, re: 1.75, im: 0, }, { - frequency: 1, amplitude: 1.03077, phase: 14.0362, re: 0.99999, im: 0.25, + frequency: 1, amplitude: 1.03077, phase: 14.03624, re: 0.99999, im: 0.25, }, { frequency: 2, amplitude: 0.25, phase: 0, re: 0.25, im: 0, }, { - frequency: 3, amplitude: 1.03077, phase: -14.0362, re: 1, im: -0.25, + frequency: 3, amplitude: 1.03077, phase: -14.03624, re: 1, im: -0.25, }, ], }, @@ -201,7 +201,7 @@ export const fourierDirectTestCases = [ frequency: 0, amplitude: 1, phase: 0, re: 1, im: 0, }, { - frequency: 1, amplitude: 1.76776, phase: 8.1301, re: 1.75, im: 0.25, + frequency: 1, amplitude: 1.76776, phase: 8.13010, re: 1.75, im: 0.25, }, { frequency: 2, amplitude: 0.5, phase: 180, re: -0.5, im: 0, @@ -240,17 +240,15 @@ export default class FourierTester { * @param {function} fourierTransform */ static testDirectFourierTransform(fourierTransform) { - fourierDirectTestCases.forEach((testCase) => { + fourierTestCases.forEach((testCase) => { const { input, output: expectedOutput } = testCase; - // Convert input into complex numbers. - const complexInput = input.map(sample => new ComplexNumber({ re: sample.amplitude })); - // Try to split input signal into sequence of pure sinusoids. - const currentOutput = fourierTransform(complexInput); + const formattedInput = input.map(sample => sample.amplitude); + const currentOutput = fourierTransform(formattedInput); // Check the signal has been split into proper amount of sub-signals. - expect(currentOutput.length).toBeGreaterThanOrEqual(complexInput.length); + expect(currentOutput.length).toBeGreaterThanOrEqual(formattedInput.length); // Now go through all the signals and check their frequency, amplitude and phase. expectedOutput.forEach((expectedSignal, frequency) => { @@ -267,4 +265,31 @@ export default class FourierTester { }); }); } + + /** + * @param {function} inverseFourierTransform + */ + static testInverseFourierTransform(inverseFourierTransform) { + fourierTestCases.forEach((testCase) => { + const { input: expectedOutput, output: inputFrequencies } = testCase; + + // Try to join frequencies into time signal. + const formattedInput = inputFrequencies.map((frequency) => { + return new ComplexNumber({ re: frequency.re, im: frequency.im }); + }); + const currentOutput = inverseFourierTransform(formattedInput); + + // Check the signal has been combined of proper amount of time samples. + expect(currentOutput.length).toBeLessThanOrEqual(formattedInput.length); + + // Now go through all the amplitudes and check their values. + expectedOutput.forEach((expectedAmplitudes, timer) => { + // Get template data we want to test against. + const currentAmplitude = currentOutput[timer]; + + // Check if current amplitude is close enough to the calculated one. + expect(currentAmplitude).toBeCloseTo(expectedAmplitudes.amplitude, 4); + }); + }); + } } diff --git a/src/algorithms/math/fourier-transform/__test__/inverseDiscreteFourierTransform.test.js b/src/algorithms/math/fourier-transform/__test__/inverseDiscreteFourierTransform.test.js new file mode 100644 index 00000000..3919b3ad --- /dev/null +++ b/src/algorithms/math/fourier-transform/__test__/inverseDiscreteFourierTransform.test.js @@ -0,0 +1,8 @@ +import inverseDiscreteFourierTransform from '../inverseDiscreteFourierTransform'; +import FourierTester from './FourierTester'; + +describe('inverseDiscreteFourierTransform', () => { + it('should calculate output signal out of input frequencies', () => { + FourierTester.testInverseFourierTransform(inverseDiscreteFourierTransform); + }); +}); diff --git a/src/algorithms/math/fourier-transform/discreteFourierTransform.js b/src/algorithms/math/fourier-transform/discreteFourierTransform.js index 0fb4666a..dcb0a014 100644 --- a/src/algorithms/math/fourier-transform/discreteFourierTransform.js +++ b/src/algorithms/math/fourier-transform/discreteFourierTransform.js @@ -3,12 +3,14 @@ import ComplexNumber from '../complex-number/ComplexNumber'; const CLOSE_TO_ZERO_THRESHOLD = 1e-10; /** - * Discrete Fourier Transform. + * Discrete Fourier Transform (DFT): time to frequencies. * * Time complexity: O(N^2) * - * @param {ComplexNumber[]} complexInputAmplitudes - Input signal amplitudes over time (complex + * @param {number[]} inputAmplitudes - Input signal amplitudes over time (complex * numbers with real parts only). + * @param {number} zeroThreshold - Threshold that is used to convert real and imaginary numbers + * to zero in case if they are smaller then this. * * @return {ComplexNumber[]} - Array of complex number. Each of the number represents the frequency * or signal. All signals together will form input signal over discrete time periods. Each signal's @@ -17,10 +19,7 @@ const CLOSE_TO_ZERO_THRESHOLD = 1e-10; * @see https://gist.github.com/anonymous/129d477ddb1c8025c9ac * @see https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/ */ -export default function dft(complexInputAmplitudes) { - // Convert complex amplitudes into real ones. - const inputAmplitudes = complexInputAmplitudes.map(complexAmplitude => complexAmplitude.re); - +export default function dft(inputAmplitudes, zeroThreshold = CLOSE_TO_ZERO_THRESHOLD) { const N = inputAmplitudes.length; const signals = []; @@ -48,11 +47,11 @@ export default function dft(complexInputAmplitudes) { } // Close to zero? You're zero. - if (Math.abs(frequencySignal.re) < CLOSE_TO_ZERO_THRESHOLD) { + if (Math.abs(frequencySignal.re) < zeroThreshold) { frequencySignal.re = 0; } - if (Math.abs(frequencySignal.im) < CLOSE_TO_ZERO_THRESHOLD) { + if (Math.abs(frequencySignal.im) < zeroThreshold) { frequencySignal.im = 0; } diff --git a/src/algorithms/math/fourier-transform/inverseDiscreteFourierTransform.js b/src/algorithms/math/fourier-transform/inverseDiscreteFourierTransform.js new file mode 100644 index 00000000..a9f06163 --- /dev/null +++ b/src/algorithms/math/fourier-transform/inverseDiscreteFourierTransform.js @@ -0,0 +1,58 @@ +import ComplexNumber from '../complex-number/ComplexNumber'; + +const CLOSE_TO_ZERO_THRESHOLD = 1e-10; + +/** + * Inverse Discrete Fourier Transform (IDFT): frequencies to time. + * + * Time complexity: O(N^2) + * + * @param {ComplexNumber[]} frequencies - Frequencies summands of the final signal. + * @param {number} zeroThreshold - Threshold that is used to convert real and imaginary numbers + * to zero in case if they are smaller then this. + * + * @return {number[]} - Discrete amplitudes distributed in time. + */ +export default function inverseDiscreteFourierTransform( + frequencies, + zeroThreshold = CLOSE_TO_ZERO_THRESHOLD, +) { + const N = frequencies.length; + const amplitudes = []; + + // Go through every discrete point of time. + for (let timer = 0; timer < N; timer += 1) { + // Compound amplitude at current time. + let amplitude = new ComplexNumber(); + + // Go through all discrete frequencies. + for (let frequency = 0; frequency < N; frequency += 1) { + const currentFrequency = frequencies[frequency]; + + // Calculate rotation angle. + const rotationAngle = (2 * Math.PI) * frequency * (timer / N); + + // Remember that e^ix = cos(x) + i * sin(x); + const frequencyContribution = new ComplexNumber({ + re: Math.cos(rotationAngle), + im: Math.sin(rotationAngle), + }).multiply(currentFrequency); + + amplitude = amplitude.add(frequencyContribution); + } + + // Close to zero? You're zero. + if (Math.abs(amplitude.re) < zeroThreshold) { + amplitude.re = 0; + } + + if (Math.abs(amplitude.im) < zeroThreshold) { + amplitude.im = 0; + } + + // Add current frequency signal to the list of compound signals. + amplitudes[timer] = amplitude.re; + } + + return amplitudes; +}