Refactor DFT and add common tests for Fourier.

This commit is contained in:
Oleksii Trekhleb 2018-08-16 12:37:06 +03:00
parent 13ed5061a3
commit 351a745f55
3 changed files with 314 additions and 164 deletions

View File

@ -0,0 +1,270 @@
import ComplexNumber from '../../complex-number/ComplexNumber';
export const fourierDirectTestCases = [
{
input: [
{ amplitude: 1 },
],
output: [
{
frequency: 0, amplitude: 1, phase: 0, re: 1, im: 0,
},
],
},
{
input: [
{ amplitude: 1 },
{ amplitude: 0 },
],
output: [
{
frequency: 0, amplitude: 0.5, phase: 0, re: 0.5, im: 0,
},
{
frequency: 1, amplitude: 0.5, phase: 0, re: 0.5, im: 0,
},
],
},
{
input: [
{ amplitude: 2 },
{ amplitude: 0 },
],
output: [
{
frequency: 0, amplitude: 1, phase: 0, re: 1, im: 0,
},
{
frequency: 1, amplitude: 1, phase: 0, re: 1, im: 0,
},
],
},
{
input: [
{ amplitude: 1 },
{ amplitude: 0 },
{ amplitude: 0 },
],
output: [
{
frequency: 0, amplitude: 0.3333, phase: 0, re: 0.3333, im: 0,
},
{
frequency: 1, amplitude: 0.3333, phase: 0, re: 0.3333, im: 0,
},
{
frequency: 2, amplitude: 0.3333, phase: 0, re: 0.3333, im: 0,
},
],
},
{
input: [
{ amplitude: 1 },
{ amplitude: 0 },
{ amplitude: 0 },
{ amplitude: 0 },
],
output: [
{
frequency: 0, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 1, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 2, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 3, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
],
},
{
input: [
{ amplitude: 0 },
{ amplitude: 1 },
{ amplitude: 0 },
{ amplitude: 0 },
],
output: [
{
frequency: 0, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 1, amplitude: 0.25, phase: -90, re: 0, im: -0.25,
},
{
frequency: 2, amplitude: 0.25, phase: 180, re: -0.25, im: 0,
},
{
frequency: 3, amplitude: 0.25, phase: 90, re: 0, im: 0.25,
},
],
},
{
input: [
{ amplitude: 0 },
{ amplitude: 0 },
{ amplitude: 1 },
{ amplitude: 0 },
],
output: [
{
frequency: 0, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 1, amplitude: 0.25, phase: 180, re: -0.25, im: 0,
},
{
frequency: 2, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 3, amplitude: 0.25, phase: 180, re: -0.25, im: 0,
},
],
},
{
input: [
{ amplitude: 0 },
{ amplitude: 0 },
{ amplitude: 0 },
{ amplitude: 2 },
],
output: [
{
frequency: 0, amplitude: 0.5, phase: 0, re: 0.5, im: 0,
},
{
frequency: 1, amplitude: 0.5, phase: 90, re: 0, im: 0.5,
},
{
frequency: 2, amplitude: 0.5, phase: 180, re: -0.5, im: 0,
},
{
frequency: 3, amplitude: 0.5, phase: -90, re: 0, im: -0.5,
},
],
},
{
input: [
{ amplitude: 0 },
{ amplitude: 1 },
{ amplitude: 0 },
{ amplitude: 2 },
],
output: [
{
frequency: 0, amplitude: 0.75, phase: 0, re: 0.75, im: 0,
},
{
frequency: 1, amplitude: 0.25, phase: 90, re: 0, im: 0.25,
},
{
frequency: 2, amplitude: 0.75, phase: 180, re: -0.75, im: 0,
},
{
frequency: 3, amplitude: 0.25, phase: -90, re: 0, im: -0.25,
},
],
},
{
input: [
{ amplitude: 4 },
{ amplitude: 1 },
{ amplitude: 0 },
{ amplitude: 2 },
],
output: [
{
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: 2, amplitude: 0.25, phase: 0, re: 0.25, im: 0,
},
{
frequency: 3, amplitude: 1.03077, phase: -14.0362, re: 1, im: -0.25,
},
],
},
{
input: [
{ amplitude: 4 },
{ amplitude: 1 },
{ amplitude: -3 },
{ amplitude: 2 },
],
output: [
{
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: 2, amplitude: 0.5, phase: 180, re: -0.5, im: 0,
},
{
frequency: 3, amplitude: 1.76776, phase: -8.13010, re: 1.75, im: -0.24999,
},
],
},
{
input: [
{ amplitude: 1 },
{ amplitude: 2 },
{ amplitude: 3 },
{ amplitude: 4 },
],
output: [
{
frequency: 0, amplitude: 2.5, phase: 0, re: 2.5, im: 0,
},
{
frequency: 1, amplitude: 0.70710, phase: 135, re: -0.5, im: 0.49999,
},
{
frequency: 2, amplitude: 0.5, phase: 180, re: -0.5, im: 0,
},
{
frequency: 3, amplitude: 0.70710, phase: -134.99999, re: -0.49999, im: -0.5,
},
],
},
];
export default class FourierTester {
/**
* @param {function} fourierTransform
*/
static testDirectFourierTransform(fourierTransform) {
fourierDirectTestCases.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);
// Check the signal has been split into proper amount of sub-signals.
expect(currentOutput.length).toBeGreaterThanOrEqual(complexInput.length);
// Now go through all the signals and check their frequency, amplitude and phase.
expectedOutput.forEach((expectedSignal, frequency) => {
// Get template data we want to test against.
const currentSignal = currentOutput[frequency];
const currentPolarSignal = currentSignal.getPolarForm(false);
// Check all signal parameters.
expect(frequency).toBe(expectedSignal.frequency);
expect(currentSignal.re).toBeCloseTo(expectedSignal.re, 4);
expect(currentSignal.im).toBeCloseTo(expectedSignal.im, 4);
expect(currentPolarSignal.phase).toBeCloseTo(expectedSignal.phase, 4);
expect(currentPolarSignal.radius).toBeCloseTo(expectedSignal.amplitude, 4);
});
});
}
}

View File

@ -1,143 +1,8 @@
import discreteFourierTransform from '../discreteFourierTransform';
/**
* Helper class of the output Signal.
*/
class Sgnl {
constructor(frequency, amplitude, phase) {
this.frequency = frequency;
this.amplitude = amplitude;
this.phase = phase;
}
}
import FourierTester from './FourierTester';
describe('discreteFourierTransform', () => {
it('should split signal into frequencies', () => {
const testCases = [
{
inputAmplitudes: [1],
outputSignals: [
new Sgnl(0, 1, 0),
],
},
{
inputAmplitudes: [1, 0],
outputSignals: [
new Sgnl(0, 0.5, 0),
new Sgnl(1, 0.5, 0),
],
},
{
inputAmplitudes: [2, 0],
outputSignals: [
new Sgnl(0, 1, 0),
new Sgnl(1, 1, 0),
],
},
{
inputAmplitudes: [1, 0, 0],
outputSignals: [
new Sgnl(0, 0.33, 0),
new Sgnl(1, 0.33, 0),
new Sgnl(2, 0.33, 0),
],
},
{
inputAmplitudes: [1, 0, 0, 0],
outputSignals: [
new Sgnl(0, 0.25, 0),
new Sgnl(1, 0.25, 0),
new Sgnl(2, 0.25, 0),
new Sgnl(3, 0.25, 0),
],
},
{
inputAmplitudes: [0, 1, 0, 0],
outputSignals: [
new Sgnl(0, 0.25, 0),
new Sgnl(1, 0.25, -90),
new Sgnl(2, 0.25, 180),
new Sgnl(3, 0.25, 90),
],
},
{
inputAmplitudes: [0, 0, 1, 0],
outputSignals: [
new Sgnl(0, 0.25, 0),
new Sgnl(1, 0.25, 180),
new Sgnl(2, 0.25, 0),
new Sgnl(3, 0.25, 180),
],
},
{
inputAmplitudes: [0, 0, 0, 2],
outputSignals: [
new Sgnl(0, 0.5, 0),
new Sgnl(1, 0.5, 90),
new Sgnl(2, 0.5, 180),
new Sgnl(3, 0.5, -90),
],
},
{
inputAmplitudes: [0, 1, 0, 2],
outputSignals: [
new Sgnl(0, 0.75, 0),
new Sgnl(1, 0.25, 90),
new Sgnl(2, 0.75, 180),
new Sgnl(3, 0.25, -90),
],
},
{
inputAmplitudes: [4, 1, 0, 2],
outputSignals: [
new Sgnl(0, 1.75, 0),
new Sgnl(1, 1.03, 14),
new Sgnl(2, 0.25, 0),
new Sgnl(3, 1.03, -14),
],
},
{
inputAmplitudes: [4, 1, -3, 2],
outputSignals: [
new Sgnl(0, 1, 0),
new Sgnl(1, 1.77, 8),
new Sgnl(2, 0.5, 180),
new Sgnl(3, 1.77, -8),
],
},
{
inputAmplitudes: [1, 2, 3, 4],
outputSignals: [
new Sgnl(0, 2.5, 0),
new Sgnl(1, 0.71, 135),
new Sgnl(2, 0.5, 180),
new Sgnl(3, 0.71, -135),
],
},
];
testCases.forEach((testCase) => {
const { inputAmplitudes, outputSignals } = testCase;
// Try to split input signal into sequence of pure sinusoids.
const signals = discreteFourierTransform(inputAmplitudes);
// Check the signal has been split into proper amount of sub-signals.
expect(signals.length).toBe(outputSignals.length);
// Now go through all the signals and check their frequency, amplitude and phase.
signals.forEach((signal, frequency) => {
// Get polar form of calculated sub-signal since it is more convenient to analyze.
const signalPolarForm = signal.getPolarForm(false);
// Get template data we want to test against.
const signalTemplate = outputSignals[frequency];
// Check all signal parameters.
expect(frequency).toBe(signalTemplate.frequency);
expect(signalPolarForm.radius).toBeCloseTo(signalTemplate.amplitude, 2);
expect(signalPolarForm.phase).toBeCloseTo(signalTemplate.phase, 0);
});
});
FourierTester.testDirectFourierTransform(discreteFourierTransform);
});
});

View File

@ -1,7 +1,15 @@
import ComplexNumber from '../complex-number/ComplexNumber';
const CLOSE_TO_ZERO_THRESHOLD = 1e-10;
/**
* @param {number[]} inputSignalAmplitudes - Input signal amplitudes over time (i.e. [1, 0, 4]).
* Discrete Fourier Transform.
*
* Time complexity: O(N^2)
*
* @param {ComplexNumber[]} complexInputAmplitudes - Input signal amplitudes over time (complex
* numbers with real parts only).
*
* @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
* complex number has radius (amplitude) and phase (angle) in polar form that describes the signal.
@ -9,47 +17,54 @@ import ComplexNumber from '../complex-number/ComplexNumber';
* @see https://gist.github.com/anonymous/129d477ddb1c8025c9ac
* @see https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/
*/
export default function discreteFourierTransform(inputSignalAmplitudes) {
const N = inputSignalAmplitudes.length;
const outputFrequencies = [];
export default function dft(complexInputAmplitudes) {
// Convert complex amplitudes into real ones.
const inputAmplitudes = complexInputAmplitudes.map(complexAmplitude => complexAmplitude.re);
// For every frequency discrete...
for (let frequencyValue = 0; frequencyValue < N; frequencyValue += 1) {
let signal = new ComplexNumber();
const N = inputAmplitudes.length;
const signals = [];
// For every discrete point in time...
for (let t = 0; t < N; t += 1) {
// Spin the signal _backwards_ at each frequency (as radians/s, not Hertz)
const rate = -1 * (2 * Math.PI) * frequencyValue;
// Go through every discrete frequency.
for (let frequency = 0; frequency < N; frequency += 1) {
// Compound signal at current frequency that will ultimately
// take part in forming input amplitudes.
let frequencySignal = new ComplexNumber();
// How far around the circle have we gone at time=t?
const time = t / N;
const distance = rate * time;
// Go through every discrete point in time.
for (let timer = 0; timer < N; timer += 1) {
const currentAmplitude = inputAmplitudes[timer];
// Data-point * e^(-i*2*pi*f) is complex, store each part.
// Calculate rotation angle.
const rotationAngle = -1 * (2 * Math.PI) * frequency * (timer / N);
// Remember that e^ix = cos(x) + i * sin(x);
const dataPointContribution = new ComplexNumber({
re: inputSignalAmplitudes[t] * Math.cos(distance),
im: inputSignalAmplitudes[t] * Math.sin(distance),
});
re: Math.cos(rotationAngle),
im: Math.sin(rotationAngle),
}).multiply(currentAmplitude);
// Add this data point's contribution.
signal = signal.add(dataPointContribution);
frequencySignal = frequencySignal.add(dataPointContribution);
}
// Close to zero? You're zero.
if (Math.abs(signal.re) < 1e-10) {
signal.re = 0;
if (Math.abs(frequencySignal.re) < CLOSE_TO_ZERO_THRESHOLD) {
frequencySignal.re = 0;
}
if (Math.abs(signal.im) < 1e-10) {
signal.im = 0;
if (Math.abs(frequencySignal.im) < CLOSE_TO_ZERO_THRESHOLD) {
frequencySignal.im = 0;
}
// Average contribution at this frequency
signal = signal.divide(N);
// Average contribution at this frequency.
// The 1/N factor is usually moved to the reverse transform (going from frequencies
// back to time). This is allowed, though it would be nice to have 1/N in the forward
// transform since it gives the actual sizes for the time spikes.
frequencySignal = frequencySignal.divide(N);
outputFrequencies[frequencyValue] = signal;
// Add current frequency signal to the list of compound signals.
signals[frequency] = frequencySignal;
}
return outputFrequencies;
return signals;
}