FFT Javascript Implementation (#599)

* implement cooley tukey in javascript

* Add Ben Chislett to CONTRIBUTORS.md

* add newline at end of file

* make approximate equality check a function in fft-js

* add words to approximation confirmation in fft-js

* attempt to fix bounds on fft-js markdown

Author: benchislett

CONTRIBUTORS.md

@@ -44,6 +44,7 @@ This file lists everyone, who contributed to this repo and wanted to show up her
 - Raven-Blue Dragon
 - Björn Heinrichs 
 - Olav Sundfør
+- Ben Chislett
 - dovisutu
 - Antetokounpo
 - Akash Dhiman

contents/cooley_tukey/code/javascript/fft.js

@@ -0,0 +1,97 @@
+const Complex = require("complex.js");
+
+function dft(x) {
+  const N = x.length;
+
+  // Initialize an array with N elements, filled with 0s
+  return Array(N)
+    .fill(new Complex(0, 0))
+    .map((temp, i) => {
+      // Reduce x into the sum of x_k * exp(-2*sqrt(-1)*pi*i*k/N)
+      return x.reduce((a, b, k) => {
+        return a.add(b.mul(new Complex(0, (-2 * Math.PI * i * k) / N).exp()));
+      }, new Complex(0, 0)); // Start accumulating from 0
+    });
+}
+
+function cooley_tukey(x) {
+  const N = x.length;
+  const half = Math.floor(N / 2);
+  if (N <= 1) {
+    return x;
+  }
+
+  // Extract even and odd indexed elements with remainder mod 2
+  const evens = cooley_tukey(x.filter((_, idx) => !(idx % 2)));
+  const odds = cooley_tukey(x.filter((_, idx) => idx % 2));
+
+  // Fill an array with null values
+  let temp = Array(N).fill(null);
+
+  for (let i = 0; i < half; i++) {
+    const arg = odds[i].mul(new Complex(0, (-2 * Math.PI * i) / N).exp());
+
+    temp[i] = evens[i].add(arg);
+    temp[i + half] = evens[i].sub(arg);
+  }
+
+  return temp;
+}
+
+function bit_reverse_idxs(n) {
+  if (!n) {
+    return [0];
+  } else {
+    const twice = bit_reverse_idxs(n - 1).map(x => 2 * x);
+    return twice.concat(twice.map(x => x + 1));
+  }
+}
+
+function bit_reverse(x) {
+  const N = x.length;
+  const indexes = bit_reverse_idxs(Math.log2(N));
+  return x.map((_, i) => x[indexes[i]]);
+}
+
+// Assumes log_2(N) is an integer
+function iterative_cooley_tukey(x) {
+  const N = x.length;
+
+  x = bit_reverse(x);
+
+  for (let i = 1; i <= Math.log2(N); i++) {
+    const stride = 2 ** i;
+    const half = stride / 2;
+    const w = new Complex(0, (-2 * Math.PI) / stride).exp();
+    for (let j = 0; j < N; j += stride) {
+      let v = new Complex(1, 0);
+      for (let k = 0; k < half; k++) {
+        // perform butterfly multiplication
+        x[k + j + half] = x[k + j].sub(v.mul(x[k + j + half]));
+        x[k + j] = x[k + j].sub(x[k + j + half].sub(x[k + j]));
+        // accumulate v as powers of w
+        v = v.mul(w);
+      }
+    }
+  }
+
+  return x;
+}
+
+// Check if two arrays of complex numbers are approximately equal
+function approx(x, y, tol = 1e-12) {
+  let diff = 0;
+  for (let i = 0; i < x.length; i++) {
+    diff += x[i].sub(y[i]).abs();
+  }
+  return diff < tol;
+}
+
+const X = Array.from(Array(8), () => new Complex(Math.random(), 0));
+const Y = cooley_tukey(X);
+const Z = iterative_cooley_tukey(X);
+const T = dft(X);
+
+// Check if the calculations are correct within a small tolerance
+console.log("Cooley tukey approximation is accurate: ", approx(Y, T));
+console.log("Iterative cooley tukey approximation is accurate: ", approx(Z, T));

contents/cooley_tukey/cooley_tukey.md

@@ -85,6 +85,8 @@ For some reason, though, putting code to this transformation really helped me fi
 [import:4-13, lang:"julia"](code/julia/fft.jl)
 {% sample lang="asm-x64" %}
 [import:15-74, lang:"asm-x64"](code/asm-x64/fft.s)
+{% sample lang="js" %}
+[import:3-15, lang:"javascript"](code/javascript/fft.js)
 {% endmethod %}
 
 In this function, we define `n` to be a set of integers from $$0 \rightarrow N-1$$ and arrange them to be a column.
@@ -101,7 +103,7 @@ M = [1.0+0.0im  1.0+0.0im           1.0+0.0im          1.0+0.0im;
 
 It was amazing to me when I saw the transform for what it truly was: an actual transformation matrix!
 That said, the Discrete Fourier Transform is slow -- primarily because matrix multiplication is slow, and as mentioned before, slow code is not particularly useful.
-So what was the trick that everyone used to go from a Discrete Fourier Transform to a *Fast* Fourier Transform?
+So what was the trick that everyone used to go from a Discrete Fourier Transform to a _Fast_ Fourier Transform?
 
 Recursion!
 
@@ -134,6 +136,8 @@ In the end, the code looks like:
 [import:16-32, lang:"julia"](code/julia/fft.jl)
 {% sample lang="asm-x64" %}
 [import:76-165, lang:"asm-x64"](code/asm-x64/fft.s)
+{% sample lang="js" %}
+[import:17-39, lang="javascript"](code/javascript/fft.js)
 {% endmethod %}
 
 As a side note, we are enforcing that the array must be a power of 2 for the operation to work.
@@ -142,6 +146,7 @@ This is a limitation of the fact that we are using recursion and dividing the ar
 The above method is a perfectly valid FFT; however, it is missing the pictorial heart and soul of the Cooley-Tukey algorithm: Butterfly Diagrams.
 
 ### Butterfly Diagrams
+
 Butterfly Diagrams show where each element in the array goes before, during, and after the FFT.
 As mentioned, the FFT must perform a DFT.
 This means that even though we need to be careful about how we add elements together, we are still ultimately performing the following operation:
@@ -214,8 +219,8 @@ I'll definitely come back to this at some point, so let me know what you liked a
 
 To be clear, the example code this time will be complicated and requires the following functions:
 
-* An FFT library (either in-built or something like FFTW)
-* An approximation function to tell if two arrays are similar
+- An FFT library (either in-built or something like FFTW)
+- An approximation function to tell if two arrays are similar
 
 As mentioned in the text, the Cooley-Tukey algorithm may be implemented either recursively or non-recursively, with the recursive method being much easier to implement.
 I would ask that you implement either the recursive or non-recursive methods (or both, if you feel so inclined).
@@ -242,9 +247,10 @@ Note: I implemented this in Julia because the code seems more straightforward in
 Some rather impressive scratch code was submitted by Jie and can be found here: https://scratch.mit.edu/projects/37759604/#editor
 {% sample lang="asm-x64" %}
 [import, lang:"asm-x64"](code/asm-x64/fft.s)
+{% sample lang="js" %}
+[import, lang:"javascript"](code/javascript/fft.js)
 {% endmethod %}
 
-
 <script>
 MathJax.Hub.Queue(["Typeset",MathJax.Hub]);
 </script>
@@ -262,6 +268,7 @@ The text of this chapter was written by [James Schloss](https://github.com/leios
 [<p><img  class="center" src="../cc/CC-BY-SA_icon.svg" /></p>](https://creativecommons.org/licenses/by-sa/4.0/)
 
 ##### Images/Graphics
+
 - The image "[FTexample](res/FT_example.png)" was created by [James Schloss](https://github.com/leios) and is licenced under the [Creative Commons Attribution-ShareAlike 4.0 International License](https://creativecommons.org/licenses/by-sa/4.0/legalcode).
 - The image "[radix2positive](res/radix-2screen_positive.jpg)" was created by [James Schloss](https://github.com/leios) and is licenced under the [Creative Commons Attribution-ShareAlike 4.0 International License](https://creativecommons.org/licenses/by-sa/4.0/legalcode).
 - The image "[radix2](res/radix-2screen.jpg)" was created by [James Schloss](https://github.com/leios) and is licenced under the [Creative Commons Attribution-ShareAlike 4.0 International License](https://creativecommons.org/licenses/by-sa/4.0/legalcode).
@@ -271,4 +278,5 @@ The text of this chapter was written by [James Schloss](https://github.com/leios
 ##### Pull Requests
 
 After initial licensing ([#560](https://github.com/algorithm-archivists/algorithm-archive/pull/560)), the following pull requests have modified the text or graphics of this chapter:
+
 - none