Skip to content

Lab 17: Divide and Conquer — From DFT to FFT

Time: ~50 minutes | Prerequisites: Lab 16 | Hardware: Pico 2 (no microphone needed)

The trick that changed signal processing

Echo waving welcome Your DFT is 530× too slow. The fix isn't a faster computer — it's noticing that most of that work is the same work done twice. Cooley and Tukey spotted it in 1965 and made real-time signal processing possible. Today you spot it too. Time to transform!

What You'll Build

A hand-rolled single split: transform the even samples and the odd samples separately, stitch them back together, and prove the answer is identical — for half the work.

Learning Objectives

  • Split a signal into even and odd samples and transform each half
  • Recombine two half-spectra with a twiddle factor
  • Explain how one complex multiply produces two output bins
  • Count the operations saved per split, and after full recursion
  • Justify why N must be a power of two

Concepts Introduced

ID Concept
363 Divide And Conquer
364 Even Odd Split
365 Recursive Decomposition
366 Subproblem
367 Recombination Step
368 Logarithmic Stages
369 Complexity Reduction
370 Power Of Two Constraint
371 Redundant Computation
372 Symmetry Exploitation

Background

The observation

Split your N samples into evens and odds. Transform each half separately, giving E[k] and O[k]. Then the full transform is:

1
2
X[k]       = E[k] + W^k · O[k]
X[k + N/2] = E[k] − W^k · O[k]

Stare at those two lines. They use the same E[k], the same O[k], and the same product W^k · O[k]. One multiplication, two output bins.

The direct DFT computes those two bins completely independently, redoing work it already had.

W^k — the twiddle factor

1
2
angle = -2 * math.pi * k / n
wr, wi = math.cos(angle), math.sin(angle)

It's a point on the unit circle — a rotation. It accounts for the odd samples being offset by one position from the evens. Lab 18 gives it a proper table.

Why stop at one split?

Echo thinking One split halves the work. But each half is itself a DFT — so split those too. And again. Keep going until each piece is a single sample, whose transform is just itself with no arithmetic at all. That's log₂(N) levels of splitting, and it's where N·log₂(N) comes from.

Procedure

Step 1 — Prove the split is exact

Open 17-divide-and-conquer.py and run it:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
# Lab 17: Divide and Conquer -- From DFT to FFT
#
# Last lab we found the DFT is 530x too slow, and that it recomputes the same
# angles thousands of times. That is a hint: there is REDUNDANT WORK in there.
#
# Here is the observation that unlocks everything. Split your samples into
# the even-numbered ones and the odd-numbered ones. Transform each half
# separately. Then the full transform can be rebuilt from those two halves:
#
#     X[k]         = E[k] + W^k * O[k]
#     X[k + N/2]   = E[k] - W^k * O[k]
#
# Look at the second line. It reuses E[k] and O[k] -- ALREADY COMPUTED for
# the first line. One multiplication buys you two output bins.
#
# Do that recursively and N^2 collapses to N log N.

import config
import math
import time


def dft(signal):
    """The brute-force DFT from Lab 14."""
    n = len(signal)
    real, imag = [], []
    for k in range(n):
        re = im = 0.0
        for t in range(n):
            angle = 2 * math.pi * k * t / n
            re += signal[t] * math.cos(angle)
            im -= signal[t] * math.sin(angle)
        real.append(re)
        imag.append(im)
    return real, imag


def dft_split(signal):
    """One level of divide and conquer, done by hand.

    Transform the even samples and the odd samples separately, then stitch
    the results together. The answer is identical to dft() -- but the two
    half-size DFTs cost 2*(N/2)^2 = N^2/2 instead of N^2.
    """
    n = len(signal)
    half = n // 2

    evens = [signal[i] for i in range(0, n, 2)]
    odds = [signal[i] for i in range(1, n, 2)]

    e_re, e_im = dft(evens)
    o_re, o_im = dft(odds)

    real = [0.0] * n
    imag = [0.0] * n
    for k in range(half):
        # W^k, the "twiddle factor" -- a rotation by -2*pi*k/n
        angle = -2 * math.pi * k / n
        wr = math.cos(angle)
        wi = math.sin(angle)

        # W^k * O[k]  (complex multiply)
        tr = wr * o_re[k] - wi * o_im[k]
        ti = wr * o_im[k] + wi * o_re[k]

        # Two outputs from one multiplication. THIS is the saving.
        real[k] = e_re[k] + tr
        imag[k] = e_im[k] + ti
        real[k + half] = e_re[k] - tr
        imag[k + half] = e_im[k] - ti

    return real, imag


def magnitudes(real, imag):
    return [math.sqrt(real[k] ** 2 + imag[k] ** 2) for k in range(len(real))]


# =========================================================================
# PART 1 -- prove the split gives the same answer
# =========================================================================
N = 32
signal = [math.sin(2 * math.pi * 3 * i / N) + 0.5 * math.sin(2 * math.pi * 7 * i / N)
          for i in range(N)]

direct = magnitudes(*dft(signal))
split = magnitudes(*dft_split(signal))

print("=== PART 1: does splitting change the answer? ===")
print("%4s %12s %12s %12s" % ("bin", "direct DFT", "split DFT", "difference"))
worst = 0.0
for k in range(N // 2 + 1):
    d = abs(direct[k] - split[k])
    if d > worst:
        worst = d
    if direct[k] > 0.01 or k < 9:
        print("%4d %12.4f %12.4f %12.2e" % (k, direct[k], split[k], d))
print()
print("largest difference anywhere: %.2e" % worst)
print("Identical (to float precision). The split is exact, not an approximation.")

# =========================================================================
# PART 2 -- count the work
# =========================================================================
print()
print("=== PART 2: how much work did we save? ===")
print()
print("%8s %14s %16s %10s" % ("N", "direct (N^2)", "one split", "saved"))
for n in (16, 32, 64, 128, 256, 512):
    direct_ops = n * n
    # two half-size DFTs, plus n/2 recombinations
    split_ops = 2 * (n // 2) ** 2 + n // 2
    print("%8d %14d %16d %9.0f%%"
          % (n, direct_ops, split_ops, 100 * (1 - split_ops / direct_ops)))

print()
print("One split saves about half the work. But why stop at one split?")

# =========================================================================
# PART 3 -- keep splitting
# =========================================================================
print()
print("=== PART 3: splitting all the way down ===")
print()
print("Each split halves the work again. Keep going until each piece is a")
print("single sample -- whose DFT is just itself, requiring no arithmetic.")
print()
print("%8s %10s %14s %16s %12s" % ("N", "splits", "direct N^2", "FFT N*log2(N)", "speedup"))
for n in (16, 64, 256, 512, 1024, 4096):
    stages = int(math.log(n, 2) + 0.5)
    fft_ops = n * stages
    print("%8d %10d %14d %16d %11.0fx"
          % (n, stages, n * n, fft_ops, (n * n) / fft_ops))

print()
print("At N=512 that is a %.0fx reduction in arithmetic." % (512 * 512 / (512 * 9)))
print("This is why the FFT matters. Same answer, a fraction of the work.")

# =========================================================================
# PART 4 -- why powers of two
# =========================================================================
print()
print("=== PART 4: why N must be a power of two ===")
print("Halving only works cleanly if the size keeps dividing by 2.")
print()
for n in (512, 500):
    print("N = %d:" % n, end=" ")
    x = n
    chain = []
    while x % 2 == 0 and x > 1:
        chain.append(x)
        x //= 2
    chain.append(x)
    print(" -> ".join(str(c) for c in chain),
          "  (clean)" if x == 1 else "  <- STUCK at %d, cannot halve" % x)

# =========================================================================
# PART 5 -- measure it
# =========================================================================
print()
print("=== PART 5: measured, not just counted ===")
n = 128
sig = [math.sin(2 * math.pi * 5 * i / n) for i in range(n)]

start = time.ticks_ms()
dft(sig)
direct_ms = time.ticks_diff(time.ticks_ms(), start)

start = time.ticks_ms()
dft_split(sig)
split_ms = time.ticks_diff(time.ticks_ms(), start)

print("N = %d" % n)
print("  direct DFT      : %5d ms" % direct_ms)
print("  ONE split       : %5d ms   (%.1fx faster)" % (split_ms, direct_ms / split_ms))
print()
print("And that is from a single split. Lab 20 does all %d of them."
      % int(math.log(512, 2) + 0.5))

1
2
3
4
5
 bin   direct DFT    split DFT   difference
   3      16.0000      16.0000     1.91e-06
   7       8.0000       8.0000     9.54e-07

largest difference anywhere: 2.39e-06

Identical to float precision. This is not an approximation — it's an exact restructuring of the same arithmetic.

Step 2 — Count what one split saves

1
2
       N   direct (N^2)        one split      saved
     512         262144           131328        50%

Two half-size DFTs cost 2·(N/2)² = N²/2. Half the work, same answer.

Step 3 — Recurse all the way down

1
2
3
       N     splits     direct N^2    FFT N*log2(N)      speedup
     512          9         262144             4608          57x
    4096         12       16777216            49152         341x

Nine splits for 512 samples. 57× less arithmetic — and the advantage grows with N, which is why the FFT matters more the bigger your problem gets.

Step 4 — Why powers of two

1
2
N = 512: 512 -> 256 -> 128 -> 64 -> 32 -> 16 -> 8 -> 4 -> 2 -> 1   (clean)
N = 500: 500 -> 250 -> 125   <- STUCK at 125, cannot halve

Halving only works if the size keeps dividing evenly. 512 splits nine times to reach 1; 500 jams at 125.

(Other FFT variants handle non-powers-of-two — mixed-radix, Bluestein's — but radix-2 is the one you can hold in your head, and 512 is a fine size for audio.)

Step 5 — Measure it

1
2
3
N = 128
  direct DFT      :  1213 ms
  ONE split       :   551 ms   (2.2x faster)

One split, measured on hardware, gives 2.2×. Lab 20 does all nine.

Expected Output

Your millisecond values will differ; the ~50% saving per split and the 57× count at N=512 should not.

Troubleshooting

Symptom Likely cause Fix
Split disagrees with direct Twiddle sign wrong The angle is -2πk/n, negative
Only the first half is right Missing the X[k+N/2] line Both outputs come from one product
Errors around 1e-6 Normal float32 behaviour See Lab 15
Crash on odd N Can't halve Use a power of two

Challenges

  1. Split twice. Apply the split to each half as well, giving four quarter-size DFTs. What fraction of the original work remains?
  2. Make it recursive. Write fft_recursive(signal) that calls itself until the input is length 1. Check it against the direct DFT.
  3. Where's the limit? The count says 57× at N=512, but Lab 20 will measure ~146×. Why might the measured gain exceed the arithmetic count? (Hint: what else did we stop doing?)

Check Your Understanding

  1. Write the two recombination equations from memory.
  2. Why does one complex multiply give two output bins?
  3. Where does the log₂(N) in N·log₂(N) come from?
  4. Why must N be a power of two for radix-2?
  5. Is the split an approximation? Justify your answer from the measured output.

You found the redundant work

Echo celebrating That's the whole insight — the rest is bookkeeping. Next lab handles the bookkeeping: what order the data ends up in, and how to stop recomputing the same angles.

Next: Lab 18: Bit Reversal and Twiddle Factors | Previous: Lab 16