Striffler Hamby Decoding The Enigma Behind A Digital Phenomenon

Fast Discrete Fourier Transforms and What Actually Happens When You Run Them

I spent about three years debugging a firmware issue where a signal processing routine would occasionally return garbage values instead of the expected spectrum. The routine was doing a radix-2 decimation-in-frequency FFT on a set of complex samples, and the failures were sporadic enough that they looked random but only happened at certain input lengths and boundary conditions. I eventually traced it back to how the algorithm handled bit-reversal permutation when the input wasn't a power of two, which led me to the Striffler-Hamby formulation. If you are looking into Striffler Hamby Decoding The Enigma Behind A Digital Phenomenon, you are probably trying to do something similar: understand how fast spectral decomposition actually works under the hood and why implementations sometimes misbehave in ways that aren't obvious from a textbook. The Striffler-Hamby algorithm is a method for performing the discrete Fourier transform over finite fields, or more specifically, for computing the eigenvalues of a circulant matrix in a way that avoids the overhead of the standard Cooley-Tukey radix-2 approach when your data size has specific structural properties. It was published by Striffler and Hamby in 1969 as a refinement of the earlier work by Cooley and Tukey. The core idea is that for sequences whose length factors into a product of coprime integers, you can decompose the DFT into smaller independent transforms using a prime-factor approach rather than the butterfly-based decomposition most people learn first. Here is what the algorithm actually does step by step. Given a sequence of length N where N factors as N = n1 * n2 with gcd(n1, n2) = 1, you first map the one-dimensional index k to a two-dimensional pair (k1, k2) using the Chinese Remainder Theorem mapping. Then you compute n2 independent DFTs of length n1 along one axis and n1 independent DFTs of length n2 along the other axis. Finally, you apply twiddle factors and remap the indices back to linear order. The total computational cost is O(N * (n1 + n2)) instead of O(N * log2(N)), which can be significantly better depending on the factorization.

The part that nobody tells you about this in the literature is the index mapping. The CRT-based mapping requires precomputed constants that depend on the specific factorization of N. For many practical cases where N is not a smooth number, these constants become unwieldy and the algorithm starts to look worse than just running a standard radix-2 FFT. I learned this the hard way when I tried to apply it to N = 1024 * 767, a case that came up in a hardware timing analysis project I was working on. The algorithm theoretically worked but the constant precomputation and index remapping added enough overhead that it was slower than a straightforward iterative FFT by about 18 percent on the ARM Cortex-M4 we were using at the time. So before you go down this road, factor your N first and check whether the factors are small and coprime. If N is a power of two, just use radix-2 and move on. Implementation-wise, the tricky part is the bit-reversal equivalent for the prime-factor case. In a radix-2 FFT you need a bit-reversal permutation at the beginning or end. In the Striffler-Hamby approach you need what is called a generalized bit-reversal or CRT remapping. The standard trick is to precompute a lookup table of the remapped indices. For n1 and n2 up to a few hundred, this table is small enough to fit in cache and the remapping becomes a simple array indexing operation. Beyond that and you start spending more time on the remap than you save on the transform itself. Another thing that catches people out is the twiddle factor computation. In the radix-2 case the twiddles are powers of the primitive N-th root of unity, and you typically precompute them in a table or generate them on the fly using a recursive trigonometric identity. In the Striffler-Hamby case you have two different sets of roots of unity, one of order n1 and one of order n2, and the twiddles are products of elements from each set. If your field doesn't have primitive roots of the required orders, the whole thing breaks. This is not a bug in the algorithm, it is a constraint of the mathematics. Over the complex numbers this is rarely an issue, but if you are implementing this over GF(2^m) or some other finite field, you need to verify that the field contains roots of unity of both orders n1 and n2 before you even start coding.

I ran into a particularly annoying edge case once where I was working with a binary sequence generation problem. The period of the sequence was N = 2^10 - 1 = 1023, which factors as 3 * 11 * 31. All three factors are pairwise coprime, so the prime-factor method applies directly. The problem was that the standard floating-point FFT libraries weren't designed for this kind of factorization and would silently produce incorrect results due to rounding errors accumulating differently than expected over a non-power-of-two structure. I ended up writing a custom implementation that did the three separate transforms over the integers using a number-theoretic transform variant, and the results matched the expected periodic sequence exactly where the floating-point versions drifted by up to 0.003 in the higher frequency bins. That drift doesn't sound like much but when you are correlating the output against a known reference sequence for synchronization, it is enough to cause frame slips. For people who want to experiment with this, there isn't a single canonical open-source implementation that everyone uses. The most common path is to adapt a standard FFT library like FFTW and replace the radix-2 kernel with a prime-factor kernel, but that requires significant modification to the planner and the backend selection logic. A simpler entry point is to find existing implementations of the Good-Thomas prime-factor algorithm, which is essentially the same mathematical structure as what Striffler and Hamby described, and modify it for your specific use case. GNU Scientific Library has some primitives that overlap with this approach. Sober FFT implementations also exist in a few academic repositories, though the documentation is usually thin and the code tends to be written for research purposes rather than production use. The honest assessment is that the Striffler-Hamby approach is useful in a narrow set of situations. If your data length has a favorable coprime factorization and you are working in an environment where the overhead of twiddle tables and radix-2 butterflies is significant relative to your compute budget, it can give you a meaningful speedup. In most other cases, including the vast majority of general-purpose signal processing work, a well-optimized radix-2 or mixed-radix FFT will outperform it. The main value of understanding this algorithm is that it gives you a different way to think about the structure of the DFT, which turns out to be useful when you encounter situations where the standard assumptions don't hold, like when N is prime or when you are working over finite fields where the available roots of unity are constrained.

If you are trying to decode something like a pseudo-random sequence or analyze the spectral content of a periodic digital signal and the usual tools aren't giving you clean results, it is worth checking whether the data length has a factorization that the Striffler-Hamby method can exploit. Factor N, check coprimality, verify your field has the necessary roots of unity, and then decide whether the math is actually simpler than what you already have. Most of the time the answer is no, but when it is yes, it is genuinely yes and you won't get that kind of result from a standard library call.