Add Inverse Discrete Fourier Transform.

This commit is contained in:
Oleksii Trekhleb 2018-08-16 13:14:40 +03:00
parent 351a745f55
commit 6e63a0548e
5 changed files with 112 additions and 22 deletions

View File

@ -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` [FisherYates Shuffle](src/algorithms/sets/fisher-yates) - random permutation of a finite sequence

View File

@ -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);
});
});
}
}

View File

@ -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);
});
});

View File

@ -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;
}

View File

@ -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;
}