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!
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.
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]
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)
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.
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.
I actually find that terms like "frequency" are a distraction here. Maybe I am odd, but for me the FFT is easier to understand when you think of the DFT in terms of roots of unity (i.e. nth roots of 1), and writing out the DFT as a matrix involving those roots. Basically, if you take only the first half of the even rows of the DFT matrix of size NxN, you get a DFT matrix of size (N/2)x(N/2); this is basically because if W is an Nth root of 1 then W^2 is an (N/2)th root. The other key observation is that the elements of the odd rows can be computed from the first half of the even rows by multiplying by powers of W (this is basic algebra), so from an (N/2)x(N/2) DFT matrix you can compute an NxN DFT matrix (hence the FFT).
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).
Roots of unity are easily related to frequencies. They are just ways of walking around the unit circle in equal steps. If we want N steps, we take the N-th power roots. These walks around the unit circle give us a discretized sine and its quadrature cosine. It's all together because rotation is multiplication, and harmonic motion is related to rotation.
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.
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
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. 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.
"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.
" 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.
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?
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.
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!