Skip to content

Lab 14: Sweeping All Frequencies — You Just Built a DFT

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

You're about to write a DFT. You just don't know it yet.

Echo waving welcome Last lab you built a detector for one frequency. Today we do the least imaginative thing possible with it — run it at every frequency and keep the answers. That list is a spectrum, and the loop that makes it has a famous name. Time to transform!

What You'll Build

A complete Discrete Fourier Transform, written by you, in about twelve lines. Then you'll use it to unmix a two-tone signal and discover why spectra are always drawn as mirror images.

Learning Objectives

  • Assemble a DFT by looping the Lab 13 detector over every frequency
  • Explain what a frequency bin is and how wide it is
  • Compute a bin's centre frequency from its index
  • Interpret the real and imaginary parts of a spectrum
  • Describe why the upper half of a real signal's spectrum is redundant
  • Identify the DC bin and the Nyquist bin

Concepts Introduced

ID Concept
332 Frequency Sweep
333 Bin Index
334 Bin Center Frequency
335 Bin Width
336 Spectrum Array
337 Eight Point DFT By Hand
338 Complex Exponential
339 Real And Imaginary Parts
340 Spectrum Symmetry
341 Negative Frequencies
342 DC Bin
343 Nyquist Bin

Background

Bins: the frequencies we can ask about

Lab 13's detector works cleanly when the test frequency fits a whole number of cycles in the window. Those special frequencies are the bins:

1
bin width = sample rate / N

With N = 64 at 12,800 Hz that's 200 Hz per bin. Bin 0 is 0 Hz, bin 1 is 200 Hz, bin 3 is 600 Hz, and so on.

A DFT doesn't report every possible frequency — it reports these N evenly-spaced ones. Want finer detail? Use a bigger N. That's the whole tradeoff, and Lab 23 makes you live with it.

The algorithm, in full

1
2
3
4
5
6
for k in range(n):              # for every bin...
    re = im = 0.0
    for t in range(n):          # ...correlate against the signal
        angle = 2 * math.pi * k * t / n
        re += signal[t] * math.cos(angle)
        im -= signal[t] * math.sin(angle)

That inner loop is Lab 13's detector, unchanged. The outer loop is the only new idea.

Real, imaginary — and why the minus sign

Echo thinking re is the cosine score, im the sine score — the two measurements you needed in Lab 13 to beat the phase problem. Mathematicians bundle them into one complex number because the algebra gets tidier, but they're just your two correlations. The minus sign on im is a sign convention for the forward transform; flip it and you get the inverse.

Procedure

Step 1 — Run it

Open 14-building-a-dft.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
# Lab 14: Sweeping All Frequencies -- You Just Built a DFT
#
# Lab 13 gave you a detector for ONE frequency. This lab does something
# almost embarrassingly simple with it:
#
#     run it at every frequency and keep all the answers.
#
# That list of answers is a SPECTRUM. And the loop that produces it is the
# Discrete Fourier Transform. Not a simplified version. The actual thing.

import config
import math

RATE = config.SAMPLE_RATE      # 12800 Hz
N = 64                         # small, so pure Python stays quick
BIN = RATE / N                 # 200 Hz per bin


def sine(freq, amp=1.0, phase=0.0, n=N):
    return [amp * math.sin(2 * math.pi * freq * (i / RATE) + phase)
            for i in range(n)]


# =========================================================================
# The DFT. This is the entire algorithm.
# =========================================================================
def dft(signal):
    """Return (real, imag) lists -- one entry per frequency bin.

    For every bin k we correlate the signal against a cosine and a sine at
    that bin's frequency. Exactly the Lab 13 detector, run in a loop.
    """
    n = len(signal)
    real = []
    imag = []
    for k in range(n):                      # for every frequency...
        re = 0.0
        im = 0.0
        for t in range(n):                  # ...correlate against the signal
            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 magnitudes(real, imag):
    return [math.sqrt(real[k] * real[k] + imag[k] * imag[k])
            for k in range(len(real))]


def show_spectrum(mags, label, limit=None):
    print()
    print("--- %s ---" % label)
    limit = limit or (len(mags) // 2 + 1)
    peak = max(mags[:limit]) or 1.0
    print("%4s %9s  %s" % ("bin", "Hz", "magnitude"))
    for k in range(limit):
        bar = "#" * int(mags[k] / peak * 44)
        print("%4d %8.0f  %s" % (k, k * BIN, bar))


# =========================================================================
# PART 1 -- an 8-point DFT you can check by hand
# =========================================================================
print("=== PART 1: a tiny DFT, small enough to verify by hand ===")
tiny_n = 8
# A signal that goes up, down, up, down: the fastest wiggle 8 samples can hold.
tiny = [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0]
print("signal:", tiny)
print()
tre = []
tim = []
for k in range(tiny_n):
    re = im = 0.0
    for t in range(tiny_n):
        angle = 2 * math.pi * k * t / tiny_n
        re += tiny[t] * math.cos(angle)
        im -= tiny[t] * math.sin(angle)
    tre.append(re)
    tim.append(im)
    print("bin %d: real %+7.3f  imag %+7.3f  magnitude %6.3f"
          % (k, re, im, math.sqrt(re * re + im * im)))
print()
print("All the energy sits in bin 4 -- exactly half of 8. That is the")
print("fastest frequency 8 samples can represent: the Nyquist bin.")

# =========================================================================
# PART 2 -- a real spectrum
# =========================================================================
print()
print("=== PART 2: the spectrum of a single tone ===")
print("Bin width = %.0f Hz, so bin k covers k * %.0f Hz." % (BIN, BIN))
sig = sine(3 * BIN)                       # 600 Hz -- exactly bin 3
mags = magnitudes(*dft(sig))
show_spectrum(mags, "600 Hz tone (should peak at bin 3)")

# =========================================================================
# PART 3 -- two tones at once
# =========================================================================
print()
print("=== PART 3: two tones ===")
mixed = [a + b for a, b in zip(sine(2 * BIN), sine(5 * BIN, amp=0.5))]
mags = magnitudes(*dft(mixed))
show_spectrum(mags, "400 Hz (full) + 1000 Hz (half) -> bins 2 and 5")
print()
print("Two peaks, and the second is half the height of the first --")
print("exactly the recipe we mixed. The DFT unmixed it.")

# =========================================================================
# PART 4 -- the mirror
# =========================================================================
print()
print("=== PART 4: why we only ever plot half ===")
full = magnitudes(*dft(sine(3 * BIN)))
print("%4s %9s %10s" % ("bin", "Hz", "magnitude"))
for k in list(range(0, 6)) + list(range(N - 5, N)):
    print("%4d %8.0f %10.3f" % (k, k * BIN, full[k]))
print()
print("Bin %d matches bin 3, bin %d matches bin 2, and so on." % (N - 3, N - 2))
print("The top half is a mirror image of the bottom half. For a real-valued")
print("signal it carries no new information, so we throw it away and plot")
print("only bins 0 to %d." % (N // 2))
print()
print("bin 0        = DC (the average of the signal)")
print("bin %d       = Nyquist (%.0f Hz), the fastest we can represent"
      % (N // 2, RATE / 2))

Step 2 — Check the 8-point case by hand

Part 1 transforms [1, -1, 1, -1, 1, -1, 1, -1] — the fastest wiggle 8 samples can hold.

1
bin 4: real  +8.000  imag  -0.000  magnitude  8.000

Everything lands in bin 4, exactly half of 8. That's the Nyquist bin (Lab 9's limit, showing up as an array index). Every other bin is zero.

Small enough that you can grind through the arithmetic on paper if you want to — and it's worth doing once.

Step 3 — Unmix a chord

Part 3 mixes 400 Hz at full amplitude with 1,000 Hz at half, then transforms it:

1
2
   2      400  ############################################
   5     1000  #####################

Two peaks, at bins 2 and 5, and the second is half the height of the first. The DFT recovered the exact recipe — amplitudes and all — from a signal where the two tones were hopelessly tangled together.

That's the superpower. A wave that looked like noise is now a list of ingredients.

Step 4 — Meet the mirror

1
2
   3      600     32.000
  61    12200     32.000

Bin 61 is an identical copy of bin 3. Bin 62 mirrors bin 2. The top half of the spectrum is a reflection of the bottom.

For a real-valued signal it carries no new information, so we plot only bins 0 to N/2:

Bin Meaning
0 DC — the signal's average
1 … N/2−1 real frequencies
N/2 Nyquist — the fastest representable
N/2+1 … N−1 the mirror; discard

Step 5 — Predict, then measure

Prediction: with N = 64 at 12,800 Hz, which bin holds a 1,400 Hz tone? What if you changed N to 128?

Work it out, then change the code and check.

Expected Output

1
2
3
4
5
6
7
8
=== PART 1: a tiny DFT, small enough to verify by hand ===
signal: [1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0]

bin 4: real  +8.000  imag  -0.000  magnitude  8.000

--- 400 Hz (full) + 1000 Hz (half) -> bins 2 and 5 ---
   2      400  ############################################
   5     1000  #####################

Troubleshooting

Symptom Likely cause Fix
Peak lands between bins Frequency isn't a bin centre Use multiples of RATE/N; Lab 22 handles the rest
Energy smeared everywhere Same cause, worse This is spectral leakage — a preview of Lab 22
Program takes ages N too large for pure Python Keep N ≤ 128 here; Lab 16 measures exactly how bad it gets
Mirror bins missing Only plotting the lower half That's intentional in Parts 2 and 3

Challenges

  1. Three tones. Mix bins 1, 4 and 9 at different amplitudes. Do all three appear at the right heights?
  2. Zero the DC. Add a constant 0.5 to your signal. Which bin changes? Now subtract the mean before transforming — this is what Lab 7's DC removal was for.
  3. Between the bins. Feed in a 500 Hz tone with 200 Hz bins. It has nowhere to go. Look at what happens to bins 2 and 3, and write down what you see — you'll explain it in Lab 22.

Check Your Understanding

  1. What is the bin width for N = 256 at 12,800 Hz?
  2. Which bin holds a 2,000 Hz tone under those settings?
  3. What do the real and imaginary parts of a bin correspond to, in Lab 13's terms?
  4. Why is the upper half of a real signal's spectrum redundant?
  5. What frequency does bin 0 represent, and what does its value tell you?

You wrote a Fourier transform

Echo celebrating Not a toy version. The real Discrete Fourier Transform, in a dozen lines, from an idea you built yourself two labs ago. Next: prove it's actually correct — because looking right and being right are different things.

Next: Lab 15: Validating Your DFT | Previous: Lab 13