Skip to content

Commit 12d649e

Browse files
trekhlebtrekhleb
committed
Implement Discrete Fourier Transform function.
1 parent 53a0b61 commit 12d649e

4 files changed

Lines changed: 210 additions & 38 deletions

File tree

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,7 @@ a set of rules that precisely define a sequence of operations.
7373
* `B` [Radian & Degree](src/algorithms/math/radian) - radians to degree and backwards conversion
7474
* `A` [Integer Partition](src/algorithms/math/integer-partition)
7575
* `A` [Liu Hui π Algorithm](src/algorithms/math/liu-hui) - approximate π calculations based on N-gons
76+
* `A` [Fourier Transform (DFT, FFT)](src/algorithms/math/fourier-transform) - decompose a function of time (a signal) into the frequencies that make it up
7677
* **Sets**
7778
* `B` [Cartesian Product](src/algorithms/sets/cartesian-product) - product of multiple sets
7879
* `B` [Fisher–Yates Shuffle](src/algorithms/sets/fisher-yates) - random permutation of a finite sequence

‎src/algorithms/math/fourier-transform/README.md‎

Lines changed: 41 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -78,12 +78,49 @@ Rather than jumping into the symbols, let's experience the key idea firsthand. H
7878
- *Why?* Recipes are easier to analyze, compare, and modify than the smoothie itself.
7979
- *How do we get the smoothie back?* Blend the ingredients.
8080

81+
**Think With Circles, Not Just Sinusoids**
82+
83+
The Fourier Transform is about circular paths (not 1-d sinusoids) and Euler's
84+
formula is a clever way to generate one:
85+
86+
![](https://betterexplained.com/wp-content/uploads/euler/equal_paths.png)
87+
88+
Must we use imaginary exponents to move in a circle? Nope. But it's convenient
89+
and compact. And sure, we can describe our path as coordinated motion in two
90+
dimensions (real and imaginary), but don't forget the big picture: we're just
91+
moving in a circle.
92+
93+
**Discovering The Full Transform**
94+
95+
The big insight: our signal is just a bunch of time spikes! If we merge the
96+
recipes for each time spike, we should get the recipe for the full signal.
97+
98+
The Fourier Transform builds the recipe frequency-by-frequency:
99+
100+
![](https://betterexplained.com/wp-content/uploads/images/fourier-explained-20121219-224649.png)
101+
102+
A few notes:
103+
104+
- N = number of time samples we have
105+
- n = current sample we're considering (0 ... N-1)
106+
- x<sub>n</sub> = value of the signal at time n
107+
- k = current frequency we're considering (0 Hertz up to N-1 Hertz)
108+
- X<sub>k</sub> = amount of frequency k in the signal (amplitude and phase, a complex number)
109+
- The 1/N factor is usually moved to the reverse transform (going from frequencies back to time). This is allowed, though I prefer 1/N in the forward transform since it gives the actual sizes for the time spikes. You can get wild and even use 1/sqrt(N) on both transforms (going forward and back creates the 1/N factor).
110+
- n/N is the percent of the time we've gone through. 2 * pi * k is our speed in radians / sec. e^-ix is our backwards-moving circular path. The combination is how far we've moved, for this speed and time.
111+
- The raw equations for the Fourier Transform just say "add the complex numbers". Many programming languages cannot handle complex numbers directly, so you convert everything to rectangular coordinates and add those.
112+
113+
Stuart Riffle has a great interpretation of the Fourier Transform:
114+
115+
![](https://betterexplained.com/wp-content/uploads/images/DerivedDFT.png)
116+
81117
## References
82118

83119
- [An Interactive Guide To The Fourier Transform](https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/)
84120
- [YouTube by Better Explained](https://www.youtube.com/watch?v=iN0VG9N2q0U&t=0s&index=77&list=PLLXdhg_r2hKA7DPDsunoDZ-Z769jWn4R8)
85121
- [YouTube by 3Blue1Brown](https://www.youtube.com/watch?v=spUNpyF58BY&t=0s&index=76&list=PLLXdhg_r2hKA7DPDsunoDZ-Z769jWn4R8)
86-
- [Wikipedia, FT](https://en.wikipedia.org/wiki/Fourier_transform)
87-
- [Wikipedia, DFT](https://www.wikiwand.com/en/Discrete_Fourier_transform)
88-
- [Wikipedia, DTFT](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform)
89-
- [Wikipedia, FFT](https://www.wikiwand.com/en/Fast_Fourier_transform)
122+
- Wikipedia
123+
- [FT](https://en.wikipedia.org/wiki/Fourier_transform)
124+
- [DFT](https://www.wikiwand.com/en/Discrete_Fourier_transform)
125+
- [DTFT](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform)
126+
- [FFT](https://www.wikiwand.com/en/Fast_Fourier_transform)

‎src/algorithms/math/fourier-transform/__test__/discreteFourierTransform.test.js‎

Lines changed: 137 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,143 @@
11
import discreteFourierTransform from '../discreteFourierTransform';
22

3+
/**
4+
* Helper class of the output Signal.
5+
*/
6+
class Sgnl {
7+
constructor(frequency, amplitude, phase) {
8+
this.frequency = frequency;
9+
this.amplitude = amplitude;
10+
this.phase = phase;
11+
}
12+
}
13+
314
describe('discreteFourierTransform', () => {
4-
it('should calculate split signal into frequencies', () => {
5-
const frequencies = discreteFourierTransform([1, 0, 0, 0]);
15+
it('should split signal into frequencies', () => {
16+
const testCases = [
17+
{
18+
inputAmplitudes: [1],
19+
outputSignals: [
20+
new Sgnl(0, 1, 0),
21+
],
22+
},
23+
{
24+
inputAmplitudes: [1, 0],
25+
outputSignals: [
26+
new Sgnl(0, 0.5, 0),
27+
new Sgnl(1, 0.5, 0),
28+
],
29+
},
30+
{
31+
inputAmplitudes: [2, 0],
32+
outputSignals: [
33+
new Sgnl(0, 1, 0),
34+
new Sgnl(1, 1, 0),
35+
],
36+
},
37+
{
38+
inputAmplitudes: [1, 0, 0],
39+
outputSignals: [
40+
new Sgnl(0, 0.33, 0),
41+
new Sgnl(1, 0.33, 0),
42+
new Sgnl(2, 0.33, 0),
43+
],
44+
},
45+
{
46+
inputAmplitudes: [1, 0, 0, 0],
47+
outputSignals: [
48+
new Sgnl(0, 0.25, 0),
49+
new Sgnl(1, 0.25, 0),
50+
new Sgnl(2, 0.25, 0),
51+
new Sgnl(3, 0.25, 0),
52+
],
53+
},
54+
{
55+
inputAmplitudes: [0, 1, 0, 0],
56+
outputSignals: [
57+
new Sgnl(0, 0.25, 0),
58+
new Sgnl(1, 0.25, -90),
59+
new Sgnl(2, 0.25, 180),
60+
new Sgnl(3, 0.25, 90),
61+
],
62+
},
63+
{
64+
inputAmplitudes: [0, 0, 1, 0],
65+
outputSignals: [
66+
new Sgnl(0, 0.25, 0),
67+
new Sgnl(1, 0.25, 180),
68+
new Sgnl(2, 0.25, 0),
69+
new Sgnl(3, 0.25, 180),
70+
],
71+
},
72+
{
73+
inputAmplitudes: [0, 0, 0, 2],
74+
outputSignals: [
75+
new Sgnl(0, 0.5, 0),
76+
new Sgnl(1, 0.5, 90),
77+
new Sgnl(2, 0.5, 180),
78+
new Sgnl(3, 0.5, -90),
79+
],
80+
},
81+
{
82+
inputAmplitudes: [0, 1, 0, 2],
83+
outputSignals: [
84+
new Sgnl(0, 0.75, 0),
85+
new Sgnl(1, 0.25, 90),
86+
new Sgnl(2, 0.75, 180),
87+
new Sgnl(3, 0.25, -90),
88+
],
89+
},
90+
{
91+
inputAmplitudes: [4, 1, 0, 2],
92+
outputSignals: [
93+
new Sgnl(0, 1.75, 0),
94+
new Sgnl(1, 1.03, 14),
95+
new Sgnl(2, 0.25, 0),
96+
new Sgnl(3, 1.03, -14),
97+
],
98+
},
99+
{
100+
inputAmplitudes: [4, 1, -3, 2],
101+
outputSignals: [
102+
new Sgnl(0, 1, 0),
103+
new Sgnl(1, 1.77, 8),
104+
new Sgnl(2, 0.5, 180),
105+
new Sgnl(3, 1.77, -8),
106+
],
107+
},
108+
{
109+
inputAmplitudes: [1, 2, 3, 4],
110+
outputSignals: [
111+
new Sgnl(0, 2.5, 0),
112+
new Sgnl(1, 0.71, 135),
113+
new Sgnl(2, 0.5, 180),
114+
new Sgnl(3, 0.71, -135),
115+
],
116+
},
117+
];
118+
119+
testCases.forEach((testCase) => {
120+
const { inputAmplitudes, outputSignals } = testCase;
121+
122+
// Try to split input signal into sequence of pure sinusoids.
123+
const signals = discreteFourierTransform(inputAmplitudes);
124+
125+
// Check the signal has been split into proper amount of sub-signals.
126+
expect(signals.length).toBe(outputSignals.length);
127+
128+
// Now go through all the signals and check their frequency, amplitude and phase.
129+
signals.forEach((signal, frequency) => {
130+
// Get polar form of calculated sub-signal since it is more convenient to analyze.
131+
const signalPolarForm = signal.getPolarForm(false);
132+
133+
// Get template data we want to test against.
134+
const signalTemplate = outputSignals[frequency];
6135

7-
expect(frequencies).toBeDefined();
136+
// Check all signal parameters.
137+
expect(frequency).toBe(signalTemplate.frequency);
138+
expect(signalPolarForm.radius).toBeCloseTo(signalTemplate.amplitude, 2);
139+
expect(signalPolarForm.phase).toBeCloseTo(signalTemplate.phase, 0);
140+
});
141+
});
8142
});
9143
});

‎src/algorithms/math/fourier-transform/discreteFourierTransform.js‎

Lines changed: 31 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,55 +1,55 @@
1+
import ComplexNumber from '../complex-number/ComplexNumber';
2+
13
/**
2-
* @param {number[]} data
3-
* @return {*[]}
4+
* @param {number[]} inputSignalAmplitudes - Input signal amplitudes over time (i.e. [1, 0, 4]).
5+
* @return {ComplexNumber[]} - Array of complex number. Each of the number represents the frequency
6+
* or signal. All signals together will form input signal over discrete time periods. Each signal's
7+
* complex number has radius (amplitude) and phase (angle) in polar form that describes the signal.
8+
*
9+
* @see https://gist.github.com/anonymous/129d477ddb1c8025c9ac
10+
* @see https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/
411
*/
5-
export default function discreteFourierTransform(data) {
6-
const N = data.length;
7-
const frequencies = [];
12+
export default function discreteFourierTransform(inputSignalAmplitudes) {
13+
const N = inputSignalAmplitudes.length;
14+
const outpuFrequencies = [];
815

9-
// for every frequency...
10-
for (let frequency = 0; frequency < N; frequency += 1) {
11-
let re = 0;
12-
let im = 0;
16+
// For every frequency discrete...
17+
for (let frequencyValue = 0; frequencyValue < N; frequencyValue += 1) {
18+
let signal = new ComplexNumber();
1319

14-
// for every point in time...
20+
// For every discrete point in time...
1521
for (let t = 0; t < N; t += 1) {
1622
// Spin the signal _backwards_ at each frequency (as radians/s, not Hertz)
17-
const rate = -1 * (2 * Math.PI) * frequency;
23+
const rate = -1 * (2 * Math.PI) * frequencyValue;
1824

1925
// How far around the circle have we gone at time=t?
2026
const time = t / N;
2127
const distance = rate * time;
2228

2329
// Data-point * e^(-i*2*pi*f) is complex, store each part.
24-
const rePart = data[t] * Math.cos(distance);
25-
const imPart = data[t] * Math.sin(distance);
30+
const dataPointContribution = new ComplexNumber({
31+
re: inputSignalAmplitudes[t] * Math.cos(distance),
32+
im: inputSignalAmplitudes[t] * Math.sin(distance),
33+
});
2634

27-
// add this data point's contribution
28-
re += rePart;
29-
im += imPart;
35+
// Add this data point's contribution.
36+
signal = signal.add(dataPointContribution);
3037
}
3138

3239
// Close to zero? You're zero.
33-
if (Math.abs(re) < 1e-10) {
34-
re = 0;
40+
if (Math.abs(signal.re) < 1e-10) {
41+
signal.re = 0;
3542
}
3643

37-
if (Math.abs(im) < 1e-10) {
38-
im = 0;
44+
if (Math.abs(signal.im) < 1e-10) {
45+
signal.im = 0;
3946
}
4047

4148
// Average contribution at this frequency
42-
re /= N;
43-
im /= N;
44-
45-
frequencies[frequency] = {
46-
re,
47-
im,
48-
frequency,
49-
amp: Math.sqrt((re ** 2) + (im ** 2)),
50-
phase: Math.atan2(im, re) * 180 / Math.PI, // in degrees
51-
};
49+
signal = signal.divide(N);
50+
51+
outpuFrequencies[frequencyValue] = signal;
5252
}
5353

54-
return frequencies;
54+
return outpuFrequencies;
5555
}

0 commit comments

Comments
 (0)