Code Lab: Fourier Analysis of Musical Sounds

Every musical sound is a sum of pure sine waves at different frequencies and amplitudes. Joseph Fourier proved this in 1822, and his insight remains the single most powerful tool for understanding timbre. In this lab you will use Python's Fast Fourier Transform (FFT) to decompose complex waveforms into their frequency components, compare the spectra of different instrument-like tones, and build a reusable function for identifying harmonics.

Prerequisites: numpy, matplotlib, scipy. Install with pip install numpy matplotlib scipy.

Setting Up the Environment

We begin by importing libraries and defining global constants that match CD-quality audio.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

SAMPLE_RATE = 44100   # samples per second (CD quality)
DURATION = 1.0        # one second of audio for good frequency resolution
N_SAMPLES = int(SAMPLE_RATE * DURATION)
t = np.linspace(0, DURATION, N_SAMPLES, endpoint=False)

The frequency resolution of an FFT equals SAMPLE_RATE / N_SAMPLES. With these settings each frequency bin is 1 Hz wide, which is more than enough to distinguish musical pitches.

Generating Three Classic Waveforms

Different instruments produce different waveform shapes. We will synthesize three archetypes: a sine wave (flute-like), a square wave (clarinet-like, odd harmonics only), and a sawtooth wave (bowed-string-like, all harmonics). Each is built from its Fourier series so the result is band-limited and alias-free.

def make_sine(freq, t):
    """Pure tone -- single frequency, no harmonics."""
    return np.sin(2 * np.pi * freq * t)

def make_square(freq, t, n_harmonics=30):
    """Square wave from odd-harmonic Fourier series: sum of sin(n*f)/n for odd n."""
    wave = np.zeros_like(t)
    for n in range(1, n_harmonics * 2, 2):  # 1, 3, 5, ...
        wave += (1.0 / n) * np.sin(2 * np.pi * n * freq * t)
    return (4.0 / np.pi) * wave

def make_sawtooth(freq, t, n_harmonics=30):
    """Sawtooth wave: all harmonics with amplitude 1/n."""
    wave = np.zeros_like(t)
    for n in range(1, n_harmonics + 1):
        wave += ((-1) ** (n + 1)) * (1.0 / n) * np.sin(2 * np.pi * n * freq * t)
    return (2.0 / np.pi) * wave

FREQ = 440.0  # A4 -- concert pitch

sine_wave = make_sine(FREQ, t)
square_wave = make_square(FREQ, t)
sawtooth_wave = make_sawtooth(FREQ, t)

Computing and Plotting the Frequency Spectrum

The function below computes a one-sided magnitude spectrum and optionally labels the strongest peaks (harmonics). This is the core analysis tool you will reuse throughout the course.

def plot_frequency_spectrum(signal, sample_rate=SAMPLE_RATE,
                            title="Frequency Spectrum",
                            max_freq=5000, peak_threshold=0.05):
    """
    Compute the FFT of *signal*, plot the magnitude spectrum up to
    *max_freq* Hz, and annotate peaks above *peak_threshold* (relative
    to the maximum).

    Returns
    -------
    freqs : numpy array -- frequency axis (Hz)
    magnitudes : numpy array -- normalized magnitudes
    """
    n = len(signal)
    fft_vals = fft(signal)
    magnitudes = np.abs(fft_vals) / n
    freqs = fftfreq(n, d=1.0 / sample_rate)

    # Keep positive frequencies only
    pos = freqs >= 0
    freqs = freqs[pos]
    magnitudes = magnitudes[pos]
    magnitudes[1:-1] *= 2  # compensate for discarded negative side

    # Plot
    fig, ax = plt.subplots(figsize=(11, 4))
    mask = freqs <= max_freq
    ax.plot(freqs[mask], magnitudes[mask], color="#2E86AB", linewidth=1.2)
    ax.set_title(title, fontsize=13)
    ax.set_xlabel("Frequency (Hz)")
    ax.set_ylabel("Amplitude")
    ax.grid(True, alpha=0.3)

    # Identify and annotate peaks
    threshold = peak_threshold * np.max(magnitudes)
    peak_indices = np.where((magnitudes > threshold) & (freqs <= max_freq))[0]
    # Keep only local maxima
    for idx in peak_indices:
        if idx == 0 or idx == len(magnitudes) - 1:
            continue
        if magnitudes[idx] > magnitudes[idx - 1] and magnitudes[idx] > magnitudes[idx + 1]:
            harmonic_num = round(freqs[idx] / freqs[peak_indices[0]])
            if harmonic_num < 1:
                continue
            ax.annotate(f"H{harmonic_num}\n{freqs[idx]:.0f} Hz",
                        xy=(freqs[idx], magnitudes[idx]),
                        xytext=(freqs[idx] + 60, magnitudes[idx]),
                        fontsize=8, color="#E84855",
                        arrowprops=dict(arrowstyle="->", color="#E84855"))

    plt.tight_layout()
    plt.show()
    return freqs, magnitudes

Analyzing the Three Waveforms

Now call the function on each waveform and observe how the harmonic content differs.

plot_frequency_spectrum(sine_wave,
                        title="Sine Wave Spectrum (A4 = 440 Hz) -- single peak")

plot_frequency_spectrum(square_wave,
                        title="Square Wave Spectrum -- odd harmonics only (1, 3, 5, ...)")

plot_frequency_spectrum(sawtooth_wave,
                        title="Sawtooth Wave Spectrum -- all harmonics (1, 2, 3, ...)")

What to notice:

  • The sine wave has a single spike at 440 Hz -- it is the definition of a "pure" tone.
  • The square wave shows peaks only at odd multiples of 440 Hz (440, 1320, 2200, ...), each falling off as 1/n. This is why a clarinet, whose cylindrical bore reinforces odd harmonics, sounds "hollow."
  • The sawtooth wave contains every harmonic (440, 880, 1320, ...) with amplitudes falling as 1/n. A bowed violin string vibrates in a roughly sawtooth pattern, which is why strings sound rich and bright.

Comparing All Three on One Plot

A direct overlay makes the differences vivid.

fig, ax = plt.subplots(figsize=(12, 5))
for signal, label, color in [
    (sine_wave,     "Sine (pure tone)",  "#2196F3"),
    (square_wave,   "Square (odd harmonics)", "#FF5722"),
    (sawtooth_wave, "Sawtooth (all harmonics)", "#4CAF50"),
]:
    n = len(signal)
    mags = 2.0 * np.abs(fft(signal))[:n // 2] / n
    freqs = fftfreq(n, 1.0 / SAMPLE_RATE)[:n // 2]
    ax.plot(freqs, mags, label=label, color=color, alpha=0.8)

ax.set_xlim(0, 5000)
ax.set_xlabel("Frequency (Hz)")
ax.set_ylabel("Amplitude")
ax.set_title("Spectral Comparison: Three Waveforms at 440 Hz", fontsize=13)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Try It Yourself

  1. Investigate the triangle wave. Write a make_triangle function using the Fourier series for a triangle wave (odd harmonics with amplitudes falling as 1/n squared, alternating sign). Plot its spectrum and compare it to the square wave. Why does the triangle wave sound softer even though both contain only odd harmonics?

  2. Build a "mystery" timbre. Create a waveform by summing harmonics 1, 3, 5, and 7 with amplitudes 1.0, 0.8, 0.3, and 0.1 respectively. Plot its spectrum and its time-domain shape. Does it look more like a square wave or a triangle wave? What would it sound like?

  3. Explore frequency resolution. Shorten DURATION to 0.05 seconds (50 ms) and recompute the FFT of the sawtooth wave. What happens to the spectral peaks? Why does a shorter analysis window reduce your ability to distinguish nearby frequencies? (Hint: compute SAMPLE_RATE / N_SAMPLES for both durations.)