Show HN: A minimal FFT code
lambdaway.free.fr
lambdaway.free.fr
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!
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.
It turns out it is {lambda talk}, a simple functional language for the site. [0]
That machine is incredible.
I can't begin to imagine the amount of work it took to design and buuld sonething like that.
Also, the video is very well produced.
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.
There are two reasons I like this approach. The first is that it is the only way to understand the number-theoretic transform, which is very useful in cryptography (among other things). The second is that you can use the same reasoning to divide the DFT in other ways (e.g. into thirds, or fifths, etc.). The downside is that it is more abstract so if all you care about is signal processing you need to somehow connect abstract "roots of unity" with complex frequencies (but that is probably not too bad if you are comfortable with complex frequencies).
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...
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?
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 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, 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.