Lab 18: Bit Reversal and Twiddle Factors

Time: ~45 minutes | Prerequisites: Lab 17 | Hardware: Pico 2 (no microphone needed)

The bookkeeping that makes it fit

Echo waving welcome Splitting works, but it leaves your data in a peculiar order and it keeps recomputing the same angles. Both have surprisingly elegant fixes — and one of them is genuinely beautiful. Let's tune in.

What You'll Build

Two tables that turn the recursive idea into fast in-place code: a bit-reversal permutation and a twiddle-factor lookup table. Plus a measurement of exactly what the table buys you.

Learning Objectives

Discover that repeated even/odd splitting produces bit-reversed order
Build a bit-reversal permutation table
Reorder an array in place using only swaps
Precompute twiddle factors as roots of unity
Measure the speedup from table lookup versus recomputation
Explain why hoisting work out of a loop is the most reliable optimization

Concepts Introduced

ID	Concept
373	Bit Reversal Permutation
374	Index Reversal
375	Permutation Table
376	Roots Of Unity
377	Twiddle Factor Table
378	Precomputation
379	Lookup Table
380	Loop Invariant Hoisting
381	In Place Reordering
382	Swap Operation
383	Interleaved Storage

Background

The strange order

Split [0..7] into evens and odds, then split again, and again:

start          : [0, 1, 2, 3, 4, 5, 6, 7]
after split 1  : [0, 2, 4, 6, 1, 3, 5, 7]
after split 2  : [0, 4, 2, 6, 1, 5, 3, 7]

Final order: 0, 4, 2, 6, 1, 5, 3, 7. Looks arbitrary. It is not.

It's bit reversal

Write each index in binary and read it backwards:

index	binary	reversed	value
0	000	000	0
1	001	100	4
2	010	010	2
3	011	110	6
4	100	001	1

0, 4, 2, 6, 1, 5, 3, 7 — exactly the split order.

That's a gift. Instead of recursively slicing lists (allocating, copying, and generally being slow), we reorder once at the start with a handful of swaps and then work in place forever after.

Why bit reversal, intuitively

Echo thinking Each split sorts by one bit: the first asks "is the lowest bit 0 or 1?" (even or odd). The next asks about the next bit up. After all the splits, the ordering is by lowest bit first, highest bit last — which is exactly the ordinary ordering with the bits reversed.

Roots of unity

The twiddle factors are evenly spaced points around the unit circle:

1	`W_N^k = cos(-2πk/N) + i·sin(-2πk/N)`

For an N-point FFT you need only N/2 of them. Compute them once; look them up forever.

Procedure

Step 1 — Watch the pattern emerge

Open 18-bit-reversal-twiddles.py and run it:

# Lab 18: Bit Reversal and Twiddle Factors
#
# Lab 17 showed that splitting into evens and odds saves half the work, and
# that doing it repeatedly saves almost everything. But repeated splitting
# leaves the data in a strange order -- and we would rather not shuffle
# lists around at every level.
#
# It turns out the final order has a beautiful description: each index is
# the ORIGINAL index with its bits written backwards. So instead of
# splitting the data over and over, we reorder ONCE at the start and then
# work in place.
#
# The other half of this lab: stop calling sin() and cos() millions of times.
# There are only a few distinct angles. Compute them once, keep them in a
# table, look them up. That is a TWIDDLE FACTOR TABLE.

import config
import math
import time


# =========================================================================
# PART 1 -- where the strange order comes from
# =========================================================================
print("=== PART 1: what repeated splitting does to the order ===")
N = 8
order = list(range(N))
print("start                :", order)

level = 1
current = [order]
while len(current[0]) > 1:
    nxt = []
    for group in current:
        # MicroPython has no slice steps, so pick out evens and odds by hand
        nxt.append([group[i] for i in range(0, len(group), 2)])   # evens
        nxt.append([group[i] for i in range(1, len(group), 2)])   # odds
    current = nxt
    flat = [x for g in current for x in g]
    print("after split %d        : %s" % (level, flat))
    level += 1

final = [x for g in current for x in g]

# =========================================================================
# PART 2 -- that order IS bit reversal
# =========================================================================
print()
print("=== PART 2: the pattern ===")
bits = int(math.log(N, 2) + 0.5)
print("%6s %10s %12s %8s" % ("index", "binary", "reversed", "value"))
computed = []
for i in range(N):
    b = "".join(str((i >> (bits - 1 - p)) & 1) for p in range(bits))
    r = "".join(str((i >> p) & 1) for p in range(bits))   # same bits, backwards
    v = int(r, 2)
    computed.append(v)
    print("%6d %10s %12s %8d" % (i, b, r, v))

print()
print("split order    :", final)
print("bit-reversed   :", computed)
print("same?          :", final == computed)
print()
print("Write the index in binary, read it backwards, and you have the")
print("position it belongs in. No recursion, no list shuffling -- just a")
print("permutation we can do once, in place.")


# =========================================================================
# PART 3 -- doing the reorder in place
# =========================================================================
def bit_reverse_table(n):
    """Precompute where each element belongs."""
    bits = 0
    while (1 << bits) < n:
        bits += 1
    table = []
    for i in range(n):
        r = 0
        x = i
        for _ in range(bits):
            r = (r << 1) | (x & 1)
            x >>= 1
        table.append(r)
    return table


def bit_reverse_in_place(data, table):
    """Swap elements into bit-reversed order without copying the list."""
    swaps = 0
    for i in range(len(data)):
        j = table[i]
        if j > i:                      # only swap each pair once
            data[i], data[j] = data[j], data[i]
            swaps += 1
    return swaps


print()
print("=== PART 3: reordering in place ===")
data = list(range(N))
table = bit_reverse_table(N)
swaps = bit_reverse_in_place(data, table)
print("after in-place reorder:", data)
print("swaps performed       : %d (not %d -- half the entries stay put or" % (swaps, N))
print("                        pair up with one already moved)")

for n in (64, 512):
    t = bit_reverse_table(n)
    s = sum(1 for i in range(n) if t[i] > i)
    print("N=%4d needs %3d swaps for %4d elements (%.0f%%)" % (n, s, n, 100 * s / n))


# =========================================================================
# PART 4 -- twiddle factors, computed once
# =========================================================================
print()
print("=== PART 4: the twiddle factor table ===")
print()
print("A twiddle factor is a point on the unit circle:")
print("    W_N^k = cos(-2*pi*k/N) + i*sin(-2*pi*k/N)")
print()

n = 8
print("The %d twiddle factors for N=%d:" % (n // 2, n))
print("%4s %10s %10s %10s" % ("k", "angle", "cos", "sin"))
for k in range(n // 2):
    angle = -2 * math.pi * k / n
    print("%4d %10.4f %10.4f %10.4f" % (k, angle, math.cos(angle), math.sin(angle)))
print()
print("They march evenly around the circle -- these are the ROOTS OF UNITY.")
print("For an N-point FFT you need only N/2 of them, ever.")


def twiddle_table(n):
    """cos and sin for every twiddle we will need, computed once."""
    half = n // 2
    tw_re = [0.0] * half
    tw_im = [0.0] * half
    for k in range(half):
        angle = -2 * math.pi * k / n
        tw_re[k] = math.cos(angle)
        tw_im[k] = math.sin(angle)
    return tw_re, tw_im


# =========================================================================
# PART 5 -- how much does a lookup table buy?
# =========================================================================
print()
print("=== PART 5: table lookup vs recomputing ===")
n = 256
reps = 20

start = time.ticks_ms()
for _ in range(reps):
    for k in range(n // 2):
        angle = -2 * math.pi * k / n
        c = math.cos(angle)
        s = math.sin(angle)
compute_ms = time.ticks_diff(time.ticks_ms(), start)

tw_re, tw_im = twiddle_table(n)
start = time.ticks_ms()
for _ in range(reps):
    for k in range(n // 2):
        c = tw_re[k]
        s = tw_im[k]
lookup_ms = time.ticks_diff(time.ticks_ms(), start)

print("N = %d, %d passes over %d twiddles" % (n, reps, n // 2))
print("  recomputing sin/cos : %5d ms" % compute_ms)
print("  table lookup        : %5d ms" % lookup_ms)
if lookup_ms > 0:
    print("  speedup             : %.1fx" % (compute_ms / lookup_ms))
print()
print("Building the table costs one pass. After that every lookup is free.")
print("Move work OUT of the inner loop -- this is the single most reliable")
print("optimization there is, and it costs only memory.")
print()
print("Bonus: remember Lab 15's precision problem? Recomputing")
print("2*pi*k*t/N gave angles up to ~390 radians, which single-precision")
print("floats handle badly. A table only ever holds angles in one turn of")
print("the circle. Faster AND more accurate.")

Part 1 splits repeatedly; Part 2 shows the same order arising from reversing bits. The program compares them and prints same? : True.

Step 2 — Reorder with swaps

after in-place reorder: [0, 4, 2, 6, 1, 5, 3, 7]
swaps performed       : 2
N= 512 needs 240 swaps for  512 elements (47%)

Only about 47% of positions need a swap — the rest either map to themselves (0 and 7 above) or have already been handled by their partner. The if j > i test is what stops us swapping every pair twice and undoing our own work.

Swap each pair exactly once

Echo warning Drop the if j > i and the loop swaps (1,4) and later (4,1), putting everything back where it started. The array ends up unchanged and your FFT produces confident nonsense. A classic, silent, deeply annoying bug.

Step 3 — Build the twiddle table

   k      angle        cos        sin
   0.0000     1.0000     0.0000
  -0.7854     0.7071    -0.7071
  -1.5708     0.0000    -1.0000
  -2.3562    -0.7071    -0.7071

Four twiddles for an 8-point FFT, marching a quarter-turn at a time around the circle.

Step 4 — Measure the payoff

N = 256, 20 passes over 128 twiddles
  recomputing sin/cos :   131 ms
  table lookup        :    24 ms
  speedup             : 5.5x

5.5× on the twiddle work alone, for the cost of one setup pass and 128 floats of memory.

Step 5 — The bonus you already met

Remember Lab 15's precision problem? Recomputing 2πkt/N gave angles up to ~390 radians, which single-precision floats handle badly.

A twiddle table only ever holds angles within one turn of the circle. So precomputation is faster and more accurate. That's rare — most optimizations trade one for the other.

Expected Output

split order    : [0, 4, 2, 6, 1, 5, 3, 7]
bit-reversed   : [0, 4, 2, 6, 1, 5, 3, 7]
same?          : True

  recomputing sin/cos :   131 ms
  table lookup        :    24 ms
  speedup             : 5.5x

Troubleshooting

Symptom	Likely cause	Fix
Array unchanged after reorder	Missing `if j > i`	Every pair swapped twice cancels out
Reversed values look wrong	Wrong bit count	8 elements needs 3 bits, 512 needs 9
`NotImplementedError` on slices	MicroPython has no slice steps	Use an explicit `range()` loop
Table lookup no faster	Table built inside the timed loop	Build it once, outside

Challenges

Reversal without a table. Write a function that bit-reverses an index directly with shifts and masks. Compare its speed against the table for N = 512.
Memory cost. How many bytes does a 512-point twiddle table use as float32? As float64? Set that against the 485 KB of RAM from Lab 3.
Interleaved storage. Store twiddles as [re0, im0, re1, im1, …] instead of two lists. Why might that be faster on real hardware? (Lab 30 answers this properly.)

Check Your Understanding

What order does repeated even/odd splitting leave the data in?
Bit-reverse index 5 in a 16-element array.
Why does the swap loop need if j > i?
How many twiddle factors does a 512-point FFT need?
Name two separate benefits of precomputing the twiddle table.

The bookkeeping is done

Echo celebrating Reorder once, look up the twiddles, work in place. All that's left is the operation at the centre of it all — and it's only four multiplies.

Next: Lab 19: The Butterfly | Previous: Lab 17