algorithm/javascript-algorithms
1
0
Fork 0
mirror of https://github.com/trekhleb/javascript-algorithms.git synced 2026-03-13 08:51:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add DFT.

Browse Source
This commit is contained in:
Oleksii Trekhleb
2018-08-14 23:12:17 +03:00
parent 8e66189eff
commit d0c4baf7b3
6 changed files with 149 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
src/algorithms/math/fourier-transform/README.md Normal file
View File

@@ -0,0 +1,89 @@
# Fourier Transform
## Definitions
The **Fourier Transform** (**FT**) decomposes a function of time (a signal) into
the frequencies that make it up, in a way similar to how a musical chord can be
expressed as the frequencies (or pitches) of its constituent notes.
The **Discrete Fourier Transform** (**DFT**) converts a finite sequence of
equally-spaced samples of a function into a same-length sequence of
equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a
complex-valued function of frequency. The interval at which the DTFT is sampled
is the reciprocal of the duration of the input sequence. An inverse DFT is a
Fourier series, using the DTFT samples as coefficients of complex sinusoids at
the corresponding DTFT frequencies. It has the same sample-values as the original
input sequence. The DFT is therefore said to be a frequency domain representation
of the original input sequence. If the original sequence spans all the non-zero
values of a function, its DTFT is continuous (and periodic), and the DFT provides
discrete samples of one cycle. If the original sequence is one cycle of a periodic
function, the DFT provides all the non-zero values of one DTFT cycle.
The Discrete Fourier transform transforms a sequence of `N` complex numbers:
{x<sub>n</sub>} = x<sub>0</sub>, x<sub>1</sub>, x<sub>2</sub> ..., x<sub>N-1</sub>
into another sequence of complex numbers:
{X<sub>k</sub>} = X<sub>0</sub>, X<sub>1</sub>, X<sub>2</sub> ..., X<sub>N-1</sub>
which is defined by:
![DFT](https://wikimedia.org/api/rest_v1/media/math/render/svg/1af0a78dc50bbf118ab6bd4c4dcc3c4ff8502223)
The **Discrete-Time Fourier Transform** (**DTFT**) is a form of Fourier analysis
that is applicable to the uniformly-spaced samples of a continuous function. The
term discrete-time refers to the fact that the transform operates on discrete data
(samples) whose interval often has units of time. From only the samples, it
produces a function of frequency that is a periodic summation of the continuous
Fourier transform of the original continuous function.
A **Fast Fourier Transform** (**FFT**) is an algorithm that samples a signal over
a period of time (or space) and divides it into its frequency components. These
components are single sinusoidal oscillations at distinct frequencies each with
their own amplitude and phase.
This transformation is illustrated in Diagram below. Over the time period measured
in the diagram, the signal contains 3 distinct dominant frequencies.
View of a signal in the time and frequency domain:
![FFT](https://upload.wikimedia.org/wikipedia/commons/6/61/FFT-Time-Frequency-View.png)
An FFT algorithm computes the discrete Fourier transform (DFT) of a sequence, or
its inverse (IFFT). Fourier analysis converts a signal from its original domain
to a representation in the frequency domain and vice versa. An FFT rapidly
computes such transformations by factorizing the DFT matrix into a product of
sparse (mostly zero) factors. As a result, it manages to reduce the complexity of
computing the DFT from O(n<sup>2</sup>), which arises if one simply applies the
definition of DFT, to O(n log n), where n is the data size.
Here a discrete Fourier analysis of a sum of cosine waves at 10, 20, 30, 40,
and 50 Hz:
![FFT](https://upload.wikimedia.org/wikipedia/commons/6/64/FFT_of_Cosine_Summation_Function.png)
## Explanation
The Fourier Transform is one of deepest insights ever made. Unfortunately, the
meaning is buried within dense equations:
![](https://betterexplained.com/wp-content/plugins/wp-latexrender/pictures/45c088dbb767150fc0bacfeb49dd49e5.png)
![](https://betterexplained.com/wp-content/plugins/wp-latexrender/pictures/faeb9c5bf2e60add63ae4a70b293c7b4.png)
Rather than jumping into the symbols, let's experience the key idea firsthand. Here's a plain-English metaphor:
- *What does the Fourier Transform do?* Given a smoothie, it finds the recipe.
- *How?* Run the smoothie through filters to extract each ingredient.
- *Why?* Recipes are easier to analyze, compare, and modify than the smoothie itself.
- *How do we get the smoothie back?* Blend the ingredients.
## References
- [An Interactive Guide To The Fourier Transform](https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/)
- [YouTube by Better Explained](https://www.youtube.com/watch?v=iN0VG9N2q0U&t=0s&index=77&list=PLLXdhg_r2hKA7DPDsunoDZ-Z769jWn4R8)
- [YouTube by 3Blue1Brown](https://www.youtube.com/watch?v=spUNpyF58BY&t=0s&index=76&list=PLLXdhg_r2hKA7DPDsunoDZ-Z769jWn4R8)
- [Wikipedia, FT](https://en.wikipedia.org/wiki/Fourier_transform)
- [Wikipedia, DFT](https://www.wikiwand.com/en/Discrete_Fourier_transform)
- [Wikipedia, DTFT](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform)
- [Wikipedia, FFT](https://www.wikiwand.com/en/Fast_Fourier_transform)

9
src/algorithms/math/fourier-transform/__test__/discreteFourierTransform.test.js Normal file
View File

@@ -0,0 +1,9 @@
import discreteFourierTransform from '../discreteFourierTransform';
describe('discreteFourierTransform', () => {
it('should calculate split signal into frequencies', () => {
const frequencies = discreteFourierTransform([1, 0, 0, 0]);
expect(frequencies).toBeDefined();
});
});

85
src/algorithms/math/fourier-transform/__test__/fastFourierTransform.test.js Normal file
View File

@@ -0,0 +1,85 @@
import fastFourierTransform from '../fastFourierTransform';
import ComplexNumber from '../../complex-number/ComplexNumber';
/**
* @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({ re: 0, im: 0 })];
const expOut1 = [new ComplexNumber({ re: 0, im: 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({ re: 1, im: 2 }),
new ComplexNumber({ re: 2, im: 3 }),
new ComplexNumber({ re: 8, im: 4 }),
];
const expOut2 = [
new ComplexNumber({ re: 11, im: 9 }),
new ComplexNumber({ re: -10, im: 0 }),
new ComplexNumber({ re: 7, im: 3 }),
new ComplexNumber({ re: -4, im: -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({ re: -83656.9359385182, im: 98724.08038374918 }),
new ComplexNumber({ re: -47537.415125808424, im: 88441.58381765135 }),
new ComplexNumber({ re: -24849.657029355192, im: -72621.79007878687 }),
new ComplexNumber({ re: 31451.27290052717, im: -21113.301128347346 }),
new ComplexNumber({ re: 13973.90836288876, im: -73378.36721594246 }),
new ComplexNumber({ re: 14981.520420492234, im: 63279.524958963884 }),
new ComplexNumber({ re: -9892.575367044381, im: -81748.44671677813 }),
new ComplexNumber({ re: -35933.00356823792, im: -46153.47157161784 }),
new ComplexNumber({ re: -22425.008561855735, im: -86284.24507370662 }),
new ComplexNumber({ re: -39327.43830818355, im: 30611.949874562706 }),
];
const expOut3 = [
new ComplexNumber({ re: -203215.3322151, im: -100242.4827503 }),
new ComplexNumber({ re: 99217.0805705, im: 270646.9331932 }),
new ComplexNumber({ re: -305990.9040412, im: 68224.8435751 }),
new ComplexNumber({ re: -14135.7758282, im: 199223.9878095 }),
new ComplexNumber({ re: -306965.6350922, im: 26030.1025439 }),
new ComplexNumber({ re: -76477.6755206, im: 40781.9078990 }),
new ComplexNumber({ re: -48409.3099088, im: 54674.7959662 }),
new ComplexNumber({ re: -329683.0131713, im: 164287.7995937 }),
new ComplexNumber({ re: -50485.2048527, im: -330375.0546527 }),
new ComplexNumber({ re: 122235.7738708, im: 91091.6398019 }),
new ComplexNumber({ re: 47625.8850387, im: 73497.3981523 }),
new ComplexNumber({ re: -15619.8231136, im: 80804.8685410 }),
new ComplexNumber({ re: 192234.0276101, im: 160833.3072355 }),
new ComplexNumber({ re: -96389.4195635, im: 393408.4543872 }),
new ComplexNumber({ re: -173449.0825417, im: 146875.7724104 }),
new ComplexNumber({ re: -179002.5662573, im: 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);
});
});

55
src/algorithms/math/fourier-transform/discreteFourierTransform.js Normal file
View File

@@ -0,0 +1,55 @@
/**
* @param {number[]} data
* @return {*[]}
*/
export default function discreteFourierTransform(data) {
const N = data.length;
const frequencies = [];
// for every frequency...
for (let frequency = 0; frequency < N; frequency += 1) {
let re = 0;
let im = 0;
// for every 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) * frequency;
// How far around the circle have we gone at time=t?
const time = t / N;
const distance = rate * time;
// Data-point * e^(-i*2*pi*f) is complex, store each part.
const rePart = data[t] * Math.cos(distance);
const imPart = data[t] * Math.sin(distance);
// add this data point's contribution
re += rePart;
im += imPart;
}
// Close to zero? You're zero.
if (Math.abs(re) < 1e-10) {
re = 0;
}
if (Math.abs(im) < 1e-10) {
im = 0;
}
// Average contribution at this frequency
re /= N;
im /= N;
frequencies[frequency] = {
re,
im,
frequency,
amp: Math.sqrt((re ** 2) + (im ** 2)),
phase: Math.atan2(im, re) * 180 / Math.PI, // in degrees
};
}
return frequencies;
}

77
src/algorithms/math/fourier-transform/fastFourierTransform.js Normal file
View File

@@ -0,0 +1,77 @@
import ComplexNumber from '../complex-number/ComplexNumber';
import bitLength from '../bits/bitLength';
/**
* Returns the number which is the flipped binary representation of input.
*
* @param {Number} [input]
* @param {Number} [bitsCount]
* @return {Number}
*/
function reverseBits(input, bitsCount) {
let reversedBits = 0;
for (let i = 0; i < bitsCount; 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 bitsCount = bitLength(inputData.length - 1);
const N = 1 << bitsCount;
while (inputData.length < N) {
inputData.push(new ComplexNumber({
real: 0,
imaginary: 0,
}));
}
const output = [];
for (let i = 0; i < N; i += 1) { output[i] = inputData[reverseBits(i, bitsCount)]; }
for (let blockLength = 2; blockLength <= N; blockLength *= 2) {
let phaseStep;
if (inverse) {
phaseStep = new ComplexNumber({
real: Math.cos(2 * Math.PI / blockLength),
imaginary: -1 * Math.sin(2 * Math.PI / blockLength),
});
} else {
phaseStep = new ComplexNumber({
real: Math.cos(2 * Math.PI / blockLength),
imaginary: Math.sin(2 * Math.PI / blockLength),
});
}
for (let blockStart = 0; blockStart < N; blockStart += blockLength) {
let phase = new ComplexNumber({
real: 1,
imaginary: 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;
}