
Show HN: A minimal FFT code - martyalain
http://lambdaway.free.fr/lambdaspeech/?view=zorg
======
gregfjohnson
I wrote an implementation of FFT based on a few choices that people may find
useful. It is at:

[https://github.com/gregfjohnson/fft](https://github.com/gregfjohnson/fft)

It is not the fastest or most full-featured, but it does have a few useful
features. 1 - written in C, small, self-contained, trivial to download and
build. 2 - uses mixed radix implementation, so data files with length other
than powers of two work ok, depending on the largest prime factor of the
length of the input. 3 - input and output formats either real scalar,
rectangular, or polar complex values. 4 - convenient to use with other shell
commands; works with pipes or files.

I am super-thankful for open source and all of the dedicated contributors out
there, and want to give back however I can!

~~~
martyalain
Thank you :)

------
jacobolus
A lisp implementation is cute and can be a nice teaching aid but is slower
than desirable for practical applications.

What I would really like to see is a reasonably short (small code size) but
reasonably optimized (for an arbitrary machine that might run a browser) Web
Assembly implementation of a power-of-2 FFT, possibly compiled down from
another language, e.g. C or Rust. Throw in an implementation of Bluestein’s
non-power-of-2 FFT for bonus points.

The existing Javascript implementations I have seen are much slower than they
could be if a serious expert took a crack at it, while the existing compiled-
to-Javascript or -wasm implementations weren’t originally designed with the
web as a target, and tend to be bloated monsters.

~~~
sheltron
have you seen this? I am not an expert in ffts so maybe it could be improved
but it is impressively fast in my experience (also very readable) :
[https://github.com/corbanbrook/dsp.js/blob/master/dsp.js](https://github.com/corbanbrook/dsp.js/blob/master/dsp.js)

------
dpwm
I'm familiar with most functional languages. This post was in a language with
slightly unusual syntax. Still, I found what this is written in to be
surprisingly quite readable.

It turns out it is {lambda talk}, a simple functional language for the site.
[0]

[0]
[http://lambdaway.free.fr/lambdaspeech/?view=start](http://lambdaway.free.fr/lambdaspeech/?view=start)

~~~
nahtnam
The syntax kind of looks like Scheme to me.

~~~
martyalain
Yes, it's a translation of
[https://www.physik.uzh.ch/~psaha/mus/fourlisp.php](https://www.physik.uzh.ch/~psaha/mus/fourlisp.php)

------
salvad0r
FFT's were very difficult for me to understand (and still are). The graphic at
the top is actually a very good way to show what a practical FFT does. I wish
I had it back when I first started trying to understand them.

~~~
giornogiovanna
That's a good visualization. I also like this one[0], where the spinning rate
of each radius is fixed, and the FFT's job is to find the size of each radius.
Another great one is imagining the DFT matrix (which is unitary) as a change
of basis.

[0]:
[https://www.youtube.com/watch?v=k8FXF1KjzY0](https://www.youtube.com/watch?v=k8FXF1KjzY0)

~~~
jmiserez
I like this one by 3Blue1Brown:
[https://youtu.be/spUNpyF58BY](https://youtu.be/spUNpyF58BY)

------
InGodsName
Would be cool if it shows code in most common programming languages

~~~
gp7
The code is mostly equivalent to the recursive form found in the wikipedia
article on the cooley-tukey algorithm. This is a good one to learn from, as
its not only a simple formulation but it forms the basis of modern optimised
FFTs such as FFTW (source[1], from the authors of FFTW).

As an aside, I also find the non-recursive, breadth-first, form easy to derive
thru a process of code transformations of the depth-first form; explanations
that start breadth-first are somewhat bewildering

[1] [https://cnx.org/contents/ulXtQbN7@15/Implementing-FFTs-in-
Pr...](https://cnx.org/contents/ulXtQbN7@15/Implementing-FFTs-in-Practice)

~~~
martyalain
"The code is mostly equivalent to the recursive form found in the wikipedia
article on the cooley-tukey algorithm" I don't think so. A JS translation of
the code shown in wikipedia can be seen in
[http://lambdaway.free.fr/lambdaspeech/?view=lispology](http://lambdaway.free.fr/lambdaspeech/?view=lispology).

~~~
gp7
You copy rather than use strided accesses. It's the same algorithm.

~~~
martyalain
" You copy rather than use strided accesses. " I don't understand "strided
accesses". « Plagiarism is stealing, copying is create. » I just translated
the code in the lambdatalk language and added missing examples, something that
I consider mandatory.

~~~
gp7
Sorry, I'm not accusing you of anything. By strided accesses I mean when you
access elements of an array sequentially using some step size. Think loops
with `i += stride` instead of `i++`. In the wikipedia pseudocode this is done
implicitly using the 's' parameter. Notice that, in the version you
implemented, you split the input into even and odd parts explicitly; you can
achieve the same end by accessing the input array in a certain order as you
are performing the mathematical operations. This is what the wikipedia
pseudocode does. If you've seen other versions of the FFT with a bit reversal
step, this is also where that comes in.

Check this out (javascript):

    
    
      function permute1(x) { 
          if (x.length == 1) return x;
          let even = [];
          let odd = []; 
          for (let i = 0; i < x.length; i += 2) { 
               even[i / 2] = x[i];
               odd[i / 2] = x[i + 1];
          }
          return [].concat(permute1(even), permute1(odd));
      }
    
      function permute2(x, offset, stride) {
          if (!offset) offset = 0;
          if (!stride) stride = 1;
          if (stride >= x.length) return [x[offset]];
          return [].concat(permute2(x, offset, stride * 2), permute2(x, offset + stride, stride * 2));
      }
      
      function permute3(x) {
          let result = [];
          for (let i = 0; i < x.length; i++) {
              let k = i;
              // pretend 32-bit ints
              k = ((k >> 1) & 0x55555555) | ((k & 0x55555555) << 1);
              k = ((k >> 2) & 0x33333333) | ((k & 0x33333333) << 2);
              k = ((k >> 4) & 0x0F0F0F0F) | ((k & 0x0F0F0F0F) << 4);
              k = ((k >> 8) & 0x00FF00FF) | ((k & 0x00FF00FF) << 8);
              k = ( k >> 16             ) | ( k               << 16); 
              k = k >> (64 - Math.log2(x.length));
              if (k < 0) k += x.length; // fix up due to signed ints
              result[i] = x[k];
          }
          return result;
      }
    

For arrays with power of two sizes, these perform the same permutation (but
fail differently for non power of two sizes). Note that, with permute1, we
effectively iterate over the entire input log2(n) times, so this is an
O(nlogn) algorithm!

edit: also, i think i may have misunderstood the relationship between your js
version and your lambdatalk version. They seem to be the same to me?

~~~
martyalain
Thank you for your clever code. This morning I improved a little the JS code
in
[http://lambdaway.free.fr/lambdaspeech/index.php?view=lispolo...](http://lambdaway.free.fr/lambdaspeech/index.php?view=lispology)
and I plan to learn, understand and insert your code in the fft() function.

Yes there is a relationship between the JS version and its translation into
lambdatalk. My project is to replace the array based version by a list based
version so that I can replace in this page
[http://lambdaway.free.fr/lambdaspeech/?view=PLR](http://lambdaway.free.fr/lambdaspeech/?view=PLR)
the inefficient unary numeration based implementation of numbers (using
standard Church numbers or just lists) by a decimal position numeration.
Standard multiplication of words seen as polynoms being O(n^2) I need to go
further and implement fast multiplication. So my interest in FFT.

As you could see in
[http://lambdaway.free.fr/lambdaspeech/meca/JS.js](http://lambdaway.free.fr/lambdaspeech/meca/JS.js),
the lambdatalk's interpreter is a regular expression window running on the
code (not an AST) and replacing in situ expressions by their values. A kind of
Turing machine. I like the idea of overcoming limits of JS numbers using
nothing but words and simple substitutions on words.

