This is part 3 of a series of posts about FFTs and the specific algorithms that
FFT-20XX uses.
- The Discrete Fourier Transform
- The Basics of the FFT
- Optimizations to the FFT (← you are here)
- The Inverse FFT
- The FFT-20XX Complex-Input Algorithm
- The FFT-20XX Real-Input Algorithm(s)
In the last post we got to a basic implementation of the radix 2 decimation-in-time Cooley-Tukey algorithm. But there are still other additions and optimizations we can make before
getting into the way FFT-20XX works. We’ll start by diving into the “radix 2” part of the algorithm: we can make
improvements there.
This is part 2 of a series of posts about FFTs and the specific algorithms that
FFT-20XX uses.
- The Discrete Fourier Transform
- The Basics of the FFT (← you are here)
- Optimizations to the FFT
- The Inverse FFT
- The FFT-20XX Complex-Input Algorithm
- The FFT-20XX Real-Input Algorithm(s)
Part 1 went through the core of how the
discrete Fourier transform (DFT) works to turn a series of
samples in time into a set of samples in frequency, and ended up with the following algoritm:
void CalculateDFT(
complex[] inputs,
complex[] outputs,
int length)
{
for (outIdx = 0; outIdx < length; outIdx++)
{
float sum = 0
float angularFreq = -2 * pi * outIdx / length
for (inIdx = 0; inIdx < length; inIdx++)
{
float radians = inIdx * angularFreq
complex rotation = cos(radians) + i*sin(radians)
sum += inputs[inIdx] * rotation
}
outputs[outIdx] = sum
}
}
A reminder: we’re going to be talking about these values as complex numbers (of the form ), but if you prefer,
just think of them as a 2D vector instead ( where is the x coordinate and is the y coordinate).
We’ll be doing complex multiplication of these vectors with others, but in all cases a multiply like will
be a rotation of point by the angle that represents (where ).
This algorithm is , but there’s a much more efficient, way: the
fast Fourier transform (FFT) – and we’re going to derive it.
Because the world clearly doesn’t already have enough of them, I recently released a new FFT library:
FFT-20XX. To the best of my ability to tell, it uses an algorithm that is
at least slightly novel, and I thought it might be interesting to write about it.
This is the first of a series of posts about FFTs and the specific algorithms that FFT-20XX uses:
- The Discrete Fourier Transform (← you are here)
- The Basics of the FFT
- Optimizations to the FFT
- The Inverse FFT
- The FFT-20XX Complex-Input Algorithm
- The FFT-20XX Real-Input Algorithm(s)
As a first step I thought it might be fun to start with “what even is the
discrete Fourier transform?”
At some point, if you’ve done multithreaded programming, you’ve probably used a mutex (or some other locking mechanism). Locks are relatively straightforward to understand and use (“lock this thing before accessing your data and unlock it when you’re done”), but they do have their issues. The most commonly-discussed issue with using locks is the dreaded deadlock, scourge of many a poor soul who needed to hold multiple locks at once for something.
But we’re not here to talk about deadlocks…instead, I’d like to focus on a problem that I have encountered much more frequently: accidentally reading or modifying protected data without having acquired the lock.
That’s right, we’re once again going to try to protect ourselves from our worst enemy: ourselves.
As with the previous entry in this semi-series, the examples and implementation are in C++, but you can probably build something similar in your language of choice (unless it doesn’t need it).
An Easy Mistake
It turns out it’s surprisingly easy to not grab a lock before reading or writing things that need to be synchronized. Here’s a toy example:
void Foo::IncrementThing(int c)
{
std::lock_guard lock { m_mutex };
m_thing += c;
}
void Foo::DecrementThing(int c)
{
m_thing -= c;
}
In the real world, “oops I accessed a thing outside of the lock” bugs tend to be more subtle - it can be especially hard to notice when a thing that should be present isn’t. Since all of the accesses are to normal members of your class/struct/global scope/whatever, it’s easy to slip and put an access to a value where it wasn’t intended, and wasn’t propertly protected.
But what if instead you could do something more like the following:
class Foo
{
MutexProtected m_state
{
int m_thing = 0;
};
};
void Foo::IncrementThing(int c)
{
lock (m_state)
{
m_thing += c;
}
}
Basically, what if you could have some sort of protected wrapper around the state that walls it off and makes it inaccessible except when the lock is actually held? With something like that it would become much more difficult to get at values when it shouldn’t be allowed.
Additionally, it would help both with organization (forcing the mutex-protected values together in the code) as well as making the intent clear: values that are protected sit inside the protected block, clarifying which values are (and are not) intended to be accessed solely through the mutex, even to folks who were not the original author (the list of which stealthily includes the original author, one month in the future).
Imagine this: you’ve got some value, x, that you want to ensure is at least 1. That is to say, you want to ensure its minimum value is 1. So, being the smart, experienced programmer that you are, you write the following:
x = Min(x, 1);
You give yourself the small, satisfied nod of a job well done and run the program and then it all goes immediately sideways because that should have been Max and not Min.
If you’ve been writing code for basically any length of time, the above was probably less an “imagine this” and more a “remember this” because if you’re anything like me, you’ve done this over. and over. and over.
Inspired by once again mistakenly using Min instead of Max to limit the minimum allowed value of something, I’ve decided to start a little series (will it have more than one entry? who knows!) called Protecting Coders From Ourselves, in which we rework some bit of API surface to make it less error-prone. We’re going to deal with the “I chose the wrong Min/Max again” problem, but first I want to talk about Clamp and Lerp.
A while back I was working on a programming language idea and while I haven’t made any progress on it in ages, I really liked the string design that I came up with. I don’t know that any of the ideas are original, but I haven’t seen anything exactly like it so I figure I’d throw the idea out into the ether in case anyone else happens to do something similar in their own personal language that they definitely shouldn’t be making 😆
(I’ll note upfront: this design uses dollar signs ($) and backticks(`). There are many languages that do so (like Javascript!) but these keys are not universally on all keyboards internationally so for those locales it may be more difficult to type these out…my language design was pretty much just for me with my standard US keyboard so I didn’t take this into account)
(This post follows Part 1: 32-bit floats and will make very little sense without having read that one first. Honestly, it might make little sense having read that one first, I dunno!)
Last time we went over how to calculate the results of the FMAdd instruction (a fused-multiply-add calculated as if it had infinite internal precision) for 32-bit single-precision float values:
- Calculate the double-precision product of
a and b
- Add this product to
c to get a double-precision sum
- Calculate the error of the sum
- Use the error to odd-round the sum
- Round the double-precision sum back down to single precision
This requires casting up to a 64-bit double-precision float to get extra bits of precision. But what if you can’t do that? What if you’re using doubles? You can’t just (in most cases) cast up to a quad-precision float. So what do you do?
Double-Precision FMAdd
(Like the 32-bit version, this is based on Emulation of FMA and correctly-rounded sums: proved algorithms using rounding to odd by Sylvie Boldo and Guillaume Melquiond, there are additional details in there)
To do this natively as doubles, we need to invent a new operation: MulWithError. This is the multiplication equivalent of the AddWithError function from the 32-bit solution:
(double prod, double err) MulWithError(double x, double y)
{
double prod = x * y;
double err =
return (prod, err);
}
We’ll get to how to implement that in a moment, but first we’ll walk through how to use that function to calculate a proper FMAdd.
A thing that I had to do at work is write an emulation of the FMAdd (fused multiply-add) instruction for hardware where it wasn’t natively supported (specifically I was writing a SIMD implementation, but the idea is the same), and so I thought I’d share a little bit about how FMAdd works, since I’ve already been posting about how float rounding works.
So, screw it, here we go with another unnecessarily technical, mathy post!
What is the FMAdd Instruction?
A fused multiply-add is basically doing a multiply and an add as a single operation, and it gives you the result as if it were computed with infinite precision and then rounded down at the final result. FMAdd computes (a * b) + c without intermediate floating-point error being introduced:
float FMAdd(float a, float b, float c)
{
return (a * b) + c;
}
Computing it normally (using the code above) for some values will get you double rounding (explained in a moment) which means you might be an extra bit off (or, more formally, one ULP) from where your actual result should be. An extra bit doesn’t sound like a lot, but it can add up over many operations.
Fused multiply-add avoids this extra rounding, making it more accurate than a multiply followed by a separate add, which is great! (It can also be faster if it’s supported by hardware but, as you’ll see, computing it without a dedicated instruction on the CPU is actually surprisingly spendy, especially once you get into doing it for 64-bit floats, but sometimes you need precision instead of performance).
C++17 added support for hex float literals, so you can put more bit-accurate floating point values into your code. They’re handy to have, and I wanted to be able to parse them from a text file in a C# app I was writing.
Some examples:
0x1.0p0
0x1.8p1
0x8.0p-3
0x0.8p1
0xAB.CDEFp-10
0x0.0000000ABp0
What Am I Even Looking At Here?
(Before reading this, if you don’t have a good feel for how floats work, maybe check my previous post about floats (and float rounding))
I had a bit of a mental block on this number format for a bit - like, what does it even mean to have fractional hex digits? But it turns out it’s a concept that we already use all the time and my brain just needed some prodding to make the connection.
With our standard base 10 numbers, moving the decimal point left one digit means dividing the number by 10:
12.3 == 1.23 * 10^1 == 0.123 * 10^2
Hex floats? Same deal, just in 16s instead:
0x1B.C8 == 0x1.BC8 * 16^1 == 0x0.1BC8 * 16^2
Okay, so now what’s the “p” part in the number? Well, that’s the start of the exponent. A standard float has an exponent starting with ‘e’:
1.3e2 == 1.3 * 10^2
But ‘e’ is a hex digit, so you can’t use ‘e’ anymore as the exponent starter, so they chose ‘p’ instead (why not ‘x’, the second letter? Probably because a hex number starts with ‘0x’, so ‘x’ also already has a use - but ‘p’ is free so it wins)
The exponent for a hex float is in powers of 2 (so it corresponds perfectly to the exponent as it is stored in the value), so:
0x1.ABp3 == 0x1.AB * 2^3
So that’s how a hex float literal works! Here’s a quick breakdown:
<hex-digits>[ '.' [fractional-hex-digits]] 'p' <exponent>
(Parsing code below the fold)
I’ve ported the blog off of Wordpress onto a static site generator (Specifically, 11ty/Eleventy). I have a bit more control over the format here, and it’s easier for me to write pages (Wordpress was fighting me on all sorts of formatting which I can just do now).
What this means is I can finally start copying over the rest of my posts that were on (the now-defunct, sadly) cohost.org (rest easy, little eggbug).
Likely there are things on the new site that aren’t set up correctly yet, so if you happen to notice anything, find me on Bluesky or Mastodon and let me know!
And now, for fun, here’s a pic I took in November at Arches National Park:
