Build a DFT Frequency Analyzer

Implement the Discrete Fourier Transform from scratch and analyze audio signals.

Introduction

30 Second Summary

Every time Shazam identifies a song in seconds, it starts by breaking the audio into individual frequencies. This frequency decomposition is the math behind every audio tool, from voice assistants to music production software.

In this project, you will build a frequency analyzer from scratch using only Python and numpy. Your analyzer generates audio signals, plays them through your speakers, and decomposes them into individual frequency components.

What You'll Build

A command-line frequency analyzer that implements the Discrete Fourier Transform, verifies it against numpy's built-in FFT, and identifies the individual frequencies hidden inside a composite audio signal.

By the end of this project, you'll have:

  • A working DFT implementation built from scratch that produces identical results to numpy's built-in FFT.
  • A frequency analyzer that decomposes a multi-frequency signal and displays exact frequencies and amplitudes in a magnitude spectrum.
  • A code-based understanding of complex numbers, Euler's formula, and how rotating phasors detect frequencies in a signal.
  • Secret Mission: Implement the Inverse DFT and prove you can perfectly reconstruct the original signal from its frequency components.

Are there any prerequisites?

This project assumes basic Python and numpy experience. No prior knowledge of complex numbers or signal processing is required.

Before We Start

Before We Start...

Set Up Your Project

Before you can analyze frequencies, you need a workspace and the right tools installed. This project uses three Python packages that handle math, plotting, and audio playback.

In this step, you'll create a project folder, install the packages, and set up the main script file that will hold all your code.

In this step, get ready to:
  • Create your project folder and open it in VS Code.
  • Install numpy, matplotlib, and sounddevice.
  • Create the dft_analyzer.py file with imports.
Create the project folder and install packages

Every project starts with a dedicated folder and the right packages. You'll create a folder on your Desktop so it's easy to find.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type VS Code and press Enter to open it.
  • Open the VS Code terminal by pressing Ctrl+` (backtick key, next to the 1 key) on both macOS and Windows.
  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder by running this command:
mkdir dft-from-scratch
  • Confirm the folder was created by running:
ls

You should see dft-from-scratch listed in the output.

  • Now open the project folder in VS Code. Click File in the top menu bar.
  • Click Open Folder.
  • Navigate to your Desktop, select the dft-from-scratch folder, then click Open.

Seeing a Restricted Mode banner?

VS Code may show a banner at the top saying the folder is in Restricted Mode. Click Trust to allow full functionality. This is safe because you just created the folder yourself.

  • Open the VS Code terminal again by pressing Ctrl+` (backtick).

The terminal should now show dft-from-scratch as your current directory.

  • Install all three project packages by running this command:
pip install numpy==2.5.0 matplotlib==3.11.0 sounddevice==0.5.5

You should see download progress bars followed by a message confirming that the packages were successfully installed.

What do these packages do?

  • numpy provides fast array math. You'll use it for every calculation in this project.
  • matplotlib handles plotting. You'll use it to visualize waveforms and frequency spectrums.
  • sounddevice plays audio through your speakers so you can hear the signals you generate.

Seeing 'command not found' for pip?

Try running the command with pip3 instead of pip. On some systems, Python 3's package manager is named pip3.

help me fix a pip command not found error when installing numpy, matplotlib, and sounddevice

Verify the installations work
  • Check that numpy installed correctly by running:
python -c "import numpy; print(numpy.__version__)"
  • Check matplotlib by running:
python -c "import matplotlib; print(matplotlib.__version__)"
  • Check sounddevice by running:
python -c "import sounddevice; print(sounddevice.__version__)"

You should see 2.5.0, 3.11.0, and 0.5.5 printed without any errors.

Seeing 'command not found' for python?

Try using python3 instead of python. On macOS especially, the python command may not be linked to your Python 3.12 installation.

If a version number appears but it is lower than 3.12, your system may be using an older Python. Try running python3.12 directly to confirm your installation is available.

help me fix my Python version or command not found error when verifying my package installations

Create the project file

With the packages installed, you need a Python file to write your code in. This single file will hold your entire frequency analyzer by the end of the project.

  • Create a new file by clicking File in the top menu bar.
  • Click New File.
  • Name the file dft_analyzer.py and press Enter.
  • Add the following imports to the top of dft_analyzer.py:
import numpy as np
import matplotlib.pyplot as plt
import sounddevice as sd

What do these imports do?

  • import numpy as np loads numpy and gives it the short alias np. You'll use this for all array math, from generating sine waves to building the DFT matrix.
  • import matplotlib.pyplot as plt loads matplotlib's plotting module. Every visualization in this project uses plt to create figures and save them as images.
  • import sounddevice as sd loads the audio playback library. It sends numpy arrays to your speakers so you can hear the signals you create.
  • Save dft_analyzer.py by pressing Cmd+S (macOS) or Ctrl+S (Windows).
  • Verify the file runs without errors by running this command in the terminal:
python dft_analyzer.py

The script should exit cleanly with no output and no errors. This confirms all three imports resolved correctly.

Seeing a ModuleNotFoundError?

Make sure you ran the pip install command from earlier in the same terminal. If you are using python3 to run the script, make sure you installed the packages with pip3 so they go to the same Python installation.

Check that the file is saved. VS Code shows a dot on the tab name when there are unsaved changes.

help me fix a ModuleNotFoundError when running my dft_analyzer.py script

✔️ Awesome, I've got everything!

Great. Double check you've saved your file before moving on.

ⓧ I'd like to double check the full code

Here is the complete dft_analyzer.py file so far. Make sure yours matches exactly.

import numpy as np
import matplotlib.pyplot as plt
import sounddevice as sd

Your workspace is ready and all three packages are installed. Next, you'll use numpy to generate sine waves from the math formula and hear them through your speakers.

Generate and Hear Sine Waves

Your development environment is set up with numpy, matplotlib, and sounddevice. Now it's time to turn math into sound.

Every sound you hear is a vibration at some frequency. In this step, you will generate sine waves from the mathematical formula, plot their waveforms, and hear them through your speakers.

In this step, get ready to:
  • Write reusable functions to generate and play sine waves.
  • Generate a 440 Hz tone and plot its waveform.
  • Play sine waves through your speakers and compare two frequencies.
Write functions to generate and play sine waves

A sine wave is defined by the formula signal(t) = sin(2π × frequency × t). To generate one digitally, you sample this formula at regular time intervals, creating an array of numbers that represents the wave.

  • In dft_analyzer.py, add the following code below the import sounddevice as sd line:
# ============================================================
# STEP 2: Generate and Hear Sine Waves
# ============================================================

def generate_sine(freq, duration, sample_rate=44100):
    """Generate a sine wave at a given frequency.

    A sine wave is the simplest possible sound: a pure tone.
    The formula is: signal(t) = sin(2 * pi * freq * t)
    """
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * freq * t)
    return t, signal

What does this code do?

  • np.linspace creates an array of evenly spaced time values from 0 to duration. At a sample rate of 44100, a 1-second signal produces 44,100 time points.
  • np.sin(2 * np.pi * freq * t) evaluates the sine function at every time point. The result is a numpy array where each value is the wave's amplitude at that instant.
  • endpoint=False excludes the final value so the samples line up correctly when repeating the wave.
  • Save dft_analyzer.py.
  • You won't see any output yet. The next chunk adds a function to play audio through your speakers.

Generating the wave is half the story. You also need a function to send that array of numbers to your speakers as audible sound.

  • Add the play_signal function directly below generate_sine:

def play_signal(signal, sample_rate=44100):
    """Play a numpy array as audio through speakers."""
    audio = signal / np.max(np.abs(signal))
    sd.play(audio.astype(np.float32), sample_rate)
    sd.wait()

What does this code do?

  • signal / np.max(np.abs(signal)) normalizes the signal so its values stay between -1 and 1. Audio hardware expects values in this range.
  • .astype(np.float32) converts the array to 32-bit float, which is the format sounddevice sends to your audio driver.
  • sd.wait() blocks the script until playback finishes, so the program doesn't exit before the sound plays.
  • Save dft_analyzer.py.
  • You won't see output yet. The next substep calls both functions to generate and plot a tone.
Generate a 440 Hz tone and visualize it

The note A4 on a piano vibrates at exactly 440 Hz. That means the air pressure oscillates 440 times per second. Plotting this wave lets you see what a pure frequency looks like before you start decomposing complex sounds.

  • Add the following code below the play_signal function to generate a 440 Hz sine wave and plot its waveform:

# Generate a 440 Hz sine wave (the note A4)
sample_rate = 44100
duration = 1.0
t, signal_440 = generate_sine(440, duration, sample_rate)

# Plot the first 10 milliseconds to see the waveform
samples_10ms = int(sample_rate * 0.01)
plt.figure(figsize=(10, 4))
plt.plot(t[:samples_10ms] * 1000, signal_440[:samples_10ms])
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.title("440 Hz Sine Wave (A4 note) - first 10 ms")
plt.grid(True)
plt.tight_layout()
plt.savefig("sine_wave_440.png", dpi=100)
plt.show()

What does this code do?

  • generate_sine(440, duration, sample_rate) creates 44,100 samples of a 440 Hz sine wave lasting 1 second.
  • samples_10ms = int(sample_rate * 0.01) calculates how many samples fit in 10 milliseconds. At 44100 Hz, that is 441 samples.
  • The plot only shows the first 10 ms because a full 1-second plot would cram 440 cycles into one figure. 10 ms gives you approximately 4.4 clean cycles to inspect.
  • plt.savefig("sine_wave_440.png", dpi=100) saves the plot as an image file before plt.show() displays it in a window.
  • Save dft_analyzer.py.
  • Run the script by executing this command in your terminal:
python dft_analyzer.py
  • A matplotlib window should open showing a clean sinusoidal wave oscillating between -1 and 1, with approximately 4.4 complete cycles.
  • Close the plot window. A file called sine_wave_440.png should now appear in your dft-from-scratch folder.

Don't see a plot window?

  • Check that your indentation is consistent. Python functions require all lines inside them to be indented with the same number of spaces.
  • Make sure plt.show() is at the bottom of the plotting section, not inside a function. It needs to be at the top level of your script.
  • If the plot window opens but appears blank, try resizing it. Some systems need a manual resize to trigger the initial render.

Still stuck? Help me debug why my matplotlib plot isn't showing when I run dft_analyzer.py.

Play the sine waves through your speakers

Seeing the waveform is one thing. Hearing it makes the connection between frequency and pitch tangible. A higher frequency means the wave oscillates faster, and your ears perceive that as a higher-pitched sound.

  • Add the following code at the bottom of dft_analyzer.py, below the plt.show() line, to play both a 440 Hz and a 220 Hz tone:

# Play the tones through your speakers
print("Playing 440 Hz tone (A4 note)...")
play_signal(signal_440, sample_rate)

_, signal_220 = generate_sine(220, duration, sample_rate)
print("Playing 220 Hz tone (A3 note - one octave lower)...")
play_signal(signal_220, sample_rate)

What does this code do?

  • play_signal(signal_440, sample_rate) sends the 440 Hz array to your speakers. The note A4 is the standard tuning pitch used by orchestras.
  • The second call generates and plays a 220 Hz tone. 220 Hz is exactly half of 440 Hz, which means it sounds one octave lower.
  • The _ in _, signal_220 discards the time array because you already have one from the 440 Hz call.
  • Save dft_analyzer.py.
  • Run the script again by executing:
python dft_analyzer.py
  • Close the plot window when it appears. After you close it, the terminal should print Playing 440 Hz tone (A4 note)... and you should hear a steady tone.
  • After the first tone finishes, you should hear a second, noticeably lower tone as the 220 Hz signal plays.

Can't hear any audio?

  • Check that your system volume is turned up and not muted.
  • Make sure no other application is exclusively holding the audio device (video calls, music players).
  • If audio does not work on your system, this is not a blocker. The core learning in this project comes from the code and visualizations, not from hearing the tones.

Need help? Help me troubleshoot why sounddevice is not playing audio on my system.

Why does 220 Hz sound lower than 440 Hz?

Frequency measures how many times the wave oscillates per second. A 440 Hz wave completes 440 cycles in one second. A 220 Hz wave completes only 220 cycles in the same time.

Your ear perceives faster oscillation as higher pitch. When you halve the frequency, the pitch drops by exactly one octave.

✔️ Awesome, I've got everything!

Great. Double check you have saved dft_analyzer.py before moving on.

ⓧ I'd like to double check the full code

Here is the complete dft_analyzer.py file as it should look at the end of this step.

import numpy as np
import matplotlib.pyplot as plt
import sounddevice as sd

# ============================================================
# STEP 2: Generate and Hear Sine Waves
# ============================================================

def generate_sine(freq, duration, sample_rate=44100):
    """Generate a sine wave at a given frequency.

    A sine wave is the simplest possible sound: a pure tone.
    The formula is: signal(t) = sin(2 * pi * freq * t)
    """
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * freq * t)
    return t, signal


def play_signal(signal, sample_rate=44100):
    """Play a numpy array as audio through speakers."""
    audio = signal / np.max(np.abs(signal))
    sd.play(audio.astype(np.float32), sample_rate)
    sd.wait()


# Generate a 440 Hz sine wave (the note A4)
sample_rate = 44100
duration = 1.0
t, signal_440 = generate_sine(440, duration, sample_rate)

# Plot the first 10 milliseconds to see the waveform
samples_10ms = int(sample_rate * 0.01)
plt.figure(figsize=(10, 4))
plt.plot(t[:samples_10ms] * 1000, signal_440[:samples_10ms])
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.title("440 Hz Sine Wave (A4 note) - first 10 ms")
plt.grid(True)
plt.tight_layout()
plt.savefig("sine_wave_440.png", dpi=100)
plt.show()

# Play the tones through your speakers
print("Playing 440 Hz tone (A4 note)...")
play_signal(signal_440, sample_rate)

_, signal_220 = generate_sine(220, duration, sample_rate)
print("Playing 220 Hz tone (A3 note - one octave lower)...")
play_signal(signal_220, sample_rate)

You have generated sine waves from the formula and heard the difference between two frequencies. Next up, you will explore the math that makes frequency detection possible: complex numbers and Euler's formula.

Visualize Complex Numbers and Euler's Formula

In the previous step you generated pure sine waves and heard them as tones. A single sine wave is simple, but real audio is made of many frequencies mixed together. To pull those individual frequencies apart, you need a mathematical tool that can "ask" a signal how much of each frequency it contains.

That tool is the Discrete Fourier Transform (DFT), and it is built on a surprisingly elegant idea: complex numbers spinning on a circle. In this step, you will visualize Euler's formula and rotating phasors to see exactly how rotation encodes frequency.

In this step, get ready to:
  • Visualize Euler's formula as a unit circle on the complex plane.
  • Animate a rotating phasor to see how spinning speed maps to frequency.
  • Connect phasor rotation to how the DFT detects frequencies in a signal.
Visualize Euler's formula on the unit circle

Euler's formula says that e^(i*theta) equals cos(theta) + i*sin(theta). In plain terms, raising e to an imaginary power traces a circle. The real part of that circle is cosine, and the imaginary part is sine.

  • In dft_analyzer.py, scroll to the bottom of the file, below the line that says:
print("Playing 220 Hz tone (A3 note - one octave lower)...")
play_signal(signal_220, sample_rate)
  • Add the following section header and function below that line:
# ============================================================
# STEP 3: Complex Numbers and Euler's Formula
# ============================================================

def plot_euler_formula():
    """Visualize Euler's formula: e^(i*theta) = cos(theta) + i*sin(theta)."""
    theta = np.linspace(0, 2 * np.pi, 200)
    z = np.exp(1j * theta)

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # Plot 1: The unit circle on the complex plane
    axes[0].plot(z.real, z.imag, "b-", linewidth=2)
    angles = [0, np.pi / 4, np.pi / 2, np.pi, 3 * np.pi / 2]
    labels = ["0", "pi/4", "pi/2", "pi", "3pi/2"]
    for a, label in zip(angles, labels):
        point = np.exp(1j * a)
        axes[0].plot(point.real, point.imag, "ro", markersize=8)
        axes[0].annotate(
            f"theta={label}",
            (point.real, point.imag),
            textcoords="offset points",
            xytext=(10, 10),
            fontsize=8,
        )
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_xlabel("Real (cosine)")
    axes[0].set_ylabel("Imaginary (sine)")
    axes[0].set_title("Unit Circle: e^(i*theta)")

    # Plot 2: Real part = cos(theta)
    axes[1].plot(theta, z.real, "b-", linewidth=2)
    axes[1].set_xlabel("theta (radians)")
    axes[1].set_ylabel("Real part")
    axes[1].set_title("cos(theta) = Real part of e^(i*theta)")
    axes[1].grid(True)

    # Plot 3: Imaginary part = sin(theta)
    axes[2].plot(theta, z.imag, "r-", linewidth=2)
    axes[2].set_xlabel("theta (radians)")
    axes[2].set_ylabel("Imaginary part")
    axes[2].set_title("sin(theta) = Imaginary part of e^(i*theta)")
    axes[2].grid(True)

    plt.tight_layout()
    plt.savefig("euler_formula.png", dpi=100)
    plt.show()

What does this code do?

This function generates 200 values of theta from 0 to 2*pi, then computes z = np.exp(1j * theta). Each value of z is a complex number sitting on the unit circle.

  • Panel 1 plots z.real vs z.imag to show the unit circle on the complex plane, with key angles marked as red dots.
  • Panel 2 plots z.real against theta, revealing that the real part of the circle IS cosine.
  • Panel 3 plots z.imag against theta, showing that the imaginary part IS sine.

What is 1j in Python?

Python uses 1j to represent the imaginary unit (what mathematicians write as i). So 1j * theta creates an imaginary number, and np.exp() knows how to raise e to a complex power using Euler's formula automatically.

  • Below the function definition, add these two lines to call the function:
print("\n--- Euler's Formula Visualization ---")
plot_euler_formula()
  • Save dft_analyzer.py.
  • Run the script in your terminal by running:
python dft_analyzer.py

After the sine wave plot and audio playback from Step 2, a new 3-panel figure appears. You should see the unit circle on the left, the cosine wave in the middle, and the sine wave on the right. The file euler_formula.png is saved in your project folder.

Plot window not appearing?

  • Make sure the plot_euler_formula() call is placed OUTSIDE and BELOW the function definition, not indented inside it.
  • If the plot saves but the window does not open, check euler_formula.png directly in your project folder. Some environments do not support interactive plot windows.
  • Close any open plot windows before running the script again. matplotlib waits for you to close each figure before showing the next one.

Still stuck? Help me debug why my euler_formula plot is not appearing when I run dft_analyzer.py.

Animate a rotating phasor to show how rotation encodes frequency

Euler's formula shows that a complex exponential traces a circle. A phasor is what happens when you add time to the equation. As time moves forward, the point on the circle spins, and the speed of that spin IS the frequency.

  • In dft_analyzer.py, add the following function below the plot_euler_formula() call:
def plot_rotating_phasor(freq):
    """Show how a phasor rotating at frequency f generates sine and cosine."""
    t = np.linspace(0, 1.0, 500, endpoint=False)
    phasor = np.exp(1j * 2 * np.pi * freq * t)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Left: phasor path on the complex plane
    axes[0].plot(phasor.real, phasor.imag, "b-", alpha=0.5)
    axes[0].plot(phasor.real[0], phasor.imag[0], "go", markersize=10, label="Start")
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_title(f"Phasor rotating at {freq} Hz")
    axes[0].set_xlabel("Real")
    axes[0].set_ylabel("Imaginary")
    axes[0].legend()

    # Right: real and imaginary parts over time
    axes[1].plot(t, phasor.real, "b-", label="Real = cos(2*pi*f*t)")
    axes[1].plot(t, phasor.imag, "r-", alpha=0.6, label="Imag = sin(2*pi*f*t)")
    axes[1].set_xlabel("Time (s)")
    axes[1].set_ylabel("Amplitude")
    axes[1].set_title(f"Components of e^(i*2*pi*{freq}*t)")
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("rotating_phasor.png", dpi=100)
    plt.show()

What does this code do?

This function computes phasor = np.exp(1j * 2 * np.pi * freq * t) for a 1-second time axis. That formula is just Euler's formula with a frequency parameter. The faster the frequency, the faster the point spins around the circle.

  • The left panel plots the phasor's path on the complex plane. At 3 Hz, it traces the unit circle 3 times in 1 second.
  • The right panel plots the real and imaginary parts over time. The real part is cosine and the imaginary part is sine, both oscillating 3 times per second.
  • Below the function, add these two lines to call it with a frequency of 3 Hz:
print("\n--- Rotating Phasor at 3 Hz ---")
plot_rotating_phasor(3)
  • Save dft_analyzer.py and run it again by running:
python dft_analyzer.py

After the previous plots close, a new 2-panel figure appears. The left panel shows a circle traced 3 times with a green "Start" dot. The right panel shows 3 complete cosine and sine oscillations over 1 second. The file rotating_phasor.png is saved in your project folder.

Only seeing 1 circle instead of 3?

  • Check that you passed 3 as the argument in plot_rotating_phasor(3). A value of 1 produces only a single rotation.
  • Make sure the function definition uses freq in the phasor formula, not a hardcoded number.

Still stuck? Help me debug why my rotating phasor plot shows the wrong number of rotations.

Connect phasors to the DFT

You now have all the visual intuition you need to understand how the DFT works. The connection is surprisingly direct.

How does the DFT use phasors to detect frequencies?

The DFT asks one question for each frequency bin k: "How much does my signal correlate with a phasor rotating at frequency k?" It answers this by multiplying the signal by e^(-i*2*pi*k*n/N) at every sample point, then summing all the products.

If the signal actually contains that frequency, the products all reinforce each other and the sum is large. If the signal does not contain that frequency, the products point in random directions around the circle and cancel out to near zero.

Every element of the DFT output X[k] is the result of multiplying the input signal by a specific rotating phasor and summing. You will implement this formula directly in the next step.

✔️ Awesome, I've got everything!

Great. Double check you have saved dft_analyzer.py and that both euler_formula.png and rotating_phasor.png are in your project folder.

ⓧ I'd like to double check the full code

import numpy as np
import matplotlib.pyplot as plt
import sounddevice as sd

# ============================================================
# STEP 2: Generate and Hear Sine Waves
# ============================================================

def generate_sine(freq, duration, sample_rate=44100):
    """Generate a sine wave at a given frequency.

    A sine wave is the simplest possible sound: a pure tone.
    The formula is: signal(t) = sin(2 * pi * freq * t)
    """
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * freq * t)
    return t, signal


def play_signal(signal, sample_rate=44100):
    """Play a numpy array as audio through speakers."""
    audio = signal / np.max(np.abs(signal))
    sd.play(audio.astype(np.float32), sample_rate)
    sd.wait()


# Generate a 440 Hz sine wave (the note A4)
sample_rate = 44100
duration = 1.0
t, signal_440 = generate_sine(440, duration, sample_rate)

# Plot the first 10 milliseconds to see the waveform
samples_10ms = int(sample_rate * 0.01)
plt.figure(figsize=(10, 4))
plt.plot(t[:samples_10ms] * 1000, signal_440[:samples_10ms])
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.title("440 Hz Sine Wave (A4 note) - first 10 ms")
plt.grid(True)
plt.tight_layout()
plt.savefig("sine_wave_440.png", dpi=100)
plt.show()

# Play the tones through your speakers
print("Playing 440 Hz tone (A4 note)...")
play_signal(signal_440, sample_rate)

_, signal_220 = generate_sine(220, duration, sample_rate)
print("Playing 220 Hz tone (A3 note - one octave lower)...")
play_signal(signal_220, sample_rate)


# ============================================================
# STEP 3: Complex Numbers and Euler's Formula
# ============================================================

def plot_euler_formula():
    """Visualize Euler's formula: e^(i*theta) = cos(theta) + i*sin(theta)."""
    theta = np.linspace(0, 2 * np.pi, 200)
    z = np.exp(1j * theta)

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # Plot 1: The unit circle on the complex plane
    axes[0].plot(z.real, z.imag, "b-", linewidth=2)
    angles = [0, np.pi / 4, np.pi / 2, np.pi, 3 * np.pi / 2]
    labels = ["0", "pi/4", "pi/2", "pi", "3pi/2"]
    for a, label in zip(angles, labels):
        point = np.exp(1j * a)
        axes[0].plot(point.real, point.imag, "ro", markersize=8)
        axes[0].annotate(
            f"theta={label}",
            (point.real, point.imag),
            textcoords="offset points",
            xytext=(10, 10),
            fontsize=8,
        )
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_xlabel("Real (cosine)")
    axes[0].set_ylabel("Imaginary (sine)")
    axes[0].set_title("Unit Circle: e^(i*theta)")

    # Plot 2: Real part = cos(theta)
    axes[1].plot(theta, z.real, "b-", linewidth=2)
    axes[1].set_xlabel("theta (radians)")
    axes[1].set_ylabel("Real part")
    axes[1].set_title("cos(theta) = Real part of e^(i*theta)")
    axes[1].grid(True)

    # Plot 3: Imaginary part = sin(theta)
    axes[2].plot(theta, z.imag, "r-", linewidth=2)
    axes[2].set_xlabel("theta (radians)")
    axes[2].set_ylabel("Imaginary part")
    axes[2].set_title("sin(theta) = Imaginary part of e^(i*theta)")
    axes[2].grid(True)

    plt.tight_layout()
    plt.savefig("euler_formula.png", dpi=100)
    plt.show()


def plot_rotating_phasor(freq):
    """Show how a phasor rotating at frequency f generates sine and cosine."""
    t = np.linspace(0, 1.0, 500, endpoint=False)
    phasor = np.exp(1j * 2 * np.pi * freq * t)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Left: phasor path on the complex plane
    axes[0].plot(phasor.real, phasor.imag, "b-", alpha=0.5)
    axes[0].plot(phasor.real[0], phasor.imag[0], "go", markersize=10, label="Start")
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_title(f"Phasor rotating at {freq} Hz")
    axes[0].set_xlabel("Real")
    axes[0].set_ylabel("Imaginary")
    axes[0].legend()

    # Right: real and imaginary parts over time
    axes[1].plot(t, phasor.real, "b-", label="Real = cos(2*pi*f*t)")
    axes[1].plot(t, phasor.imag, "r-", alpha=0.6, label="Imag = sin(2*pi*f*t)")
    axes[1].set_xlabel("Time (s)")
    axes[1].set_ylabel("Amplitude")
    axes[1].set_title(f"Components of e^(i*2*pi*{freq}*t)")
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("rotating_phasor.png", dpi=100)
    plt.show()


print("\n--- Euler's Formula Visualization ---")
plot_euler_formula()

print("\n--- Rotating Phasor at 3 Hz ---")
plot_rotating_phasor(3)

You now have a concrete, visual understanding of how complex exponentials encode frequency as rotation. Next up, you will translate this directly into code and implement the DFT formula from scratch.

Implement the DFT from Scratch

In the previous step, you explored Euler's formula and saw how rotating phasors on the complex plane generate sine and cosine. Now it's time to turn that insight into a working Discrete Fourier Transform.

The DFT multiplies your signal by phasors at every possible frequency and sums the results. If a frequency is present, the products reinforce. If it's absent, they cancel out. You'll implement this idea twice: first as a direct loop that maps line-for-line to the math formula, then as a vectorized matrix version that runs much faster.

In this step, get ready to:
  • Implement a loop-based DFT that directly translates the DFT formula.
  • Implement a matrix-based DFT using np.outer for speed.
  • Verify both implementations match numpy's FFT to machine precision.
Implement the DFT with explicit Python loops

The DFT formula says: for each frequency bin k, multiply every sample in the signal by a phasor rotating at that frequency, then sum the products. Two nested loops translate this formula directly into Python.

  • In dft_analyzer.py, scroll to the bottom of the file and find this line:
plot_rotating_phasor(3)
  • Add the following code below that line to create a section header and the loop-based DFT function:


# ============================================================
# STEP 4: Implement the DFT from Scratch
# ============================================================

def dft_loop(x):
    """Compute the DFT using explicit loops.

    Formula: X[k] = sum_{n=0}^{N-1} x[n] * e^(-i * 2*pi * k * n / N)

    This is O(N^2) - slow but maps directly to the math.
    """
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X

What does this code do?

  • The outer loop iterates over frequency bins k from 0 to N-1. Each bin represents one frequency the DFT checks for.
  • The inner loop multiplies every sample x[n] by the phasor e^(-2j*pi*k*n/N) at that frequency, then adds the result to X[k].
  • If frequency k is present in the signal, the phasor products line up and reinforce. If it's absent, they spin in circles and cancel to zero.
  • Save dft_analyzer.py.

You won't see output from this function yet. The verification code in the final substep will test it.

Why is this O(N^2)?

For each of the N frequency bins, the inner loop visits all N samples. That means N x N = N^2 total operations. For N=1000, that's a million multiplications. The loop version is intentionally slow because it maps directly to the math formula. The matrix version in the next substep computes the same result much faster.

Implement the DFT using matrix multiplication

The loop version computes one phasor product at a time. The matrix version computes all of them at once by building a matrix W where every entry is a precomputed phasor value. One matrix-vector multiplication then produces the entire DFT result.

  • In dft_analyzer.py, add the matrix-based DFT function directly below dft_loop:


def dft_matrix(x):
    """Compute the DFT using matrix multiplication.

    X = W @ x, where W[k,n] = e^(-i * 2*pi * k * n / N)

    Uses np.outer to build the exponent matrix. Same result, much faster.
    """
    N = len(x)
    n = np.arange(N)
    k = np.arange(N)
    W = np.exp(-2j * np.pi * np.outer(k, n) / N)
    return W @ x

What does this code do?

  • np.arange(N) creates arrays [0, 1, ..., N-1] for both the frequency indices k and the sample indices n.
  • np.outer(k, n) produces an NxN matrix where entry (i, j) equals k[i] * n[j]. This gives every k*n product needed for the DFT exponents in a single call.
  • np.exp(-2j * np.pi * ... / N) applies the exponential element-wise, building the full DFT matrix W.
  • W @ x performs one matrix-vector multiplication to compute all frequency bins simultaneously.
  • Save dft_analyzer.py.

You won't see output from this function yet. The next substep adds verification code to test both implementations.

Verify both implementations against numpy's FFT

Your two DFT functions should produce identical results. To prove it, you'll run both on test signals and compare their output against np.fft.fft(). Numpy's FFT is a battle-tested reference implementation.

  • Add the verification code below the dft_matrix function:


# Test with the impulse signal [1, 0, 0, 0]
test_signal = np.array([1.0, 0.0, 0.0, 0.0])

X_loop = dft_loop(test_signal)
X_matrix = dft_matrix(test_signal)
X_numpy = np.fft.fft(test_signal)

print("\n--- DFT Verification ---")
print("DFT of impulse [1, 0, 0, 0]:")
print(f"  Loop DFT:   {X_loop}")
print(f"  Matrix DFT: {X_matrix}")
print(f"  NumPy FFT:  {X_numpy}")
print(f"  Max error (loop vs numpy):   {np.max(np.abs(X_loop - X_numpy)):.2e}")
print(f"  Max error (matrix vs numpy): {np.max(np.abs(X_matrix - X_numpy)):.2e}")

# Test with a second signal
test_signal_2 = np.array([1.0, 2.0, 3.0, 4.0])
X_loop_2 = dft_loop(test_signal_2)
X_numpy_2 = np.fft.fft(test_signal_2)
print(f"\nDFT of [1, 2, 3, 4]:")
print(f"  Loop DFT:  {X_loop_2}")
print(f"  NumPy FFT: {X_numpy_2}")
print(f"  Max error: {np.max(np.abs(X_loop_2 - X_numpy_2)):.2e}")

What does this verification do?

  • The first test uses an impulse signal [1, 0, 0, 0]. An impulse has equal energy at all frequencies, so every DFT bin should be 1+0j.
  • The second test uses [1, 2, 3, 4] to confirm the math works on a non-trivial signal with varying sample values.
  • np.max(np.abs(...)) computes the largest difference between any two corresponding complex values. A result on the order of 1e-14 means the difference is at the limit of floating-point precision. Your DFT is mathematically correct.
  • Save dft_analyzer.py.
  • Run the script by entering this command in the VS Code terminal:
python dft_analyzer.py

Previous steps run first

The script runs everything from Steps 2 and 3 before reaching the DFT verification code. You'll see the sine wave plot, Euler's formula plot, and the phasor plot pop up. Close each plot window to continue. The audio tones will also play.

Once all the previous steps finish, look at the terminal for the DFT verification output at the bottom.

In the terminal, you should see output showing all three DFT values for the impulse signal are [1+0j, 1+0j, 1+0j, 1+0j]. The max error lines should show values on the order of 1e-14 or 0.00e+00. Both of your DFT implementations match numpy's FFT to machine precision.

Seeing errors or unexpected values?

  • If you see a NameError for dft_loop or dft_matrix, make sure the function definitions appear above the test code in the file. Python reads top to bottom.
  • If the max error is larger than 1e-10, double check the exponent in your functions. The sign must be negative (-2j) and the denominator must be N.
  • If you see an IndentationError, make sure your function bodies are indented with consistent spaces (4 spaces per level).

Why isn't my DFT matching numpy's FFT?

✔️ Awesome, I've got everything!

Great. Double check you've saved dft_analyzer.py before moving on.

ⓧ I'd like to double check the full code

Here is the complete dft_analyzer.py file up to this point. Compare it against yours and fix any differences.

import numpy as np
import matplotlib.pyplot as plt
import sounddevice as sd

# ============================================================
# STEP 2: Generate and Hear Sine Waves
# ============================================================

def generate_sine(freq, duration, sample_rate=44100):
    """Generate a sine wave at a given frequency.

    A sine wave is the simplest possible sound: a pure tone.
    The formula is: signal(t) = sin(2 * pi * freq * t)
    """
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * freq * t)
    return t, signal


def play_signal(signal, sample_rate=44100):
    """Play a numpy array as audio through speakers."""
    audio = signal / np.max(np.abs(signal))
    sd.play(audio.astype(np.float32), sample_rate)
    sd.wait()


# Generate a 440 Hz sine wave (the note A4)
sample_rate = 44100
duration = 1.0
t, signal_440 = generate_sine(440, duration, sample_rate)

# Plot the first 10 milliseconds to see the waveform
samples_10ms = int(sample_rate * 0.01)
plt.figure(figsize=(10, 4))
plt.plot(t[:samples_10ms] * 1000, signal_440[:samples_10ms])
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.title("440 Hz Sine Wave (A4 note) - first 10 ms")
plt.grid(True)
plt.tight_layout()
plt.savefig("sine_wave_440.png", dpi=100)
plt.show()

# Play the tones through your speakers
print("Playing 440 Hz tone (A4 note)...")
play_signal(signal_440, sample_rate)

_, signal_220 = generate_sine(220, duration, sample_rate)
print("Playing 220 Hz tone (A3 note - one octave lower)...")
play_signal(signal_220, sample_rate)


# ============================================================
# STEP 3: Complex Numbers and Euler's Formula
# ============================================================

def plot_euler_formula():
    """Visualize Euler's formula: e^(i*theta) = cos(theta) + i*sin(theta)."""
    theta = np.linspace(0, 2 * np.pi, 200)
    z = np.exp(1j * theta)

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # Plot 1: The unit circle on the complex plane
    axes[0].plot(z.real, z.imag, "b-", linewidth=2)
    angles = [0, np.pi / 4, np.pi / 2, np.pi, 3 * np.pi / 2]
    labels = ["0", "pi/4", "pi/2", "pi", "3pi/2"]
    for a, label in zip(angles, labels):
        point = np.exp(1j * a)
        axes[0].plot(point.real, point.imag, "ro", markersize=8)
        axes[0].annotate(
            f"theta={label}",
            (point.real, point.imag),
            textcoords="offset points",
            xytext=(10, 10),
            fontsize=8,
        )
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_xlabel("Real (cosine)")
    axes[0].set_ylabel("Imaginary (sine)")
    axes[0].set_title("Unit Circle: e^(i*theta)")

    # Plot 2: Real part = cos(theta)
    axes[1].plot(theta, z.real, "b-", linewidth=2)
    axes[1].set_xlabel("theta (radians)")
    axes[1].set_ylabel("Real part")
    axes[1].set_title("cos(theta) = Real part of e^(i*theta)")
    axes[1].grid(True)

    # Plot 3: Imaginary part = sin(theta)
    axes[2].plot(theta, z.imag, "r-", linewidth=2)
    axes[2].set_xlabel("theta (radians)")
    axes[2].set_ylabel("Imaginary part")
    axes[2].set_title("sin(theta) = Imaginary part of e^(i*theta)")
    axes[2].grid(True)

    plt.tight_layout()
    plt.savefig("euler_formula.png", dpi=100)
    plt.show()


def plot_rotating_phasor(freq):
    """Show how a phasor rotating at frequency f generates sine and cosine."""
    t = np.linspace(0, 1.0, 500, endpoint=False)
    phasor = np.exp(1j * 2 * np.pi * freq * t)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Left: phasor path on the complex plane
    axes[0].plot(phasor.real, phasor.imag, "b-", alpha=0.5)
    axes[0].plot(phasor.real[0], phasor.imag[0], "go", markersize=10, label="Start")
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_title(f"Phasor rotating at {freq} Hz")
    axes[0].set_xlabel("Real")
    axes[0].set_ylabel("Imaginary")
    axes[0].legend()

    # Right: real and imaginary parts over time
    axes[1].plot(t, phasor.real, "b-", label="Real = cos(2*pi*f*t)")
    axes[1].plot(t, phasor.imag, "r-", alpha=0.6, label="Imag = sin(2*pi*f*t)")
    axes[1].set_xlabel("Time (s)")
    axes[1].set_ylabel("Amplitude")
    axes[1].set_title(f"Components of e^(i*2*pi*{freq}*t)")
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("rotating_phasor.png", dpi=100)
    plt.show()


print("\n--- Euler's Formula Visualization ---")
plot_euler_formula()

print("\n--- Rotating Phasor at 3 Hz ---")
plot_rotating_phasor(3)


# ============================================================
# STEP 4: Implement the DFT from Scratch
# ============================================================

def dft_loop(x):
    """Compute the DFT using explicit loops.

    Formula: X[k] = sum_{n=0}^{N-1} x[n] * e^(-i * 2*pi * k * n / N)

    This is O(N^2) - slow but maps directly to the math.
    """
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X


def dft_matrix(x):
    """Compute the DFT using matrix multiplication.

    X = W @ x, where W[k,n] = e^(-i * 2*pi * k * n / N)

    Uses np.outer to build the exponent matrix. Same result, much faster.
    """
    N = len(x)
    n = np.arange(N)
    k = np.arange(N)
    W = np.exp(-2j * np.pi * np.outer(k, n) / N)
    return W @ x


# Test with the impulse signal [1, 0, 0, 0]
test_signal = np.array([1.0, 0.0, 0.0, 0.0])

X_loop = dft_loop(test_signal)
X_matrix = dft_matrix(test_signal)
X_numpy = np.fft.fft(test_signal)

print("\n--- DFT Verification ---")
print("DFT of impulse [1, 0, 0, 0]:")
print(f"  Loop DFT:   {X_loop}")
print(f"  Matrix DFT: {X_matrix}")
print(f"  NumPy FFT:  {X_numpy}")
print(f"  Max error (loop vs numpy):   {np.max(np.abs(X_loop - X_numpy)):.2e}")
print(f"  Max error (matrix vs numpy): {np.max(np.abs(X_matrix - X_numpy)):.2e}")

# Test with a second signal
test_signal_2 = np.array([1.0, 2.0, 3.0, 4.0])
X_loop_2 = dft_loop(test_signal_2)
X_numpy_2 = np.fft.fft(test_signal_2)
print(f"\nDFT of [1, 2, 3, 4]:")
print(f"  Loop DFT:  {X_loop_2}")
print(f"  NumPy FFT: {X_numpy_2}")
print(f"  Max error: {np.max(np.abs(X_loop_2 - X_numpy_2)):.2e}")

Your DFT is verified and mathematically correct. Next up, you'll use it to analyze a real multi-frequency signal and watch it decompose into individual frequency peaks.

Analyze Multi-Frequency Signals

Your DFT is verified and matches numpy's FFT to machine precision. Now it is time to put it to work on a signal that actually sounds like something.

Real-world audio is never a single pure tone. It is always a mix of frequencies layered on top of each other. In this step, you will create a composite signal from three sine waves, run your DFT on it, and see if it can pull the three original frequencies back out of the mix.

In this step, get ready to:
  • Create a composite signal by adding three sine waves of different frequencies and amplitudes.
  • Run your DFT and plot a magnitude spectrum showing the frequency content.
  • Detect and report the exact frequency peaks from the spectrum.
Create a composite signal from three frequencies

A single sine wave is easy to identify by ear. The interesting challenge is when multiple frequencies are mixed together into one signal. Your DFT should be able to separate them.

  • In dft_analyzer.py, scroll to the bottom of the file (after the Step 4 verification code) and add the Step 5 section comment and composite signal creation code:
# ============================================================
# STEP 5: Analyze Multi-Frequency Signals
# ============================================================

# Create a composite signal: 100 Hz + 250 Hz + 440 Hz
analysis_sr = 2000
analysis_dur = 0.5  # N = 1000 samples
t_analysis = np.linspace(
    0, analysis_dur, int(analysis_sr * analysis_dur), endpoint=False
)

composite = (
    1.0 * np.sin(2 * np.pi * 100 * t_analysis)
    + 0.5 * np.sin(2 * np.pi * 250 * t_analysis)
    + 0.7 * np.sin(2 * np.pi * 440 * t_analysis)
)

print("\n--- Frequency Analysis ---")
print("Composite signal: 1.0*sin(100Hz) + 0.5*sin(250Hz) + 0.7*sin(440Hz)")

Why 2000 Hz sample rate instead of 44100?

Your DFT builds an NxN matrix in memory. At 44100 Hz for 0.5 seconds, N would be 22050, creating a matrix of about 3.9 GB. At 2000 Hz for 0.5 seconds, N is 1000, and the matrix is only about 16 MB.

The frequency resolution is sample_rate / N = 2000 / 1000 = 2 Hz. Since 100, 250, and 440 are all exact multiples of 2, each frequency lands in a single bin with zero spectral leakage.

What does this code do?

  • analysis_sr = 2000 sets the sample rate low enough for your DFT matrix to fit comfortably in memory.
  • t_analysis creates 1000 evenly spaced time points across 0.5 seconds.
  • composite adds three sine waves together. Each has a different amplitude (1.0, 0.5, 0.7) so you can verify your DFT recovers them correctly.

The composite signal at 2000 Hz is too low a sample rate to hear properly through your speakers. To make it audible, you need to regenerate the same three-frequency mix at 44100 Hz.

  • Add the audible playback code directly below the print statements you just wrote:
# Play the composite signal at audible rate
t_audible = np.linspace(0, 1.0, 44100, endpoint=False)
composite_audible = (
    1.0 * np.sin(2 * np.pi * 100 * t_audible)
    + 0.5 * np.sin(2 * np.pi * 250 * t_audible)
    + 0.7 * np.sin(2 * np.pi * 440 * t_audible)
)
print("\nPlaying composite signal (100 + 250 + 440 Hz)...")
play_signal(composite_audible, 44100)

What does this code do?

This creates the same three-frequency signal at 44100 Hz so your speakers can play it. The play_signal function from Step 2 normalizes the signal and sends it to your audio output.

  • Save dft_analyzer.py.
  • Run the script in your terminal by running:
python dft_analyzer.py
  • After the earlier plots and DFT verification output, you should hear a complex chord-like tone. It sounds richer and messier than the pure 440 Hz sine wave from Step 2.

Audio not working?

If your system does not play audio, that is perfectly fine. Audio playback is optional. The core learning in this project comes from the code and visualizations.

Don't hear any sound?

  • Check that your system volume is turned up and not muted.
  • Make sure no other application is exclusively holding your audio output.
  • If the script crashes with an audio error, you can comment out the play_signal calls and continue. The frequency analysis works without audio.

Why can't I hear anything from my script?

Run your DFT and plot the magnitude spectrum

You have a composite signal. Now you need a function that feeds it through your DFT, converts the raw complex output into a magnitude spectrum, and plots the result.

  • In dft_analyzer.py, scroll up to the line that says # STEP 5: Analyze Multi-Frequency Signals.
  • Add the analyze_signal function directly below the section comment block, before the line # Create a composite signal: 100 Hz + 250 Hz + 440 Hz:
def analyze_signal(signal, sample_rate):
    """Run your DFT on a signal and plot the magnitude spectrum."""
    N = len(signal)

    # Run YOUR DFT (not numpy's)
    X = dft_matrix(signal)

    # Map DFT bins to frequencies in Hz
    freqs = np.fft.fftfreq(N, d=1.0 / sample_rate)

    # Only keep positive frequencies
    positive_mask = freqs >= 0
    freqs_pos = freqs[positive_mask]
    magnitudes = np.abs(X[positive_mask]) * 2 / N

    fig, axes = plt.subplots(2, 1, figsize=(10, 8))

    # Top: time domain (first 50 ms)
    samples_50ms = min(int(0.05 * sample_rate), N)
    t_plot = np.arange(samples_50ms) / sample_rate
    axes[0].plot(t_plot * 1000, signal[:samples_50ms])
    axes[0].set_xlabel("Time (ms)")
    axes[0].set_ylabel("Amplitude")
    axes[0].set_title("Time Domain (first 50 ms)")
    axes[0].grid(True)

    # Bottom: frequency domain
    axes[1].stem(freqs_pos, magnitudes, linefmt="b-", markerfmt="bo", basefmt="r-")
    axes[1].set_xlabel("Frequency (Hz)")
    axes[1].set_ylabel("Magnitude")
    axes[1].set_title("Magnitude Spectrum (computed by YOUR DFT)")
    axes[1].set_xlim(0, 600)
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("frequency_analysis.png", dpi=100)
    plt.show()

    return freqs_pos, magnitudes

What does this function do?

  • dft_matrix(signal) runs your hand-built DFT on the input signal. This is the same function you verified against numpy in Step 4.
  • np.fft.fftfreq(N, d=1.0 / sample_rate) maps each DFT bin index to a physical frequency in Hz. This is just a labeling helper, not a transform.
  • np.abs(X[positive_mask]) * 2 / N converts complex DFT output to magnitude. The 2 / N factor recovers the original signal amplitude from the DFT coefficients.
  • plt.stem draws a lollipop chart for the spectrum. Each stem shows the magnitude at one frequency bin.

Why multiply by 2/N?

The DFT spreads a real signal's energy across both positive and negative frequency bins (they are mirror images). Dividing by N undoes the DFT's summation scaling. Multiplying by 2 accounts for the energy in the negative-frequency mirror. The result is the original amplitude you put into each sine wave.

Now you need to call this function on your composite signal.

  • Scroll down to find the two print statements you added earlier. Find this line:
print("Composite signal: 1.0*sin(100Hz) + 0.5*sin(250Hz) + 0.7*sin(440Hz)")
  • Add the function call directly below that line, before the # Play the composite signal at audible rate comment:
freqs_result, mags_result = analyze_signal(composite, analysis_sr)
  • Save dft_analyzer.py and run the script again:
python dft_analyzer.py
  • A 2-panel figure should appear. The top panel shows the messy-looking composite waveform over the first 50 ms. The bottom panel shows the magnitude spectrum with three clear spikes at 100, 250, and 440 Hz.
  • Confirm that frequency_analysis.png has been saved in your dft-from-scratch folder.

Don't see three peaks?

  • Make sure analyze_signal is defined above the composite signal code. Python reads top-to-bottom, so the function must exist before it is called.
  • Check that analyze_signal calls dft_matrix(signal), not np.fft.fft(signal). You want to use your own DFT.
  • If you see peaks but the x-axis goes way beyond 600, make sure you have axes[1].set_xlim(0, 600) in your plotting code.

Help me debug my magnitude spectrum plot.

Detect and report the frequency peaks

The plot visually confirms three peaks. Now you will write code to automatically find them and print their exact frequencies and magnitudes.

  • In dft_analyzer.py, find the line:
freqs_result, mags_result = analyze_signal(composite, analysis_sr)
  • Add the peak detection code directly below that line, before the # Play the composite signal at audible rate comment:
# Detect peaks above threshold
peak_indices = np.where(mags_result > 0.1)[0]
print("\nDetected frequency peaks:")
for idx in peak_indices:
    print(f"  {freqs_result[idx]:.1f} Hz, magnitude: {mags_result[idx]:.3f}")

What does this code do?

  • np.where(mags_result > 0.1) finds every frequency bin whose magnitude exceeds the threshold of 0.1. This filters out the near-zero noise floor and keeps only the real signal peaks.
  • The loop prints each detected peak's frequency in Hz and its magnitude. These should match the amplitudes you put into the composite: 1.0, 0.5, and 0.7.
  • Save dft_analyzer.py and run the script:
python dft_analyzer.py
  • In your terminal output, look for the "Detected frequency peaks" section. You should see three lines: 100.0 Hz with magnitude 1.000, 250.0 Hz with magnitude 0.500, and 440.0 Hz with magnitude 0.700.

These magnitudes match the exact amplitudes you used when creating the composite signal. Your DFT has successfully decomposed a mixed signal back into its individual frequency components.

Seeing unexpected peak values?

  • If you see extra peaks or wrong magnitudes, check that your analyze_signal function uses np.abs(X[positive_mask]) * 2 / N for the magnitude calculation. Missing the * 2 / N factor will give raw DFT values instead of signal amplitudes.
  • Make sure analysis_sr is 2000 and analysis_dur is 0.5. Different values change N and the frequency resolution, which can shift where peaks land.

Why are my peak magnitudes wrong?

✔️ Awesome, I've got everything!

Great. Double check you have saved dft_analyzer.py.

ⓧ I'd like to double check the full code

Here is the complete dft_analyzer.py file including all code from Steps 2 through 5. Compare your file against this:

import numpy as np
import matplotlib.pyplot as plt
import sounddevice as sd

# ============================================================
# STEP 2: Generate and Hear Sine Waves
# ============================================================

def generate_sine(freq, duration, sample_rate=44100):
    """Generate a sine wave at a given frequency.

    A sine wave is the simplest possible sound: a pure tone.
    The formula is: signal(t) = sin(2 * pi * freq * t)
    """
    t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    signal = np.sin(2 * np.pi * freq * t)
    return t, signal


def play_signal(signal, sample_rate=44100):
    """Play a numpy array as audio through speakers."""
    audio = signal / np.max(np.abs(signal))
    sd.play(audio.astype(np.float32), sample_rate)
    sd.wait()


# Generate a 440 Hz sine wave (the note A4)
sample_rate = 44100
duration = 1.0
t, signal_440 = generate_sine(440, duration, sample_rate)

# Plot the first 10 milliseconds to see the waveform
samples_10ms = int(sample_rate * 0.01)
plt.figure(figsize=(10, 4))
plt.plot(t[:samples_10ms] * 1000, signal_440[:samples_10ms])
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.title("440 Hz Sine Wave (A4 note) - first 10 ms")
plt.grid(True)
plt.tight_layout()
plt.savefig("sine_wave_440.png", dpi=100)
plt.show()

# Play the tones through your speakers
print("Playing 440 Hz tone (A4 note)...")
play_signal(signal_440, sample_rate)

_, signal_220 = generate_sine(220, duration, sample_rate)
print("Playing 220 Hz tone (A3 note - one octave lower)...")
play_signal(signal_220, sample_rate)


# ============================================================
# STEP 3: Complex Numbers and Euler's Formula
# ============================================================

def plot_euler_formula():
    """Visualize Euler's formula: e^(i*theta) = cos(theta) + i*sin(theta)."""
    theta = np.linspace(0, 2 * np.pi, 200)
    z = np.exp(1j * theta)

    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # Plot 1: The unit circle on the complex plane
    axes[0].plot(z.real, z.imag, "b-", linewidth=2)
    angles = [0, np.pi / 4, np.pi / 2, np.pi, 3 * np.pi / 2]
    labels = ["0", "pi/4", "pi/2", "pi", "3pi/2"]
    for a, label in zip(angles, labels):
        point = np.exp(1j * a)
        axes[0].plot(point.real, point.imag, "ro", markersize=8)
        axes[0].annotate(
            f"theta={label}",
            (point.real, point.imag),
            textcoords="offset points",
            xytext=(10, 10),
            fontsize=8,
        )
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_xlabel("Real (cosine)")
    axes[0].set_ylabel("Imaginary (sine)")
    axes[0].set_title("Unit Circle: e^(i*theta)")

    # Plot 2: Real part = cos(theta)
    axes[1].plot(theta, z.real, "b-", linewidth=2)
    axes[1].set_xlabel("theta (radians)")
    axes[1].set_ylabel("Real part")
    axes[1].set_title("cos(theta) = Real part of e^(i*theta)")
    axes[1].grid(True)

    # Plot 3: Imaginary part = sin(theta)
    axes[2].plot(theta, z.imag, "r-", linewidth=2)
    axes[2].set_xlabel("theta (radians)")
    axes[2].set_ylabel("Imaginary part")
    axes[2].set_title("sin(theta) = Imaginary part of e^(i*theta)")
    axes[2].grid(True)

    plt.tight_layout()
    plt.savefig("euler_formula.png", dpi=100)
    plt.show()


def plot_rotating_phasor(freq):
    """Show how a phasor rotating at frequency f generates sine and cosine."""
    t = np.linspace(0, 1.0, 500, endpoint=False)
    phasor = np.exp(1j * 2 * np.pi * freq * t)

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Left: phasor path on the complex plane
    axes[0].plot(phasor.real, phasor.imag, "b-", alpha=0.5)
    axes[0].plot(phasor.real[0], phasor.imag[0], "go", markersize=10, label="Start")
    axes[0].set_xlim(-1.5, 1.5)
    axes[0].set_ylim(-1.5, 1.5)
    axes[0].set_aspect("equal")
    axes[0].grid(True)
    axes[0].set_title(f"Phasor rotating at {freq} Hz")
    axes[0].set_xlabel("Real")
    axes[0].set_ylabel("Imaginary")
    axes[0].legend()

    # Right: real and imaginary parts over time
    axes[1].plot(t, phasor.real, "b-", label="Real = cos(2*pi*f*t)")
    axes[1].plot(t, phasor.imag, "r-", alpha=0.6, label="Imag = sin(2*pi*f*t)")
    axes[1].set_xlabel("Time (s)")
    axes[1].set_ylabel("Amplitude")
    axes[1].set_title(f"Components of e^(i*2*pi*{freq}*t)")
    axes[1].legend()
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("rotating_phasor.png", dpi=100)
    plt.show()


print("\n--- Euler's Formula Visualization ---")
plot_euler_formula()

print("\n--- Rotating Phasor at 3 Hz ---")
plot_rotating_phasor(3)


# ============================================================
# STEP 4: Implement the DFT from Scratch
# ============================================================

def dft_loop(x):
    """Compute the DFT using explicit loops.

    Formula: X[k] = sum_{n=0}^{N-1} x[n] * e^(-i * 2*pi * k * n / N)

    This is O(N^2) - slow but maps directly to the math.
    """
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X


def dft_matrix(x):
    """Compute the DFT using matrix multiplication.

    X = W @ x, where W[k,n] = e^(-i * 2*pi * k * n / N)

    Uses np.outer to build the exponent matrix. Same result, much faster.
    """
    N = len(x)
    n = np.arange(N)
    k = np.arange(N)
    W = np.exp(-2j * np.pi * np.outer(k, n) / N)
    return W @ x


# Test with the impulse signal [1, 0, 0, 0]
test_signal = np.array([1.0, 0.0, 0.0, 0.0])

X_loop = dft_loop(test_signal)
X_matrix = dft_matrix(test_signal)
X_numpy = np.fft.fft(test_signal)

print("\n--- DFT Verification ---")
print("DFT of impulse [1, 0, 0, 0]:")
print(f"  Loop DFT:   {X_loop}")
print(f"  Matrix DFT: {X_matrix}")
print(f"  NumPy FFT:  {X_numpy}")
print(f"  Max error (loop vs numpy):   {np.max(np.abs(X_loop - X_numpy)):.2e}")
print(f"  Max error (matrix vs numpy): {np.max(np.abs(X_matrix - X_numpy)):.2e}")

# Test with a second signal
test_signal_2 = np.array([1.0, 2.0, 3.0, 4.0])
X_loop_2 = dft_loop(test_signal_2)
X_numpy_2 = np.fft.fft(test_signal_2)
print(f"\nDFT of [1, 2, 3, 4]:")
print(f"  Loop DFT:  {X_loop_2}")
print(f"  NumPy FFT: {X_numpy_2}")
print(f"  Max error: {np.max(np.abs(X_loop_2 - X_numpy_2)):.2e}")


# ============================================================
# STEP 5: Analyze Multi-Frequency Signals
# ============================================================

def analyze_signal(signal, sample_rate):
    """Run your DFT on a signal and plot the magnitude spectrum."""
    N = len(signal)

    # Run YOUR DFT (not numpy's)
    X = dft_matrix(signal)

    # Map DFT bins to frequencies in Hz
    freqs = np.fft.fftfreq(N, d=1.0 / sample_rate)

    # Only keep positive frequencies
    positive_mask = freqs >= 0
    freqs_pos = freqs[positive_mask]
    magnitudes = np.abs(X[positive_mask]) * 2 / N

    fig, axes = plt.subplots(2, 1, figsize=(10, 8))

    # Top: time domain (first 50 ms)
    samples_50ms = min(int(0.05 * sample_rate), N)
    t_plot = np.arange(samples_50ms) / sample_rate
    axes[0].plot(t_plot * 1000, signal[:samples_50ms])
    axes[0].set_xlabel("Time (ms)")
    axes[0].set_ylabel("Amplitude")
    axes[0].set_title("Time Domain (first 50 ms)")
    axes[0].grid(True)

    # Bottom: frequency domain
    axes[1].stem(freqs_pos, magnitudes, linefmt="b-", markerfmt="bo", basefmt="r-")
    axes[1].set_xlabel("Frequency (Hz)")
    axes[1].set_ylabel("Magnitude")
    axes[1].set_title("Magnitude Spectrum (computed by YOUR DFT)")
    axes[1].set_xlim(0, 600)
    axes[1].grid(True)

    plt.tight_layout()
    plt.savefig("frequency_analysis.png", dpi=100)
    plt.show()

    return freqs_pos, magnitudes


# Create a composite signal: 100 Hz + 250 Hz + 440 Hz
analysis_sr = 2000
analysis_dur = 0.5  # N = 1000 samples
t_analysis = np.linspace(
    0, analysis_dur, int(analysis_sr * analysis_dur), endpoint=False
)

composite = (
    1.0 * np.sin(2 * np.pi * 100 * t_analysis)
    + 0.5 * np.sin(2 * np.pi * 250 * t_analysis)
    + 0.7 * np.sin(2 * np.pi * 440 * t_analysis)
)

print("\n--- Frequency Analysis ---")
print("Composite signal: 1.0*sin(100Hz) + 0.5*sin(250Hz) + 0.7*sin(440Hz)")

freqs_result, mags_result = analyze_signal(composite, analysis_sr)

# Detect peaks above threshold
peak_indices = np.where(mags_result > 0.1)[0]
print("\nDetected frequency peaks:")
for idx in peak_indices:
    print(f"  {freqs_result[idx]:.1f} Hz, magnitude: {mags_result[idx]:.3f}")

# Play the composite signal at audible rate
t_audible = np.linspace(0, 1.0, 44100, endpoint=False)
composite_audible = (
    1.0 * np.sin(2 * np.pi * 100 * t_audible)
    + 0.5 * np.sin(2 * np.pi * 250 * t_audible)
    + 0.7 * np.sin(2 * np.pi * 440 * t_audible)
)
print("\nPlaying composite signal (100 + 250 + 440 Hz)...")
play_signal(composite_audible, 44100)

Secret mission

Implement the Inverse DFT

Your DFT can decompose a signal into frequencies. But can you reverse the process? Implement the Inverse DFT to perfectly reconstruct the original signal from its frequency components, proving that no information is lost in the transform.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your files for future use or delete them entirely.

Resources you used:

  • The dft-from-scratch folder on your Desktop (contains dft_analyzer.py and all plot images).
  • Python packages installed via pip: numpy, matplotlib, sounddevice.

Keep everything running

No action needed. Choose this if you plan to continue building on this project.

  • Your dft_analyzer.py file and DFT functions are the foundation for Project 2, where you will add windowing, framing, and build a full spectrogram.
  • The plot images (sine_wave_440.png, euler_formula.png, rotating_phasor.png, frequency_analysis.png, idft_reconstruction.png) are useful references for reviewing Euler's formula and the DFT.
  • The installed packages (numpy, matplotlib, sounddevice) will remain available for other Python projects.

Pause - I'll come back to this later

Since this project is entirely local, there are no running processes to pause. Your files will stay on disk until you return.

  • Close the dft-from-scratch folder in VS Code.
  • When you are ready to continue, reopen the folder in VS Code and pick up where you left off.

Delete - I don't want to use this again

Remove all project files and installed packages. There are no cloud resources or persistent state to clean up.

  • Close the dft-from-scratch folder in VS Code.
  • Delete the project folder by running this command in your terminal.

macOS

rm -rf ~/Desktop/dft-from-scratch

Windows

rmdir /s /q %USERPROFILE%\Desktop\dft-from-scratch
  • Uninstall the Python packages by running:
pip uninstall numpy matplotlib sounddevice -y

Keep packages for other projects

If you use numpy or matplotlib in other Python projects, skip the uninstall step. Only uninstall packages you no longer need.

That's a Wrap!

That's a Wrap!

Nice work! You just built a frequency analyzer from the ground up, starting with pure math and ending with a tool that pulls apart a complex audio signal into its individual frequencies.

You've learned how to:

  • Implemented the Discrete Fourier Transform from scratch in two forms: an explicit loop version mapping directly to the math formula, and an optimized matrix version using np.outer. Both matched numpy's FFT to machine precision.
  • Visualized Euler's formula and rotating phasors on the complex plane, connecting complex exponentials to sine and cosine to understand how the DFT detects frequencies.
  • Analyzed a composite 3-frequency signal and recovered the exact frequencies and amplitudes from the magnitude spectrum, proving your DFT can decompose real signals into their individual components.
  • Secret Mission: Implemented the Inverse DFT and proved perfect round-trip signal reconstruction by running DFT followed by IDFT to recover the original signal.

Ready to quiz yourself?