9 ms·

Show HN: A minimal FFT code

100 points by martyalain 8y ago

InGodsName 8y ago
Would be cool if it shows code in most common programming languages
- Avshalom 8y ago
  Most implementations on http://rosettacode.org/wiki/Fast_Fourier_transform http://rosettacode.org/wiki/Fast_Fourier_transform use the same algorithm.
  - martyalain 8y ago
    I didn't see one alike, in a pure functional way. (http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambdatalk http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambd...)
- gp7 8y ago
  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-Practice https://cnx.org/contents/ulXtQbN7@15/Implementing-FFTs-in-Pr...
  - martyalain 8y ago
    Actually my quest is to write the FFT code at the lambda-calculus level. Why? First for the fun and also because I consider that lambda-calculus is for the mind what assembly language is for the computer. See http://lambdaway.free.fr/lambdaspeech/?view=PLR http://lambdaway.free.fr/lambdaspeech/?view=PLR. In this page I would like to replace the inefficient unary numeration by a more efficient decimal position numeration and I guess that FFT, the Katsuba or any divide & conquer algorithm could be useful. It should be pleasant for the mind and could overcome limits of the JS numbers implementation.
  - martyalain 8y ago
    "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 8y ago
      You copy rather than use strided accesses. It's the same algorithm.
      - martyalain 8y ago
        " 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 8y ago
        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 8y ago
  You are right. You can find a Javascript version in this page http://lambdaway.free.fr/lambdaspeech/?view=lispology http://lambdaway.free.fr/lambdaspeech/?view=lispology and the initial LISP code in https://www.physik.uzh.ch/~psaha/mus/fourlisp.php https://www.physik.uzh.ch/~psaha/mus/fourlisp.php on which the {lambda talk} code has been built with some missing examples.
- basementcat 8y ago
  Here are a few implementations in various languages. https://codegolf.stackexchange.com/questions/12420/too-fast-too-fourier-fft-code-golf https://codegolf.stackexchange.com/questions/12420/too-fast-...
dpwm 8y ago
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
- martyalain 8y ago
  And {lambda talk} works in a small wiki (25kb zipped folder) working as a light overlay on any modern web browsers, taking benefit of their powerfull capabilities (HTML/CSS/DOM/Javascript)
- nahtnam 8y ago
  The syntax kind of looks like Scheme to me.
  - martyalain 8y ago
    Yes, it's a translation of https://www.physik.uzh.ch/~psaha/mus/fourlisp.php https://www.physik.uzh.ch/~psaha/mus/fourlisp.php
gregfjohnson 8y ago
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 8y ago
  Thank you :)
jacobolus 8y ago
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 8y ago
  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
  - deleted 8y ago
    [deleted]
salvad0r 8y ago
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 8y ago
  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
  - martyalain 8y ago
    This one is fine too: https://bl.ocks.org/jinroh/7524988 https://bl.ocks.org/jinroh/7524988
  - jmiserez 8y ago
    I like this one by 3Blue1Brown: https://youtu.be/spUNpyF58BY https://youtu.be/spUNpyF58BY
- kazinator 8y ago
  Do you understand the DFT? The FFT just calculates that faster. There is an intuition to it, something like this: suppose we want to calculate a 2048 element DFT. If we instead calculate a pair of 1024 element DFT's over halves of data, we have all the same high frequency information there in the two windows. What we're missing is the lower frequency: the one that goes through just one cycle over the 2048 window. But we don't need 2048 points for that; the low frequency doesn't contain that much information. The FFT reveals how it can be obtained from the two halves. The two halves have the necessary information in their lowest frequency component; we don't need to sample 2048 points of the signal; we just look into the half-sized FFT results and put that together.
  - betterunix2 8y ago

100 points by martyalain 8y ago

InGodsName 8y ago

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

Avshalom 8y ago
Most implementations on http://rosettacode.org/wiki/Fast_Fourier_transform http://rosettacode.org/wiki/Fast_Fourier_transform use the same algorithm.
- martyalain 8y ago
  I didn't see one alike, in a pure functional way. (http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambdatalk http://www.rosettacode.org/wiki/Fast_Fourier_transform#lambd...)
gp7 8y ago
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-Practice https://cnx.org/contents/ulXtQbN7@15/Implementing-FFTs-in-Pr...
- martyalain 8y ago
  Actually my quest is to write the FFT code at the lambda-calculus level. Why? First for the fun and also because I consider that lambda-calculus is for the mind what assembly language is for the computer. See http://lambdaway.free.fr/lambdaspeech/?view=PLR http://lambdaway.free.fr/lambdaspeech/?view=PLR. In this page I would like to replace the inefficient unary numeration by a more efficient decimal position numeration and I guess that FFT, the Katsuba or any divide & conquer algorithm could be useful. It should be pleasant for the mind and could overcome limits of the JS numbers implementation.
- martyalain 8y ago
  "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 8y ago
    You copy rather than use strided accesses. It's the same algorithm.
    - martyalain 8y ago
      " 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 8y ago
        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 8y ago
You are right. You can find a Javascript version in this page http://lambdaway.free.fr/lambdaspeech/?view=lispology http://lambdaway.free.fr/lambdaspeech/?view=lispology and the initial LISP code in https://www.physik.uzh.ch/~psaha/mus/fourlisp.php https://www.physik.uzh.ch/~psaha/mus/fourlisp.php on which the {lambda talk} code has been built with some missing examples.
basementcat 8y ago
Here are a few implementations in various languages. https://codegolf.stackexchange.com/questions/12420/too-fast-too-fourier-fft-code-golf https://codegolf.stackexchange.com/questions/12420/too-fast-...

dpwm 8y ago

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

martyalain 8y ago
And {lambda talk} works in a small wiki (25kb zipped folder) working as a light overlay on any modern web browsers, taking benefit of their powerfull capabilities (HTML/CSS/DOM/Javascript)
nahtnam 8y ago
The syntax kind of looks like Scheme to me.
- martyalain 8y ago
  Yes, it's a translation of https://www.physik.uzh.ch/~psaha/mus/fourlisp.php https://www.physik.uzh.ch/~psaha/mus/fourlisp.php

gregfjohnson 8y ago

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 8y ago
Thank you :)

jacobolus 8y ago

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 8y ago
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
- deleted 8y ago
  [deleted]

salvad0r 8y ago

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 8y ago
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
- martyalain 8y ago
  This one is fine too: https://bl.ocks.org/jinroh/7524988 https://bl.ocks.org/jinroh/7524988
- jmiserez 8y ago
  I like this one by 3Blue1Brown: https://youtu.be/spUNpyF58BY https://youtu.be/spUNpyF58BY
kazinator 8y ago
Do you understand the DFT? The FFT just calculates that faster. There is an intuition to it, something like this: suppose we want to calculate a 2048 element DFT. If we instead calculate a pair of 1024 element DFT's over halves of data, we have all the same high frequency information there in the two windows. What we're missing is the lower frequency: the one that goes through just one cycle over the 2048 window. But we don't need 2048 points for that; the low frequency doesn't contain that much information. The FFT reveals how it can be obtained from the two halves. The two halves have the necessary information in their lowest frequency component; we don't need to sample 2048 points of the signal; we just look into the half-sized FFT results and put that together.
- betterunix2 8y ago