Code Lab: Nyquist Theorem and Sampling Demonstration

The Nyquist-Shannon sampling theorem is one of the most important results in digital audio. It states that a continuous signal can be perfectly reconstructed from its samples if the sampling rate is at least twice the highest frequency present in the signal. When this condition is violated, aliasing occurs: high-frequency components masquerade as lower frequencies, introducing artifacts that cannot be removed after the fact. In this lab, you will build intuition for the theorem by sampling signals at different rates, visualizing aliasing, and reconstructing a signal from discrete samples using sinc interpolation.

Setting Up the Signal

We begin by creating a "continuous" reference signal (approximated at a very high sample rate) and then sampling it at three different rates: well above Nyquist, just above Nyquist, and below Nyquist.

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

# Parameters
SIGNAL_FREQ = 800.0       # Hz -- the frequency we want to capture
CONTINUOUS_SR = 96000      # Hz -- approximates a continuous waveform
DURATION = 0.006           # seconds (enough for ~5 cycles at 800 Hz)

# Three sampling rates to compare
FS_HIGH = 8000             # Well above Nyquist (Nyquist = 4000 Hz)
FS_MARGINAL = 1800         # Just above Nyquist (Nyquist = 900 Hz)
FS_LOW = 1200              # Below Nyquist (Nyquist = 600 Hz) -- will alias

# Build the "continuous" signal
t_cont = np.linspace(0, DURATION, int(CONTINUOUS_SR * DURATION), endpoint=False)
x_cont = np.sin(2 * np.pi * SIGNAL_FREQ * t_cont)

The variable x_cont now holds a densely sampled sine wave that serves as our ground truth. The three sampling rates let us observe what happens when the Nyquist condition is satisfied comfortably, barely satisfied, and violated.

Sampling at Different Rates

def sample_at_rate(t_cont, x_cont, fs):
    """Sample the continuous signal at rate fs using linear interpolation."""
    t_samples = np.arange(0, t_cont[-1], 1.0 / fs)
    interpolator = interp1d(t_cont, x_cont, kind='linear')
    x_samples = interpolator(t_samples)
    return t_samples, x_samples

t_high, x_high = sample_at_rate(t_cont, x_cont, FS_HIGH)
t_marg, x_marg = sample_at_rate(t_cont, x_cont, FS_MARGINAL)
t_low, x_low   = sample_at_rate(t_cont, x_cont, FS_LOW)

Visualizing Aliasing

When we sample at FS_LOW = 1200 Hz, the Nyquist frequency is only 600 Hz, which is below our 800 Hz signal. The sampled points trace out a wave at the alias frequency rather than the true frequency. The alias frequency is given by f_alias = |f_signal - n * f_sample| for the integer n that minimizes the result.

# Compute alias frequency
n = round(SIGNAL_FREQ / FS_LOW)
f_alias = abs(SIGNAL_FREQ - n * FS_LOW)
print(f"True frequency: {SIGNAL_FREQ} Hz")
print(f"Sampling rate:  {FS_LOW} Hz (Nyquist = {FS_LOW / 2} Hz)")
print(f"Alias frequency: {f_alias} Hz  <-- this is what we hear instead!")

fig, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True)
fig.suptitle(f"Sampling a {int(SIGNAL_FREQ)} Hz Sine Wave at Three Rates",
             fontsize=13, fontweight='bold')

configs = [
    (t_high, x_high, FS_HIGH, "Well above Nyquist", "darkgreen"),
    (t_marg, x_marg, FS_MARGINAL, "Just above Nyquist", "darkorange"),
    (t_low,  x_low,  FS_LOW,  f"Below Nyquist -- aliases to {f_alias:.0f} Hz", "firebrick"),
]

for ax, (ts, xs, fs, label, color) in zip(axes, configs):
    ax.plot(t_cont * 1000, x_cont, color='lightblue', linewidth=1.0,
            alpha=0.7, label="True signal")
    ax.stem(ts * 1000, xs, linefmt=color, markerfmt='o', basefmt='gray',
            label=f"Samples (fs={fs} Hz)")
    ax.set_ylabel("Amplitude")
    ax.set_title(f"{label}  |  fs = {fs} Hz, Nyquist = {fs // 2} Hz", fontsize=9)
    ax.legend(fontsize=8, loc='upper right')
    ax.grid(True, alpha=0.3)

axes[-1].set_xlabel("Time (ms)")
plt.tight_layout()
plt.show()

In the bottom panel you can see that the sample points no longer follow the true sine wave. Instead, they trace a slower oscillation at the alias frequency. This artifact is permanent -- once the samples are taken, no amount of processing can recover the original 800 Hz tone.

Reconstructing a Signal with Sinc Interpolation

The Nyquist theorem guarantees that when the sampling condition is met, the original signal can be perfectly reconstructed using sinc interpolation. Each sample contributes a shifted and scaled sinc function, and the sum of all these sinc functions recovers the continuous waveform exactly.

def sinc_reconstruct(t_samples, x_samples, t_recon, fs):
    """
    Reconstruct a continuous signal from samples using the Whittaker-Shannon
    (sinc) interpolation formula:
        x(t) = sum_n x[n] * sinc(fs * (t - n/fs))
    """
    x_recon = np.zeros_like(t_recon)
    T = 1.0 / fs
    for n, x_n in enumerate(x_samples):
        x_recon += x_n * np.sinc(fs * (t_recon - n * T))
    return x_recon

# Reconstruct from the well-sampled case
t_recon = np.linspace(0, DURATION, 2000, endpoint=False)
x_recon = sinc_reconstruct(t_high, x_high, t_recon, FS_HIGH)

# Measure reconstruction error
x_true_at_recon = np.sin(2 * np.pi * SIGNAL_FREQ * t_recon)
max_error = np.max(np.abs(x_recon - x_true_at_recon))
print(f"Maximum reconstruction error (sinc): {max_error:.6f}")

plt.figure(figsize=(10, 4))
plt.plot(t_recon * 1000, x_true_at_recon, 'b-', linewidth=1.5, label="True signal")
plt.plot(t_recon * 1000, x_recon, 'r--', linewidth=1.5, label="Sinc reconstruction")
plt.stem(t_high * 1000, x_high, linefmt='green', markerfmt='go',
         basefmt='gray', label=f"Samples (fs={FS_HIGH} Hz)")
plt.xlabel("Time (ms)")
plt.ylabel("Amplitude")
plt.title("Sinc Interpolation: Perfect Reconstruction from Samples")
plt.legend(fontsize=9)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

The red dashed reconstruction line lies directly on top of the blue true signal. This is the Nyquist theorem in action: given enough samples per cycle, the original waveform is fully recoverable.

Try It Yourself

Explore the boundary. Change FS_MARGINAL to values very close to 2 * SIGNAL_FREQ (e.g., 1601, 1610, 1650 Hz). At what sampling rate does the reconstruction start to visibly degrade? How does the number of samples per cycle relate to the visual quality?
Multi-frequency aliasing. Replace the single sine wave with a composite signal containing two frequencies (e.g., 800 Hz and 3500 Hz). Sample at 4000 Hz. Which component aliases, and to what frequency? Verify your prediction by plotting the samples alongside both the original and the alias waveform.
Compare interpolation methods. Modify the reconstruction step to use linear interpolation (interp1d with kind='linear') instead of sinc interpolation. Plot both reconstructions on the same axes. Where does linear interpolation introduce the most error, and why? How does increasing the sampling rate reduce the difference between the two methods?