Tutorial 2: Fourier Analysis#
Week 2, Day 2: Time series and signal processing
By Neuromatch Academy (community contribution)
Content Creators: Sricharan Sunder, Selva Maran, Ilia Leonov
Content Reviewers: Andreas Neef, Richard Gao
Production editors: Konstantine Tsafatinos
Tutorial Objectives#
Estimated timing of tutorial: 45 min
This tutorial builds the core intuition behind Fourier analysis from the ground up.
In the previous tutorial, we learned about sine waves and how their 3 parameters (amplitude, frequency, phase) shape them.
In this tutorial, we will use sines (and cosines) as the elementary building block of every signal: we will see how to decompose a signal into a weighted sum of sine waves and how a handful of components can combine to reproduce a target. This can be interpreted as a change of basis, similar to what you already learned in the Dimensionality Reduction tutorial.
We will use this idea to build the Fourier basis: decomposing periodic signals (like a square or sawtooth wave) as sums of orthogonal sine and cosine functions, watching the approximation converge as we add more terms, and even building an aperiodic or non-rhythmic Gaussian pulse the same way — where the phases decides whether the result is a sharp spike or noise.
Finally, we will introduce and implement the (Discrete) Fourier Transform (DFT), which is how we compute the Fourier coefficients for any signal.
The fundamental concept of Fourier analysis is that any signal can be represented as a sum of sine waves. This is the starting point for many signal processing techniques used across neuroscience, including filtering and spectral analysis, as we will see in later tutorials.
By the end of this tutorial, you will be able to:
Reconstruct a target signal from its Fourier components, i.e. a weighted sum of sine waves
Read a magnitude / power spectrum and relate it to the time-domain components, and visualize a complex exponential
See how an aperiodic pulse is built from sine waves, and how phase alone separates a sharp spike from noise
Write down the Discrete Fourier Transform (DFT), as well as how to determine the number of components and their coefficients for any signal
Setup#
Install and import feedback gadget#
Show code cell source
# @title Install and import feedback gadget
!pip3 install vibecheck datatops --quiet
!pip3 install remfile h5py --quiet
from vibecheck import DatatopsContentReviewContainer
def content_review(notebook_section: str):
return DatatopsContentReviewContainer(
"",
notebook_section,
{
"url": "https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab",
"name": "neuromatch_cn",
"user_key": "y1x3mpx5",
},
).render()
feedback_prefix = "W2D2_T2"
Imports#
Show code cell source
# @title Imports
# General
import io
import time
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import ipywidgets as widgets
import warnings
warnings.filterwarnings('ignore')
# Sine waves
from scipy.integrate import quad
# Spectral Analysis
import remfile, h5py
from scipy.fft import fft, fftfreq
Figure Settings#
Show code cell source
# @title Figure Settings
import logging
logging.getLogger('matplotlib.font_manager').disabled = True
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
plt.style.use("https://raw.githubusercontent.com/NeuromatchAcademy/course-content/main/nma.mplstyle")
Section 1: Signal decomposition#
Towards the Fourier Transform#
You previously learned about the rules of sampling using the example of individual sine waves. This might appear to be an arbitrary and specific choice of signals. However, it can be trivially extended to any sum of sine waves, where the rules of sampling (and the Nyquist theorem) still apply.
As a warm-up exercise, here we create a “mystery signal” and then attempt to find the correct way to approximate it as a sum of two sine waves.
Video 1: Finding the right sine waves#
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Video_1_Fourier_Intro")
Interactive Demo 1: Decomposing the mystery signal#
A mystery signal has been constructed from a constant offset plus two sine waves:
You know the frequencies — \(f_1\) = 3 Hz and \(f_2\) = 7 Hz — but the weights \(w_0\), \(w_1\) and \(w_2\), and phases \(\phi_1\) and \(\phi_2\) are hidden.
Your goal: use the sliders in the interactive demo to adjust the weights until your reconstructed signal (blue) matches the mystery signal (gray) as closely as possible. The quality of your match is measured by the Mean Squared Error (MSE) shown on the plot — lower is better. Try to drive it as close to 0 as you can!
The demo below shows multiple “views” of the same signal:
Top — Reconstruction (the sum): your reconstructed signal (blue) overlaid on the ground truth mystery signal (gray dashed). This is the sum of the three weighted components. The MSE badge in the corner reports how close your sum is to the true signal.
Middle — Components: the three weighted pieces drawn separately before they are added. The dotted cyan line is the DC offset \(w_0\) (a flat line at height \(w_0\)), and the orange and purple curves are \(w_1\sin(2\pi f_1 t + \phi_1)\) and \(w_2\sin(2\pi f_2 t + \phi_2)\). Adding these three curves point-by-point gives you the blue curve above.
Bottom row:
Left: coefficients draw in the complex domain, where the radius of the circle is the amplitude / weight of that frequency, and the angle where the dot currently sits is the phase. We’ll see why this is useful in a bit.
Middle: magnitude or amplitude spectrum, where the height of each bar is the amplitude of that frequency. The 0 Hz bar is the DC weight, and the bars at \(f_1\) and \(f_2\) are the amplitudes of the two sines, \(w_1\) and \(w_2\).
Right: phase spectrum, where the y-axis is between -π and π.
These views provide the same information, indexing either by time or frequency — the bottom row is the picture a Fourier analysis would hand you.
As you move the sliders, watch all of them update together: the spectrum bars change height, the component sines rescale, and the sum on top morphs toward the ground truth signal.
A note on the weights \(w\):
The weights here are amplitudes — they represent how much of each frequency component is present in the signal. Each weight independently controls the contribution of its frequency. \(w_0\) is the DC offset — the 0 Hz sinusoidal component, i.e. the constant average level the whole signal sits on. It is a Fourier coefficient just like the others, and you tune it with its own slider. It also has a phase, but note what happens when you change it?
Decompose the signal#
Show code cell source
# @title Decompose the signal
def make_signal(t, w0, w1, w2, f1, f2, phi1=0.0, phi2=0.0):
"""Two-tone signal with a DC offset.
S(t) = w0 + w1 · sin(2π f1 t + φ1) + w2 · sin(2π f2 t + φ2)
Args:
t (ndarray): time samples in seconds, shape (n_samples,)
w0 (float): DC offset (0 Hz component)
w1, w2 (float): amplitudes of the two sine components
f1, f2 (float): frequencies in Hz
phi1, phi2 (float): phase offsets in radians (default 0)
Returns:
ndarray: S(t), shape (n_samples,)
"""
# Compute S(t) = w0 + w1·sin(2π f1 t + φ1) + w2·sin(2π f2 t + φ2)
s = (w0
+ w1 * np.sin(2 * np.pi * f1 * t + phi1)
+ w2 * np.sin(2 * np.pi * f2 * t + phi2))
return s
def build_decomposition_demo():
# ── Instructor config ──
F1, F2 = 3, 7
W0_TRUE, W1_TRUE, W2_TRUE = 0.8, 0.6, 1.4
PHI0_TRUE = 0.0
PHI1_TRUE, PHI2_TRUE = 0.25*np.pi, 0.15*np.pi
DURATION, FS = 1.0, 1000
t = np.linspace(0, DURATION, int(FS * DURATION), endpoint=False)
y_true = (W0_TRUE
+ W1_TRUE * np.sin(2 * np.pi * F1 * t + PHI1_TRUE)
+ W2_TRUE * np.sin(2 * np.pi * F2 * t + PHI2_TRUE))
rng = np.random.default_rng()
PHI1_RAND, PHI2_RAND = rng.uniform(-np.pi, np.pi, 2)
# Harmonic structure for the spectrum panel (0 Hz prepended for the DC term)
_F0_spec = math.gcd(int(F1), int(F2))
_all_h = [0] + list(range(1, int(max(F1, F2)) // _F0_spec + 1))
_h1 = int(F1) // _F0_spec
_h2 = int(F2) // _F0_spec
# Current displayed phases (used for random-phase animation)
_phases = [0.0, 0.0]
img_out = widgets.Image(format='png')
def draw(phases):
phi1, phi2 = phases
w0 = w0_slider.value
w1 = w1_slider.value
w2 = w2_slider.value
y_student = make_signal(t, w0, w1, w2, F1, F2, phi1, phi2)
mse = np.mean((y_true - y_student) ** 2)
# Phase-independent y-limits: use DC offset + sum of |coefficients| as upper
# bound so the axis only moves when magnitudes change, not during animation
y_abs_max = max(abs(W0_TRUE) + abs(W1_TRUE) + abs(W2_TRUE),
abs(w0) + abs(w1) + abs(w2), 0.5) * 1.2
ylim = (-y_abs_max, y_abs_max)
spec_h, spec_c = [], []
for h in _all_h:
if h == 0:
spec_h.append(abs(w0)); spec_c.append('#0891b2')
elif h == _h1:
spec_h.append(abs(w1)); spec_c.append('#e85d04')
elif h == _h2:
spec_h.append(abs(w2)); spec_c.append('#7c3aed')
else:
spec_h.append(0.0); spec_c.append('#dddddd')
spec_labels = [f'{h * _F0_spec} Hz' for h in _all_h]
fig = plt.figure(figsize=(14, 9))
gs = fig.add_gridspec(3, 3, height_ratios=[1, 1, 1.5])
ax_main = fig.add_subplot(gs[0, :])
ax_comp = fig.add_subplot(gs[1, :])
ax_complex = fig.add_subplot(gs[2, 0])
ax_mag_spec = fig.add_subplot(gs[2, 1])
ax_ph_spec = fig.add_subplot(gs[2, 2])
# ── Top: reconstruction ───────────────────────────────
ax_main.plot(t, y_true, color='#888888', linestyle='--', linewidth=2,
label='Mystery signal', alpha=0.85)
ax_main.plot(t, y_student, color='#2563EB', linewidth=2,
label='Your signal')
ax_main.set_title('Reconstruction (sum of components)', fontsize=13)
ax_main.set_xlabel('Time (s)')
ax_main.set_ylabel('Amplitude')
ax_main.legend(loc='upper right')
ax_main.set_xlim(0, DURATION)
ax_main.set_ylim(ylim)
ax_main.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
mse_color = '#16a34a' if mse < 0.01 else '#dc2626'
ax_main.text(0.02, 0.95, f'MSE: {mse:.4f}',
transform=ax_main.transAxes, fontsize=12, verticalalignment='top',
color=mse_color,
bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
edgecolor=mse_color, alpha=0.8))
phi0 = phi0_slider.value * np.pi
# ── Middle: components overlaid ──────────────────────
phi1_pi = phi1 / np.pi
phi2_pi = phi2 / np.pi
phi1_txt = f' + {phi1_pi:.2f}π' if abs(phi1) > 1e-3 else ''
phi2_txt = f' + {phi2_pi:.2f}π' if abs(phi2) > 1e-3 else ''
c1 = w1 * np.sin(2 * np.pi * F1 * t + phi1)
c2 = w2 * np.sin(2 * np.pi * F2 * t + phi2)
ax_comp.plot(t, np.full_like(t, w0), color='#0891b2', linewidth=2,
linestyle=':', label=f'w₀={w0:.2f} (DC)')
ax_comp.plot(t, c1, color='#e85d04', linewidth=2,
label=f'w₁={w1:.2f}·sin(2π·{F1}Hz·t{phi1_txt})')
ax_comp.plot(t, c2, color='#7c3aed', linewidth=2,
label=f'w₂={w2:.2f}·sin(2π·{F2}Hz·t{phi2_txt})')
ax_comp.set_title('Components', fontsize=12)
ax_comp.set_xlabel('Time (s)')
ax_comp.set_ylabel('Amplitude')
ax_comp.legend(loc='upper right', fontsize=9)
ax_comp.set_xlim(0, DURATION)
ax_comp.set_ylim(ylim)
ax_comp.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
# ── Bottom-left: complex-domain coefficient view ─────
rmax = max(abs(w0), abs(w1), abs(w2), 0.5) * 1.25
theta = np.linspace(-np.pi, np.pi, 300)
# Circles (radius = coefficient magnitude)
ax_complex.plot(abs(w0) * np.cos(theta), abs(w0) * np.sin(theta),
color='#0891b2', linewidth=1.8, alpha=0.9)
ax_complex.plot(abs(w1) * np.cos(theta), abs(w1) * np.sin(theta),
color='#e85d04', linewidth=1.8, alpha=0.9)
ax_complex.plot(abs(w2) * np.cos(theta), abs(w2) * np.sin(theta),
color='#7c3aed', linewidth=1.8, alpha=0.9)
# Phase dots on each circle
p0x, p0y = abs(w0) * np.cos(phi0), abs(w0) * np.sin(phi0)
p1x, p1y = abs(w1) * np.cos(phi1), abs(w1) * np.sin(phi1)
p2x, p2y = abs(w2) * np.cos(phi2), abs(w2) * np.sin(phi2)
ax_complex.scatter([p0x], [p0y], color='#0891b2', s=55, zorder=3)
ax_complex.scatter([p1x], [p1y], color='#e85d04', s=55, zorder=3)
ax_complex.scatter([p2x], [p2y], color='#7c3aed', s=55, zorder=3)
if show_answer_checkbox.value:
# Show ground-truth coefficients only in the complex-domain panel
a0x, a0y = abs(W0_TRUE) * np.cos(PHI0_TRUE), abs(W0_TRUE) * np.sin(PHI0_TRUE)
a1x, a1y = abs(W1_TRUE) * np.cos(PHI1_TRUE), abs(W1_TRUE) * np.sin(PHI1_TRUE)
a2x, a2y = abs(W2_TRUE) * np.cos(PHI2_TRUE), abs(W2_TRUE) * np.sin(PHI2_TRUE)
ax_complex.scatter([a0x], [a0y], marker='*', s=220, color='#0891b2',
edgecolors='black', linewidths=0.8, zorder=5)
ax_complex.scatter([a1x], [a1y], marker='*', s=220, color='#e85d04',
edgecolors='black', linewidths=0.8, zorder=5)
ax_complex.scatter([a2x], [a2y], marker='*', s=220, color='#7c3aed',
edgecolors='black', linewidths=0.8, zorder=5)
# Radius guides
ax_complex.plot([0, p0x], [0, p0y], color='#0891b2', linewidth=1.2, alpha=0.7)
ax_complex.plot([0, p1x], [0, p1y], color='#e85d04', linewidth=1.2, alpha=0.7)
ax_complex.plot([0, p2x], [0, p2y], color='#7c3aed', linewidth=1.2, alpha=0.7)
ax_complex.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_complex.axvline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_complex.set_xlim(-rmax, rmax)
ax_complex.set_ylim(-rmax, rmax)
ax_complex.set_aspect('equal', adjustable='box')
ax_complex.set_title('Complex Coefficients', fontsize=12)
ax_complex.set_xlabel('Real')
ax_complex.set_ylabel('Imag')
ax_complex.grid(alpha=0.2)
# ── Bottom-middle: magnitude spectrum ─────────────────
ax_mag_spec.bar(_all_h, spec_h, color=spec_c, width=0.55, alpha=0.9)
ax_mag_spec.set_xticks(_all_h)
ax_mag_spec.set_xticklabels(spec_labels, fontsize=9)
ax_mag_spec.set_title('Amplitude Spectrum', fontsize=12)
ax_mag_spec.set_ylabel('Magnitude (w)')
ax_mag_spec.set_ylim(0, max(3.2, max(spec_h) * 1.3))
ax_mag_spec.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
# ── Bottom-right: phase spectrum ─────────────────────
phase_x = [0, _h1, _h2]
phase_vals = [phi0, phi1, phi2]
phase_colors = ['#0891b2', '#e85d04', '#7c3aed']
ax_ph_spec.vlines(phase_x, 0, phase_vals, colors=phase_colors, linewidth=1.6, alpha=0.8)
ax_ph_spec.scatter(phase_x, phase_vals, color=phase_colors, s=45, zorder=3)
ax_ph_spec.set_xticks(_all_h)
ax_ph_spec.set_xticklabels(spec_labels, fontsize=9)
ax_ph_spec.set_title('Phase Spectrum', fontsize=12)
ax_ph_spec.set_ylabel('Phase (rad)')
ax_ph_spec.set_ylim(-np.pi, np.pi)
ax_ph_spec.set_yticks([-np.pi, -np.pi/2, 0, np.pi/2, np.pi])
ax_ph_spec.set_yticklabels(['-π', '-π/2', '0', 'π/2', 'π'])
ax_ph_spec.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
fig.tight_layout()
buf = io.BytesIO()
fig.savefig(buf, format='png', dpi=100, bbox_inches='tight')
plt.close(fig)
buf.seek(0)
img_out.value = buf.read()
def update(_change=None):
if phase_radio.value.startswith('Random'):
draw(_phases)
else:
draw([phi1_slider.value * np.pi, phi2_slider.value * np.pi])
def on_phase_change(change):
is_random = change['new'].startswith('Random')
phi1_slider.disabled = is_random
phi2_slider.disabled = is_random
if is_random:
target = [PHI1_RAND, PHI2_RAND]
start = [phi1_slider.value * np.pi, phi2_slider.value * np.pi]
else:
target = [phi1_slider.value * np.pi, phi2_slider.value * np.pi]
start = list(_phases)
n_frames, duration = 20, 1.6
dt = duration / n_frames
for i in range(1, n_frames + 1):
alpha = i / n_frames
_phases[0] = start[0] + alpha * (target[0] - start[0])
_phases[1] = start[1] + alpha * (target[1] - start[1])
draw(_phases)
if i < n_frames:
time.sleep(dt)
def on_phi_slider_change(_change):
if not phase_radio.value.startswith('Random'):
_phases[0] = phi1_slider.value * np.pi
_phases[1] = phi2_slider.value * np.pi
draw(_phases)
w0_slider = widgets.FloatSlider(
value=0.0, min=0.0, max=3.0, step=0.05,
description='w₀ (DC, f=0 Hz)',
style={'description_width': 'initial'},
layout=widgets.Layout(width='400px')
)
w1_slider = widgets.FloatSlider(
value=0.0, min=0.0, max=3.0, step=0.05,
description=f'w₁ (f={F1} Hz)',
style={'description_width': 'initial'},
layout=widgets.Layout(width='400px')
)
w2_slider = widgets.FloatSlider(
value=0.0, min=0.0, max=3.0, step=0.05,
description=f'w₂ (f={F2} Hz)',
style={'description_width': 'initial'},
layout=widgets.Layout(width='400px')
)
phi1_slider = widgets.FloatSlider(
value=0.0, min=-1.0, max=1.0, step=0.05,
description='φ₁ (×π)',
style={'description_width': 'initial'},
layout=widgets.Layout(width='400px'),
readout_format='.2f'
)
phi0_slider = widgets.FloatSlider(
value=0.0, min=-1.0, max=1.0, step=0.05,
description='φ₀ (×π)',
style={'description_width': 'initial'},
layout=widgets.Layout(width='400px'),
readout_format='.2f'
)
phi2_slider = widgets.FloatSlider(
value=0.0, min=-1.0, max=1.0, step=0.05,
description='φ₂ (×π)',
style={'description_width': 'initial'},
layout=widgets.Layout(width='400px'),
readout_format='.2f'
)
phase_radio = widgets.RadioButtons(
options=['Manual phase (use φ sliders)', 'Random phase (fixed per session)'],
value='Manual phase (use φ sliders)',
description='Phase:',
style={'description_width': 'initial'}
)
show_answer_checkbox = widgets.Checkbox(
value=False,
description='Show answer (in complex coefficients plot)',
indent=False
)
w0_slider.observe(update, names='value')
w1_slider.observe(update, names='value')
w2_slider.observe(update, names='value')
phi0_slider.observe(update, names='value')
phi1_slider.observe(on_phi_slider_change, names='value')
phi2_slider.observe(on_phi_slider_change, names='value')
phase_radio.observe(on_phase_change, names='value')
show_answer_checkbox.observe(update, names='value')
ui = widgets.VBox([
widgets.HBox([w0_slider, w1_slider, w2_slider]),
widgets.HBox([
phi0_slider,
phi1_slider,
phi2_slider
]),
phase_radio,
show_answer_checkbox,
img_out
])
display(ui)
update()
build_decomposition_demo()
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Demo_1_Decomposition")
Section 2: The Fourier series#
Estimated timing to here from start of tutorial: 10 min
Video 2: Building complex signals from sine waves#
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Video_2_FourierSeries")
Using the Fourier basis to approximate any signal#
In the previous demo, we made a ground truth signal from sine waves and reconstructed it with sine waves. Maybe not super impressive.
But, as it turns out, any signal — even ones that look nothing like sine waves, like square waves and sawtooth waves — can also be represented as a sum of sine waves:
where the frequencies \(f_1, f_2, \ldots\) are fixed harmonics of a fundamental frequency \(f_0\), and the \(a_n\) are the Fourier coefficients that determine how much of each frequency is present.
This works because sine waves that satisfy certain conditions compose a set of orthogonal basis functions, conceptually similar to the orthogonal basis vectors you learned about in the Dimensionality Reduction tutorials.
Check out the Bonus section for a deeper dive into the sinusoidal basis and its orthogonality.
Interactive Demo 2: Building signals from (a few) sine waves#
Here we look only at the 5 components — enough to recognise the shape of each signal.
Your task: choose the coefficient for each component using the sliders and try to match the target signal (gray dashed) — the actual square / sawtooth wave.
With only 5 terms you cannot drive the MSE all the way to 0; instead, aim for the best possible 5-term fit. When your coefficients reach it the badge turns green, and the MSE it shows is the leftover truncation error — the gap between 5 sine waves and the true wave. (You’ll close that gap by adding more terms in the next section.)
To make it easier, there is no phase shift, i.e., all components have \(\phi = 0\).
But notice that not all magnitude coefficients are positive — some components are subtracted rather than added. Look at what happens to the shape in the bottom panel as you go negative.
A negative coefficient is equivalent to a positive one with a phase shift of \(\pi\) — the component is simply flipped upside-down. So the sign of a slider tells you the phase: positive means \(\phi = 0\), negative means \(\phi = \pi\). In the magnitude spectrum the bar height is always the positive magnitude \(|a_n|\); the sign/phase is shown by the hatch.
Cycle through the two signal types (square vs. sawtooth) and pay attention to the frequency spectra (bottom row) and the signs of the slider values you needed to match each target.
Building signals from sine waves#
Show code cell source
# @title Building signals from sine waves
def build_fourier_builder_demo():
# -- Instructor config --
F0 = 2.0
DURATION, FS = 1.0, 2000
t = np.linspace(0, DURATION, int(FS * DURATION), endpoint=False)
rng = np.random.default_rng()
RAND_PHASES = rng.uniform(0, 2 * np.pi, 5)
COMP_COLORS = ['#e85d04', '#7c3aed', '#059669', '#dc2626', "#ffee00"]
SIGNAL_DEFS = {
'Square wave': {
'harmonics': [1, 3, 5, 7, 9],
'coeffs': [4 / (n * np.pi) for n in [1, 3, 5, 7, 9]],
},
'Sawtooth wave': {
'harmonics': [1, 2, 3, 4, 5],
'coeffs': [2 * (-1)**(n + 1) / (n * np.pi) for n in [1, 2, 3, 4, 5]],
},
}
def true_signal(signal_type):
"""Analytic ground-truth waveform (amplitude +/-1) -- the *actual* wave,
not a Fourier sum. This is what the student's reconstruction is scored
against, so the displayed MSE never quite reaches 0 with only 5 terms."""
x = 2 * np.pi * F0 * t
if signal_type == 'Square wave':
return np.where(np.sin(x) >= 0, 1.0, -1.0)
else: # Sawtooth wave
return 2 * ((F0 * t + 0.5) % 1.0) - 1.0
_phases = [0.0] * 5
img_out = widgets.Image(format='png')
def build_signal(coeffs, freqs, phases):
return sum(a * np.sin(2 * np.pi * f * t + phi)
for a, f, phi in zip(coeffs, freqs, phases))
def draw(phases):
signal_type = signal_dropdown.value
defn = SIGNAL_DEFS[signal_type]
harmonics = defn['harmonics']
freqs = [n * F0 for n in harmonics]
true_coeffs = defn['coeffs']
student_coeffs = [s.value for s in sliders]
for slider, n, f in zip(sliders, harmonics, freqs):
slider.description = f'a{n} ({f:.0f} Hz)'
# y_target: the actual analytic wave (gray dashed) -- what we score against.
# y_best: the best achievable 5-term Fourier fit (true coefficients,
# zero phase). Matching it is the goal, and it drives the
# green/red trigger exactly as before.
y_target = true_signal(signal_type)
y_best = build_signal(true_coeffs, freqs, [0.0] * 5)
y_student = build_signal(student_coeffs, freqs, phases)
components = [a * np.sin(2 * np.pi * f * t + phi)
for a, f, phi in zip(student_coeffs, freqs, phases)]
# Displayed MSE: distance to the *actual* wave (floors at the truncation
# residual, so it stays >0 even for a perfect 5-term fit).
mse_display = np.mean((y_target - y_student) ** 2)
# Trigger MSE: distance to the best 5-term fit -- "the way it is now".
mse_trigger = np.mean((y_best - y_student) ** 2)
is_matched = mse_trigger < 0.01
# Phase-independent y-limits: use sum of |coefficients| as upper bound
# so the axis only moves when magnitudes change, not during phase animation
y_abs_max = max(
sum(abs(c) for c in true_coeffs),
sum(abs(c) for c in student_coeffs),
float(np.max(np.abs(y_target))),
0.5
) * 1.2
ylim = (-y_abs_max, y_abs_max)
# Spectrum/complex panel data
max_h = max(harmonics)
all_h = list(range(1, max_h + 1))
active = dict(zip(harmonics, range(len(harmonics))))
spec_labels = [f'{h*F0:.0f} Hz' for h in all_h]
# Effective phase includes a pi shift for negative coefficients
eff_phases = []
for coeff, phi in zip(student_coeffs, phases):
phi_eff = phi + (np.pi if coeff < 0 else 0.0)
eff_phases.append(np.angle(np.exp(1j * phi_eff)))
spec_heights = []
spec_colors = []
spec_hatches = []
for h in all_h:
if h in active:
coeff = student_coeffs[active[h]]
spec_heights.append(abs(coeff))
spec_colors.append(COMP_COLORS[active[h]])
spec_hatches.append('///' if coeff < 0 else '')
else:
spec_heights.append(0.0)
spec_colors.append('#dddddd')
spec_hatches.append('')
# FFT-derived reference coefficients at the active harmonics.
# Convert FFT complex coefficients to the sine-form used in this demo:
# x(t) ~= sum a_n sin(2*pi*f_n*t + phi_n), where a_n = 2|X(f_n)|.
answer_amps = []
answer_phases = []
yf = np.fft.rfft(y_target) / len(y_target)
fft_freqs = np.fft.rfftfreq(len(y_target), d=1.0 / FS)
for f in freqs:
idx = int(np.argmin(np.abs(fft_freqs - f)))
c = yf[idx]
amp = 2.0 * np.abs(c)
phi = np.angle(c) + np.pi / 2
answer_amps.append(amp)
answer_phases.append(np.angle(np.exp(1j * phi)))
fig = plt.figure(figsize=(14, 9))
gs = fig.add_gridspec(3, 3, height_ratios=[1, 1, 1.5])
ax_main = fig.add_subplot(gs[0, :])
ax_comp = fig.add_subplot(gs[1, :])
ax_complex = fig.add_subplot(gs[2, 0])
ax_mag_spec = fig.add_subplot(gs[2, 1])
ax_phase_spec = fig.add_subplot(gs[2, 2])
# -- Top: reconstruction vs the true wave --
ax_main.plot(t, y_target, color='#888888', linestyle='--', linewidth=2,
label=f'Target (true {signal_type.lower()})', alpha=0.85)
ax_main.plot(t, y_student, color='#2563EB', linewidth=2,
label='Your reconstruction')
ax_main.set_title(f'{signal_type} -- Reconstruction', fontsize=13)
ax_main.set_xlabel('Time (s)')
ax_main.set_ylabel('Amplitude')
ax_main.legend(loc='upper right')
ax_main.set_xlim(0, DURATION)
ax_main.set_ylim(ylim)
ax_main.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
mse_color = '#16a34a' if is_matched else '#dc2626'
ax_main.text(0.02, 0.95, f'MSE: {mse_display:.4f}',
transform=ax_main.transAxes, fontsize=12, verticalalignment='top',
color=mse_color,
bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
edgecolor=mse_color, alpha=0.8))
# -- Middle: components overlaid --
phi_mode = phase_radio.value.startswith('Random')
for comp, f, coeff, phi_eff, n in zip(components, freqs, student_coeffs, eff_phases, harmonics):
phi_lbl = f' + {phi_eff/np.pi:.2f}pi' if phi_mode and abs(phi_eff) > 1e-3 else ''
ax_comp.plot(t, comp, color=COMP_COLORS[active[n]], linewidth=1.8,
label=f'a{n}={coeff:.2f}*sin(2pi*{f:.0f}Hz*t{phi_lbl})')
ax_comp.set_title('Components', fontsize=12)
ax_comp.set_xlabel('Time (s)')
ax_comp.set_ylabel('Amplitude')
ax_comp.legend(loc='upper right', fontsize=8)
ax_comp.set_xlim(0, DURATION)
ax_comp.set_ylim(ylim)
ax_comp.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
# -- Bottom-left: complex coefficient view --
rmax = max(max(abs(c) for c in student_coeffs), max(answer_amps), 0.5) * 1.25
theta = np.linspace(-np.pi, np.pi, 300)
for coeff, phi_eff, color in zip(student_coeffs, eff_phases, COMP_COLORS):
r = abs(coeff)
ax_complex.plot(r * np.cos(theta), r * np.sin(theta),
color=color, linewidth=1.6, alpha=0.9)
px, py = r * np.cos(phi_eff), r * np.sin(phi_eff)
ax_complex.scatter([px], [py], color=color, s=50, zorder=3)
ax_complex.plot([0, px], [0, py], color=color, linewidth=1.1, alpha=0.7)
if show_answer_checkbox.value:
for amp, phi, color in zip(answer_amps, answer_phases, COMP_COLORS):
ax, ay = amp * np.cos(phi), amp * np.sin(phi)
ax_complex.scatter([ax], [ay], marker='*', s=220, color=color,
edgecolors='black', linewidths=0.8, zorder=6)
ax_complex.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_complex.axvline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_complex.set_xlim(-rmax, rmax)
ax_complex.set_ylim(-rmax, rmax)
ax_complex.set_aspect('equal', adjustable='box')
ax_complex.set_title('Complex Coefficients', fontsize=12)
ax_complex.set_xlabel('Real')
ax_complex.set_ylabel('Imag')
ax_complex.grid(alpha=0.2)
# -- Bottom-middle: magnitude spectrum --
for h, height, color, hatch in zip(all_h, spec_heights, spec_colors, spec_hatches):
ax_mag_spec.bar(h, height, color=color, width=0.55, alpha=0.85,
hatch=hatch, edgecolor='#444444' if hatch else color,
linewidth=0.8)
legend_handles = [
Patch(facecolor='#aaaaaa', edgecolor='#444444', label='positive amplitude'),
Patch(facecolor='#aaaaaa', edgecolor='#444444', hatch='///',
label='negative amplitude (phase +pi)'),
]
ax_mag_spec.legend(handles=legend_handles, fontsize=8, loc='upper right')
ax_mag_spec.set_xticks(all_h)
ax_mag_spec.set_xticklabels(spec_labels, fontsize=8)
ax_mag_spec.set_title('Magnitude Spectrum', fontsize=12)
ax_mag_spec.set_ylabel('Magnitude |w|')
ax_mag_spec.set_ylim(0, max(2.2, max(spec_heights) * 1.3))
ax_mag_spec.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
# -- Bottom-right: phase spectrum --
phase_x = harmonics
phase_vals = eff_phases
phase_colors = [COMP_COLORS[active[h]] for h in harmonics]
ax_phase_spec.vlines(phase_x, 0, phase_vals, colors=phase_colors, linewidth=1.6, alpha=0.8)
ax_phase_spec.scatter(phase_x, phase_vals, color=phase_colors, s=45, zorder=3)
ax_phase_spec.set_xticks(all_h)
ax_phase_spec.set_xticklabels(spec_labels, fontsize=8)
ax_phase_spec.set_title('Phase Spectrum', fontsize=12)
ax_phase_spec.set_ylabel('Phase (rad)')
ax_phase_spec.set_ylim(-np.pi, np.pi)
ax_phase_spec.set_yticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi])
ax_phase_spec.set_yticklabels(['-pi', '-pi/2', '0', 'pi/2', 'pi'])
ax_phase_spec.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
fig.tight_layout()
buf = io.BytesIO()
fig.savefig(buf, format='png', dpi=100, bbox_inches='tight')
plt.close(fig)
buf.seek(0)
img_out.value = buf.read()
def update(_change=None):
draw(_phases)
def on_phase_change(change):
target = list(RAND_PHASES[:5]) if change['new'].startswith('Random') else [0.0] * 5
n_frames, duration = 20, 1.6
dt = duration / n_frames
start = list(_phases)
for i in range(1, n_frames + 1):
alpha = i / n_frames
for j in range(5):
_phases[j] = start[j] + alpha * (target[j] - start[j])
draw(_phases)
if i < n_frames:
time.sleep(dt)
def on_signal_change(change):
for slider in sliders:
slider.unobserve(update, names='value')
for slider in sliders:
slider.value = 0.0
for slider in sliders:
slider.observe(update, names='value')
update()
def on_reset(_btn):
for slider in sliders:
slider.unobserve(update, names='value')
for slider in sliders:
slider.value = 0.0
for slider in sliders:
slider.observe(update, names='value')
update()
signal_dropdown = widgets.Dropdown(
options=list(SIGNAL_DEFS.keys()),
value='Square wave',
description='Signal type:',
style={'description_width': 'initial'},
layout=widgets.Layout(width='280px')
)
sliders = [
widgets.FloatSlider(
value=0.0, min=-2.0, max=2.0, step=0.01,
description='',
style={'description_width': 'initial'},
layout=widgets.Layout(width='550px'),
readout_format='.2f'
)
for i in range(5)
]
phase_radio = widgets.RadioButtons(
options=['No phase (phi = 0)', 'Random phase (fixed per session)'],
value='No phase (phi = 0)',
description='Phase:',
style={'description_width': 'initial'}
)
show_answer_checkbox = widgets.Checkbox(
value=False,
description='Show answer (in complex coefficients)',
indent=False
)
reset_btn = widgets.Button(
description='Reset coefficients',
button_style='warning',
layout=widgets.Layout(width='180px')
)
signal_dropdown.observe(on_signal_change, names='value')
phase_radio.observe(on_phase_change, names='value')
show_answer_checkbox.observe(update, names='value')
reset_btn.on_click(on_reset)
for slider in sliders:
slider.observe(update, names='value')
ui = widgets.VBox([
widgets.HBox([signal_dropdown, reset_btn]),
phase_radio,
show_answer_checkbox,
widgets.VBox(sliders),
img_out
])
display(ui)
update()
build_fourier_builder_demo()
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Demo_2_FourierBuilder")
Think! 1: Look at the spectra#
1. What do you observe about the pattern of signs on the weights? Are all the coefficients positive, or do some need to be negative? Is there a regular pattern (e.g. + + + + or + − + −)? Remember that a negative coefficient is not a different kind of component — it is just the same sine wave with a phase of π (i.e. flipped upside-down). So the sign pattern is really a phase pattern.
2. What do you observe about the magnitudes of the weights? Are the weights all the same size, or do they shrink as the frequency increases? If they shrink, do they decrease linearly (like 1, 1/2, 1/3, 1/4, …) or faster (like 1, 1/9, 1/25, …)? This rate of decay is what determines how quickly the Fourier series converges — and it’s the reason some waveforms need many more harmonics than others to look “right.”
Keep your answers in mind as you move into the next section, where we add more and more components and watch the approximation improve.
Section 3: More Fourier components#
Estimated timing to here from start of tutorial: 20 min
So far you’ve been working with just 5 Fourier components. But in theory you can use many more — so how much does adding more actually help?
Interactive Demo 3: Convergence of the Fourier series#
Use the slider below to increase N, the number of sine wave components included in the approximation. The coefficients are set automatically to the correct Fourier coefficients (magnitude and phases) — you’re just controlling how many terms to keep.
The approximation (blue) is compared against the true signal (gray dashed) — the actual square, sawtooth, or Gaussian pulse.
A few things to watch for as you increase N:
Sawtooth wave converges slowly — its coefficients shrink only as 1/n.
Square wave is the most stubborn: you’ll notice the “ringing” overshoot near each edge (Gibbs phenomenon) that never fully disappears, no matter how many terms you add — it just gets narrower.
Gaussian pulse is something new — an aperiodic signal: a single localized bump rather than a repeating wave. Yet it too is just a sum of plain sine waves. A sharp, narrow pulse needs many harmonics (a broad band of frequencies) to build up; with only a few terms it stays a low, wide blob.
In this demo, we skip the complex coefficient plot and phase spectrum plot, and just show you the power spectrum:
From amplitude to power#
In the previous demos, we have always shown an amplitude (or magnitude) spectrum — the size \(|a_n|\) of each component. Here we switch to the power spectrum, which is just the square of the amplitude: the power carried by each component, \(P_n = |a_n|^2/2\) (the average power of a sine of amplitude \(a_n\)).
Power is what most measurements — and our ears — actually register, and squaring exaggerates the tall components and flattens the small ones, so the shape of the spectrum stands out.
Power is what the convergence metric has been measuring: By Parseval’s theorem, the total power of a signal equals the sum of the power of its Fourier components — so the mean-squared error between the N-term approximation and the true signal is exactly the total power of the components you left out.
Interactive Demo 3: Fourier series convergence#
Show code cell source
# @title Interactive Demo 3: Fourier series convergence
def build_convergence_demo():
# ── Instructor config ──
F0 = 2.0 # fundamental of the square / sawtooth (Hz)
DURATION, FS = 1.0, 2000
t = np.linspace(0, DURATION, int(FS * DURATION), endpoint=False)
GAUSS_F0 = 1.0 / DURATION # one pulse per window → a single localised bump
# Fixed random phases so that raising N just adds more scrambled terms.
_rng = np.random.default_rng(0)
RAND_PHASES = _rng.uniform(0, 2 * np.pi, 256)
# Current Gaussian phases (interpolated during the pulse ↔ noise animation)
_phases = []
def gauss_sigma():
"""Gaussian pulse width in seconds, read live from the width slider."""
return gauss_width_slider.value / 1000.0
def true_signal(signal_type):
"""Analytic ground-truth waveform."""
x = 2 * np.pi * F0 * t
if signal_type == 'Square wave':
return np.where(np.sin(x) >= 0, 1.0, -1.0)
elif signal_type == 'Sawtooth wave':
return 2 * ((F0 * t + 0.5) % 1.0) - 1.0
else: # Gaussian pulse — a single narrow, centred bump (an aperiodic transient)
sigma = gauss_sigma()
return np.exp(-((t - DURATION / 2) ** 2) / (2 * sigma ** 2))
def spectrum(signal_type, N):
"""First-N-harmonic description of the signal.
Every signal is written in the single form
x(t) = DC + Σ |aₙ|·sin(2π fₙ t + φₙ),
so the phases can later be swapped for random values *without touching
the magnitudes* — the heart of the Gaussian phase demo. Returns the
present harmonics, their frequencies, magnitudes |aₙ|, the aligned
phases, the DC term, the fundamental, and the full harmonic grid (so the
square's missing even harmonics can be marked as zeros).
"""
if signal_type == 'Square wave':
harm = [2 * k + 1 for k in range(N)]
mags = [4 / (h * np.pi) for h in harm]
phases = [0.0 for _ in harm]
dc, fundamental = 0.0, F0
grid = list(range(1, harm[-1] + 1))
elif signal_type == 'Sawtooth wave':
harm = list(range(1, N + 1))
mags = [2 / (h * np.pi) for h in harm]
phases = [0.0 if h % 2 == 1 else np.pi for h in harm] # sign (-1)^(h+1)
dc, fundamental = 0.0, F0
grid = list(harm)
else: # Gaussian pulse — true coefficients straight from the FFT
g = true_signal('Gaussian pulse')
coeffs = np.fft.rfft(g) / len(g)
harm = list(range(1, N + 1))
mags = [2 * np.abs(coeffs[h]) for h in harm]
phases = [np.angle(coeffs[h]) + np.pi / 2 for h in harm] # cos → sin
dc, fundamental = float(coeffs[0].real), GAUSS_F0
grid = list(harm)
freqs = [h * fundamental for h in harm]
return harm, freqs, mags, phases, dc, fundamental, grid
def envelope_power(signal_type, nu):
"""Continuous power envelope |X(ν)|²/2 that the harmonics sample, as a
function of (continuous) harmonic number ν. This is the transform of one
period: a sinc² for the square (boxcar edges → nulls + side lobes), a
smooth 1/ν² roll-off for the sawtooth ramp, and a Gaussian for the
Gaussian pulse. The discrete harmonic powers land exactly on it."""
if signal_type == 'Square wave':
return 2.0 * np.sinc(nu / 2.0) ** 2 # np.sinc(x)=sin(πx)/(πx)
elif signal_type == 'Sawtooth wave':
return 2.0 / (np.pi ** 2 * nu ** 2)
else: # Gaussian pulse — width set by the slider
sigma = gauss_sigma()
return 4 * np.pi * sigma ** 2 * np.exp(-4 * np.pi ** 2 * sigma ** 2 * nu ** 2)
def reconstruct(freqs, mags, phases, dc):
y = np.full_like(t, dc)
for f, m, p in zip(freqs, mags, phases):
y = y + m * np.sin(2 * np.pi * f * t + p)
return y
def target_phases(N):
"""The Gaussian phase vector for the *current* radio mode at this N."""
_, _, _, aligned, _, _, _ = spectrum('Gaussian pulse', N)
if phase_radio.value.startswith('Random'):
return list(RAND_PHASES[:N])
return list(aligned)
conv_dropdown = widgets.Dropdown(
options=['Square wave', 'Sawtooth wave', 'Gaussian pulse'],
value='Square wave',
description='Signal type:',
style={'description_width': 'initial'},
layout=widgets.Layout(width='280px')
)
n_slider = widgets.IntSlider(
value=1, min=1, max=200, step=1,
description='N components:',
style={'description_width': 'initial'},
layout=widgets.Layout(width='500px')
)
phase_radio = widgets.RadioButtons(
options=['Aligned phase → pulse', 'Random phase → noise'],
value='Aligned phase → pulse',
description='Gaussian phase:',
style={'description_width': 'initial'},
disabled=True
)
scale_toggle = widgets.ToggleButtons(
options=['dB', 'Log-log'],
value='dB',
description='Power and frequency axis:',
style={'description_width': 'initial'}
)
gauss_width_slider = widgets.IntSlider(
value=8, min=4, max=30, step=1,
description='Gaussian width σ (ms):',
style={'description_width': 'initial'},
layout=widgets.Layout(width='500px'),
disabled=True
)
conv_img = widgets.Image(format='png')
def draw(gauss_phases):
signal_type = conv_dropdown.value
N = n_slider.value
is_gauss = signal_type == 'Gaussian pulse'
harm, freqs, mags, phases, dc, fundamental, grid = spectrum(signal_type, N)
recon_phases = gauss_phases if is_gauss else phases
randomized = is_gauss and phase_radio.value.startswith('Random')
y_ref = true_signal(signal_type)
y_approx = reconstruct(freqs, mags, recon_phases, dc)
mse = np.mean((y_ref - y_approx) ** 2)
max_h = max(grid)
fig = plt.figure(figsize=(10, 7.5), constrained_layout=True)
gs = fig.add_gridspec(2, 1, height_ratios=[3, 2])
ax = fig.add_subplot(gs[0])
ax_spec = fig.add_subplot(gs[1])
# ── Top: reconstruction vs the true signal ──
approx_lbl = (f'{N}-term sum (phases scrambled)' if randomized
else f'{N}-term approximation')
ax.plot(t, y_ref, color='#888888', linestyle='--', linewidth=2,
label=f'Ground truth {signal_type.lower()}', alpha=0.85)
ax.plot(t, y_approx, color='#2563EB', linewidth=2, label=approx_lbl)
gauss_tail = f' — σ = {gauss_width_slider.value} ms' if is_gauss else ''
title_tail = ' — phases scrambled → noise' if randomized else gauss_tail
ax.set_title(f'{signal_type} — Fourier Series (N = {N}){title_tail}', fontsize=13)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.legend(loc='upper right')
ax.set_xlim(0, DURATION)
ax.set_ylim((-0.6, 1.25) if is_gauss else (-1.35, 1.35))
ax.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
mse_color = '#16a34a' if mse < 0.005 else '#dc2626'
ax.text(0.02, 0.95, f'MSE vs ground truth: {mse:.5f}',
transform=ax.transAxes, fontsize=11, verticalalignment='top',
color=mse_color,
bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
edgecolor=mse_color, alpha=0.8))
# ── Bottom: POWER spectrum — discrete harmonics sampling a continuous
# envelope (the transform of one period) ──
powers = [m ** 2 / 2 for m in mags] # Pₙ = |aₙ|²/2
missing_h = [h for h in grid if h not in set(harm)]
missing_freq = [h * fundamental for h in missing_h]
nu = np.linspace(0.5, max_h + 1, 800)
env = envelope_power(signal_type, nu)
env_freq = nu * fundamental
pmax = max(powers) if powers else 1.0
mode = scale_toggle.value
if mode == 'dB':
FLOOR = -40.0
to_db = lambda P: np.clip(10 * np.log10(np.maximum(np.asarray(P, float), 1e-12) / pmax),
FLOOR, None)
env_y = to_db(env)
present_y = to_db(powers)
missing_y = np.full(len(missing_freq), FLOOR)
base = FLOOR
ylabel = 'Power (dB)'
xlim = (0, (max_h + 1) * fundamental)
ylim = (FLOOR, 5)
else: # Log-log
floor = pmax * 1e-5
env_y = np.maximum(env, floor)
present_y = np.maximum(np.asarray(powers, float), floor)
missing_y = np.full(len(missing_freq), floor)
base = floor
ylabel = 'Power (log)'
xlim = (min(freqs) * 0.8, (max_h + 1) * fundamental)
ylim = (floor, pmax * 2)
ax_spec.plot(env_freq, env_y, color='#9aa0a6', linewidth=1.6, zorder=1,
label='envelope (transform of one period)')
ax_spec.vlines(freqs, base, present_y, color='#2563EB', linewidth=1.2,
alpha=0.55, zorder=2)
ax_spec.scatter(freqs, present_y, color='#2563EB', s=16, zorder=3,
label='harmonic')
if missing_freq:
ax_spec.scatter(missing_freq, missing_y, marker='x', color='#888888',
s=26, linewidths=1.2, zorder=3,
label='zero power (between harmonics)')
spec_note = ' (identical for pulse and noise)' if is_gauss else ''
ax_spec.set_title('Power Spectrum of the Components' + spec_note, fontsize=12)
ax_spec.set_xlabel('Frequency (Hz)')
ax_spec.set_ylabel(ylabel)
if mode == 'Log-log':
ax_spec.set_xscale('log')
ax_spec.set_yscale('log')
ax_spec.set_xlim(xlim)
ax_spec.set_ylim(ylim)
ax_spec.axhline(base, color='#cccccc', linewidth=0.8, zorder=0)
ax_spec.legend(loc='upper right', fontsize=8)
if mode == 'dB':
# Thinned frequency ticks so large N stays readable (linear-freq view)
step = max(1, len(grid) // 10)
tick_h = grid[::step]
ax_spec.set_xticks([h * fundamental for h in tick_h])
ax_spec.set_xticklabels([f'{h*fundamental:.0f}' for h in tick_h], fontsize=8)
buf = io.BytesIO()
fig.savefig(buf, format='png', dpi=100, bbox_inches='tight')
plt.close(fig)
buf.seek(0)
conv_img.value = buf.read()
def update_conv(_change=None):
# Static redraw (slider / dropdown / scale / width). Refresh the Gaussian
# phases to match the current radio mode at this N; others ignore them.
if conv_dropdown.value == 'Gaussian pulse':
_phases[:] = target_phases(n_slider.value)
draw(_phases)
def on_phase_change(change):
# Animate the Gaussian smoothly between aligned-phase pulse and
# random-phase noise (both directions). The power spectrum never changes
# — only the phases morph.
if conv_dropdown.value != 'Gaussian pulse':
return
N = n_slider.value
target = target_phases(N) # destination (new mode)
start = list(_phases) if len(_phases) == N else target_phases(N)
n_frames, duration = 20, 1.6
dt = duration / n_frames
for i in range(1, n_frames + 1):
a = i / n_frames
_phases[:] = [s + a * (tt - s) for s, tt in zip(start, target)]
draw(_phases)
if i < n_frames:
time.sleep(dt)
def on_signal_change(change):
# The phase + width controls only make sense for the Gaussian pulse
is_gauss = (change['new'] == 'Gaussian pulse')
phase_radio.disabled = not is_gauss
gauss_width_slider.disabled = not is_gauss
update_conv()
conv_dropdown.observe(on_signal_change, names='value')
n_slider.observe(update_conv, names='value')
phase_radio.observe(on_phase_change, names='value')
scale_toggle.observe(update_conv, names='value')
gauss_width_slider.observe(update_conv, names='value')
display(widgets.VBox([
widgets.HBox([conv_dropdown, n_slider]),
widgets.HBox([phase_radio, scale_toggle]),
gauss_width_slider,
conv_img
]))
update_conv()
build_convergence_demo()
Think! 2: Fourier series and convergence#
After playing around with the number of components, we want to point out a few things:
1. Power and frequency axes#
Use the axis toggle to switch how the horizontal and vertical scales are drawn:
dB (decibels — log power against linear frequency), usually used in engineering contexts
Log-log (log power against log frequency), the natural view for power laws, which often appear in neural signals.
2. Reading the envelope in frequency domain#
The components are drawn as stems sampling a smooth gray envelope — the Fourier Transform of one period of the signal. Notice the different shapes of the envelopes for the three signals:
Square wave → a sinc² envelope: a fat main lobe, then nulls (every even harmonic lands exactly in a null, with zero power) and decaying side lobes carrying the odd harmonics. (You’ll see this same sinc again when we get to filtering.)
Sawtooth wave → a smooth 1/f² roll-off, no nulls — every harmonic present, power shrinking steadily.
Gaussian pulse → a Gaussian envelope: smooth, no nulls, falling off fastest of all. Try changing the width slider.
3. Time-frequency trade-off#
Changing the width slider of the Gaussian signal highlights the time–frequency trade-off:
larger σ and the pulse spreads out in time while its spectrum narrows — a broad, smooth bump is built from just a few low frequencies
smaller σ and the pulse sharpens in time while its spectrum broadens toward higher frequencies
Time-width and frequency-width move in opposite directions (\(\Delta t \cdot \Delta f \approx\) constant) — the same reciprocity behind the uncertainty principle.
4. Phase: same power spectrum, very different signal#
Select Gaussian pulse, then toggle Aligned phase → pulse against Random phase → noise. This changes only the phases of the components, and it animates: you’ll watch the sharp pulse smoothly melt into noise (and morph back) while the power spectrum stays exactly the same — the bottom panel never moves, yet the signal completely changes:
Aligned: all the sine waves peak at the same instant, add up constructively, and cancel elsewhere → a single sharp pulse.
Random: the same amplitudes now peak at scattered times, so the energy smears across the whole window → noise.
The power spectrum alone does not determine a signal. It tells you which frequencies and how much — but it’s the phase that says when things happen. A spike and a burst of noise can have an identical spectrum and differ only in their phases.
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Demo_3_Convergence")
Section 4: Discrete Fourier Transform (DFT) and Fast Fourier Transform (FFT)#
Estimated timing to here from start of tutorial: 25 min
We wrap up the theory of Fourier analysis with the Discrete Fourier Transform (DFT), and build it from scratch:
“Fourier analysis” actually refers to multiple methods that convert a function or signal into its frequency-domain representation using the sinusoidal basis functions, with different assumptions and conditions.
But in almost all cases, you will be working with the Discrete Fourier Transform (DFT), which operates on discrete or sampled signals (like we saw in Tutorial 1). All signals in the computer, neural or otherwise, are inherently discrete due to the nature of data acquisition.
The Fast Fourier Transform (FFT) is a specific implementation of the DFT using a more complicated but efficient algorithm (Cooley–Tukey or butterfly algorithm) to compute the Fourier coefficients, which is widely used in practice in any digital signal processing applications.
In the first and second demo, we asked you to manually find the best-fitting Fourier coefficients, and later showed that increasing the number of components improves the approximation in Demo 3. Of course, you would never do this by hand. But we haven’t shown you:
How to determine the number of Fourier components to use, and
How to compute the coefficients for each component, i.e., the amplitude (or power) and phase of each frequency.
The DFT equation answers both of these questions, and we will implement the naive formulation from scratch in the final exercise.
Video 3: Building the DFT from scratch#
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Video_3_DFT")
Discrete Fourier Transform (DFT)#
Definition: for a sequence of N numbers \(x_0, x_1, …, x_{N-1}\) (like a signal with N samples), the DFT transforms it into another sequence of N numbers \(X_0, X_1, ..., X_{N-1}\), where:
This might look scary, but you’ve seen everything in it before:
\(X_k\) are the Fourier coefficients, which is a complex number that encodes both the amplitude and phase of the \(k\)-th component, exactly as we have been seeing in the (bottom left) complex coefficient panel in Demos 1 and 2. This set of \(X_k\) is our target to compute.
\(x_n\) are the values of the signal over time, \(N\) is the total number of samples, and \(i\) is the imaginary unit.
To answer the two questions we posed above:
The number of Fourier components for a sampled signal is equal to the number of samples in the signal, \(N\), as integer multiples \(k = 0, 1, ..., N-1\) of the lowest frequency \(f_0 = \frac{1}{N}\). \(k\) is often taken to be \(-N/2, ... 0, …, N/2-1\) instead, centering on the zero component and ending at the Nyquist frequency. For a real-valued signal, the negative frequency components are redundant, so we usually only look at the first half of the spectrum (i.e., \(k = 0, 1, ..., N/2\) up to the Nyquist frequency).
The coefficients \(X_k\) are computed using the DFT equation above, which gives you the amplitude and phase of each frequency component. It’s a dot product between the signal and another vector of the same length.
To compute the DFT coefficients (and the power spectrum) in practice, we use the Fast Fourier Transform (FFT) by calling numpy.fft.fft or scipy.fft.fft in Python. This decomposition is exact and invertible, meaning that you can reconstruct the original signal exactly from its Fourier coefficients, using the inverse DFT by calling numpy.fft.ifft or scipy.fft.ifft.
Complex exponentials#
The weird-looking exponential term \(e^{-i 2\pi \frac{k}{N}n}\) is the fundamental building block of the DFT — this represents the complex sinusoidal basis functions we discussed earlier.
Euler’s formula states that \(e^{i\theta} = \cos(\theta) + i\sin(\theta)\), which means that this complex exponential in the equation can be thought of as a combination of a cosine and (imaginary) sine wave.
Replacing the argument \(\theta\) with \(-2\pi \frac{k}{N}n\), the DFT formula can be rewritten as: $\(X_k = \sum_{n=0}^{N-1} x_n \left[ \cos\left(2\pi \frac{k}{N}n\right) - i\sin\left(2\pi \frac{k}{N}n\right) \right]\)$
You have seen exactly this in the complex coefficient plots: instead of writing the sine wave with an amplitude and phase, we can write it as a sine and cosine pair with the same amplitude. You might also recognize the polar form of a complex number, \(z = r e^{i\phi}\), where \(r\) is the magnitude and \(\phi\) is the phase. The DFT coefficients \(X_k\) are complex numbers that encode both the amplitude and phase of each frequency component:
Code Exercise 1: Implement the DFT from scratch#
We can implement the DFT equation directly in code, using the expanded formula with cosine and sine:
There are just 4 steps:
create_cos_sin(N, k): For each frequency \(k\), compute the complex exponential \(e^{-i 2\pi \frac{k}{N}n}\) for all time points \(n\). We do this by making the cosine and sine terms.dot_product(x, y): For each frequency \(k\), multiply the complex exponential by the signal values \(x_n\) and sum over all time points \(n\) to get the Fourier coefficient \(X_k\). This is just a dot product between the signal and the complex exponential of frequency \(k\).compute_Xk(signal, k): now you can compute the DFT coefficient \(X_k\) for any frequency \(k\) by calling the two functions above to make inner term and take the sum overn.compute_DFT(signal): Finally, you can compute the full DFT by iterating over all frequency indices \(k = 0, 1, …, N-1\) and storing the results in an array.
def create_cos_sin(N, k):
"""
Create the k-th cosine and sine basis vectors of length N.
N: number of samples in your signal.
k: frequency index or "wave number" (0 <= k < N).
"""
####################################################################################
## TODO for students: create cosine and sine basis vectors.
# Fill out the function and remove the line below.
raise NotImplementedError("Student exercise: create cosine and sine basis vectors")
####################################################################################
n = np.arange(N)
cos_k = ...
sin_k = ...
return cos_k, sin_k
def dot_product(x, y):
"""
Compute the dot product of two vectors x and y.
Identical to np.dot(x,y) but implemented manually.
x: first vector.
y: second vector.
"""
####################################################################################
## TODO for students: compute dot product.
# Fill out the function and remove the line below.
raise NotImplementedError("Student exercise: compute dot product")
####################################################################################
return ...
def compute_Xk(x, k):
"""
Compute the k-th DFT coefficient X[k] for a signal x.
x: input signal (1D array).
k: frequency index (0 <= k < N).
"""
####################################################################################
## TODO for students: compute X[k].
# Fill out the function and remove the line below.
raise NotImplementedError("Student exercise: compute X[k]")
####################################################################################
N = len(x)
cos_k, sin_k = create_cos_sin(...)
# Compute the real and imaginary parts of X[k] separately, each with a dot product
real_part = ...
imag_part = ...
# Combine into a complex number, note the negative sign for the imaginary part
X_k = ... - 1j * ...
return X_k
def compute_DFT(x):
"""
Compute the Discrete Fourier Transform (DFT) of a signal x.
x: input signal (1D array).
"""
####################################################################################
## TODO for students: compute DFT.
# Fill out the function and remove the line below.
raise NotImplementedError("Student exercise: compute DFT")
####################################################################################
N = ...
X = np.zeros(N, dtype=complex)
# Loop over all frequency indices k
for k in ...:
X[k] = ...
return X
To test your functions individually as you write each one, uncomment each assert as you go.#
# Make 1 second of square wave of 5 Hz, with 1000 Hz rate as example:
FS = 1000 # signal sampling rate
F0 = 5.0 # square wave frequency
DURATION = 1.0 # seconds
N = int(FS * DURATION) # number of samples in the signal
t = np.arange(0, DURATION, 1/FS)
signal = np.where(np.sin(2 * np.pi * F0 * t + 0.25 * np.pi) >= 0, 1.0, -1.0)
# Compute "ground truth reference" DFT using numpy fft
X_fft = np.fft.fft(signal)
# Test by computing the 5th DFT coefficient (k=5) of the signal
k = 2
# 1. Test create_cos_sin function
cos_basis, sin_basis = create_cos_sin(N, k)
assert np.isclose(cos_basis[2], np.cos(2 * np.pi * k * 2 / N))
assert np.isclose(sin_basis[2], np.sin(2 * np.pi * k * 2 / N))
# 2. Test dot_product function
assert np.isclose(dot_product(cos_basis, signal), np.dot(cos_basis, signal))
# 3. Test compute_Xk function
X_k = compute_Xk(signal, k)
assert np.isclose(X_k, X_fft[k])
# 4. Test compute_DFT function by computing all X_ks
X_dft = compute_DFT(signal)
assert np.allclose(X_dft, X_fft)
print("All tests passed! Your DFT implementation is correct.")
If your DFT implementation is correct, you can visualize it here exactly in the same way we did in the previous demos, and compare it against the built-in numpy.fft.fft function.#
Interactive Demo: Visualize your DFT implementation#
Show code cell source
# @title Interactive Demo: Visualize your DFT implementation
def visualize_single_k_dft_demo():
FS = 1000 # signal sampling rate (Hz)
F0 = 2.0 # square-wave frequency (Hz)
DURATION = 1.0 # seconds
N = int(FS * DURATION)
t = np.arange(0, DURATION, 1 / FS)
signal = np.where(np.sin(2 * np.pi * F0 * t + 0.25 * np.pi) >= 0, 1.0, -1.0)
# Compute both user-implemented DFT and NumPy FFT for comparison.
X_user = compute_DFT(signal)
X_fft = np.fft.fft(signal)
# For a real-valued signal, positive-frequency bins are unique up to Nyquist.
k_max = N // 2
def draw(k=5):
k = int(k)
cos_k, sin_k = create_cos_sin(N, k)
real_k = np.real(X_user[k])
imag_proj_k = -np.imag(X_user[k]) # equals sum_n x[n] * sin(2*pi*k*n/N)
# Component weights from the single complex coefficient X[k] in the IDFT.
# (No factor of 2 here: that factor appears only when combining +/-k pairs.)
cos_comp = (real_k / N) * cos_k
sin_comp = (imag_proj_k / N) * sin_k
freqs = np.arange(0, k_max + 1) * FS / N
X_pos_user = X_user[:k_max + 1]
amp = np.abs(X_pos_user)
phase = np.angle(X_pos_user)
fig = plt.figure(figsize=(14, 8))
gs = fig.add_gridspec(2, 3, height_ratios=[1.2, 1.0])
ax_top = fig.add_subplot(gs[0, :])
ax_complex = fig.add_subplot(gs[1, 0])
ax_amp = fig.add_subplot(gs[1, 1])
ax_phase = fig.add_subplot(gs[1, 2])
# Top row: time-domain signal and the k-th basis contributions.
ax_top.plot(t, signal, color='#888888', linestyle='--', linewidth=2, label='Original signal')
ax_top.plot(t, cos_comp, color='#e85d04', linewidth=2,
label=f'k={k} cosine part ({real_k / N:.3f} * cos)')
ax_top.plot(t, sin_comp, color='#7c3aed', linewidth=2,
label=f'k={k} sine part ({imag_proj_k / N:.3f} * sin)')
ax_top.set_title(f'Time Domain: original signal and k={k} components (f={k * FS / N:.1f} Hz)', fontsize=13)
ax_top.set_xlabel('Time (s)')
ax_top.set_ylabel('Amplitude')
ax_top.set_xlim(0, DURATION)
y_lim = 1.2 * max(np.max(np.abs(signal)), np.max(np.abs(cos_comp)), np.max(np.abs(sin_comp)), 0.5)
ax_top.set_ylim(-y_lim, y_lim)
ax_top.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_top.legend(loc='upper right')
# Bottom-left: complex coefficient comparison at selected k.
xk_user = X_user[k]
xk_fft = X_fft[k]
r_user = np.abs(xk_user)
rmax = max(1.0, 1.25 * max(r_user, np.abs(xk_fft)))
theta = np.linspace(-np.pi, np.pi, 300)
ax_complex.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_complex.axvline(0, color='#cccccc', linewidth=0.8, zorder=0)
# FFT: star marker at the FFT coefficient (drawn behind user overlays).
ax_complex.scatter([np.real(xk_fft)], [np.imag(xk_fft)],
marker='*', s=220, color='#dc2626', edgecolors='black',
linewidths=0.8, zorder=2, label='NumPy FFT X[k]')
# User DFT: circle (magnitude) + component lines + dot (coefficient endpoint).
ax_complex.plot(r_user * np.cos(theta), r_user * np.sin(theta),
color='#1f2937', linewidth=1.8, alpha=0.9, zorder=4,
label=None)
# User coefficient radius and endpoint.
ax_complex.plot([0, np.real(xk_user)], [0, np.imag(xk_user)],
color='#111111', linewidth=2, alpha=0.8, zorder=4)
ax_complex.scatter([np.real(xk_user)], [np.imag(xk_user)],
color='#111111', s=10, alpha=0.75, zorder=5, label='Your DFT X[k]')
# Cosine component (real part): horizontal line in orange.
ax_complex.plot([0, np.real(xk_user)], [0, 0],
color='#e85d04', linewidth=2.0, zorder=4,
label='Cos component (Re)')
# Sine component (imag part): vertical line in purple.
ax_complex.plot([np.real(xk_user), np.real(xk_user)], [0, np.imag(xk_user)],
color='#7c3aed', linewidth=2.0, zorder=4,
label='Sin component (Im)')
ax_complex.set_xlim(-rmax, rmax)
ax_complex.set_ylim(-rmax, rmax)
ax_complex.set_aspect('equal', adjustable='box')
ax_complex.grid(alpha=0.2)
ax_complex.set_title(f'Complex coefficient comparison at k={k}', fontsize=12)
ax_complex.set_xlabel('Real')
ax_complex.set_ylabel('Imag')
ax_complex.legend(loc='upper right', fontsize=8)
# Bottom-middle: raw amplitude spectrum from user DFT coefficients.
ax_amp.plot(freqs, amp, color='#6b7280', linewidth=1.5)
ax_amp.scatter([freqs[k]], [amp[k]], color='#dc2626', s=55, zorder=3, label=f'k={k}')
ax_amp.set_title('Amplitude Spectrum', fontsize=12)
ax_amp.set_xlabel('Frequency (Hz)')
ax_amp.set_ylabel('|X[k]|')
ax_amp.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_amp.set_xlim(freqs[0], 100)
ax_amp.legend(loc='upper right', fontsize=9)
# Bottom-right: phase spectrum with selected k highlighted.
ax_phase.plot(freqs, phase, color='#6b7280', linewidth=1.0)
ax_phase.scatter([freqs[k]], [phase[k]], color='#dc2626', s=55, zorder=3, label=f'k={k}')
ax_phase.set_title('Phase Spectrum', fontsize=12)
ax_phase.set_xlabel('Frequency (Hz)')
ax_phase.set_ylabel('Phase (rad)')
ax_phase.set_ylim(-np.pi, np.pi)
ax_phase.set_yticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi])
ax_phase.set_yticklabels(['-pi', '-pi/2', '0', 'pi/2', 'pi'])
ax_phase.axhline(0, color='#cccccc', linewidth=0.8, zorder=0)
ax_phase.set_xlim(freqs[0], 100)
ax_phase.legend(loc='upper right', fontsize=9)
fig.tight_layout()
plt.show()
k_slider = widgets.IntSlider(
value=5,
min=0,
max=30,
step=1,
description='k index:',
style={'description_width': 'initial'},
layout=widgets.Layout(width='420px'),
)
output = widgets.interactive_output(draw, {'k': k_slider})
ui = widgets.VBox([k_slider, output])
display(ui)
visualize_single_k_dft_demo()
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Exercise_1_DFT")
Congratulations!#
You’ve just implemented the DFT from scratch. You can now compute the Fourier coefficients for any signal, and visualize its amplitude / power spectrum and phase spectrum.
In practice, you will almost always use the optimized numpy.fft.fft or scipy.fft.fft functions, which are much faster than the naive dot-product implementation we wrote.
Physical frequencies#
There is one more thing: the DFT coefficients are indexed by integers \(k=0, 1, \dots, N-1\), where \(N\) is the number of samples in the signal. But they also correspond to actual, physical frequencies in Hz, like we plot in the amplitude and phase spectra. The mapping from \(k\) to frequency is: $\(f_k = \frac{k}{N} f_s\)\( where \)f_k\( is the frequency in Hz and \)f_s$ is the sampling frequency of the signal.
This means that the DFT coefficients are:
evenly spaced integer multiples of the minimum frequency \(f_1 = \frac{f_s}{N}\),
the minimum (slowest) frequency we can resolve is \(f_1 = \frac{f_s}{N}\) (the inverse of the total duration of the signal in seconds), which is the slowest sine wave that completes one full cycle over the duration of the signal
and the maximum (fastest) frequency we can resolve is \(f_{N/2} = f_s / 2\), which is the Nyquist frequency and the fastest sine wave that can be accurately represented between 2 samples per cycle.
np.fft.fftfreq is a convenient function that computes the physical frequencies given the number of samples and the sampling frequency.
Section 5 : Fourier Analysis of Neural Data#
Estimated timing to here from start of tutorial: 40 min
Video 4: Real data and wrap-up#
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Video_4_WrapUp")
Load IBL Neuropixels Data — Streams directly from DANDI (run me first!)#
We finish this tutorial by looking at a real neural signal and its spectrum.
We will download a snippet of data from the International Brain Laboratory (IBL) dataset. Streams 2 seconds from subject DY-008, channel 60 (best LFP+spike channel). No large download needed — remfile fetches only the bytes we ask for.
url = 'https://api.dandiarchive.org/api/assets/604b7457-430e-4160-b7ff-594562d942e3/download/'
f = h5py.File(remfile.File(url), 'r')
ap = f['acquisition']['ElectricalSeriesProbe00AP']
ap_fs = 30000.0 # Hz, confirmed from timestamps
conv = 2.34375e-06 # to volts
BEST_CH = 60 # channel with clearest LFP structure + spikes
T_DUR = 2.0 # seconds to stream
i1 = int(T_DUR * ap_fs)
t_ap = np.arange(i1) / ap_fs
print(f'Streaming {T_DUR}s from channel {BEST_CH} at {ap_fs} Hz...')
# Contiguous slice required by remfile — [:, CH:CH+1] not [:, CH]
raw_ap = ap['data'][:i1, BEST_CH:BEST_CH+1].astype(float)[:, 0] * conv * 1e6 # µV
raw_ap = raw_ap - np.mean(raw_ap) # zero-centre
print(f'Done. Shape: {raw_ap.shape}, duration: {t_ap[-1]:.2f}s')
print(f'Signal range: {raw_ap.min():.1f} to {raw_ap.max():.1f} µV')
Streaming 2.0s from channel 60 at 30000.0 Hz...
Done. Shape: (60000,), duration: 2.00s
Signal range: -344.4 to 185.3 µV
Spectral analysis allows us to see which frequencies are present in the data, and how strong they are. This is particularly useful in neurophysiology to identify action potentials, brain rhythms, or noise components.
Here’s a basic overview of the process:
Time Domain Signal: We start with a signal recorded over time (e.g.,
raw_apin our notebook).Fast Fourier Transform (FFT): We apply the FFT to this time-domain signal (
scipy.fft.fft). The output of the FFT is a complex-valued array (which you implemented in the previous exercise), where each element corresponds to a specific frequency. It tells us the amplitude and phase of each frequency component present in the original signal.Frequency Axis (
xf): Thefftfreqfunction (fromscipy.fft) is used to generate the corresponding frequencies for each point in the FFT output. Otherwise, the components would only correspond to indices, not actual frequencies with physical meaning.Spectrum Calculation: From the FFT output, we can derive different types of spectra:
Amplitude Spectrum: This shows the magnitude (amplitude) of each frequency component. It’s calculated as
2.0/N * np.abs(yf[0:N//2]), whereNis the number of samples andyfis the FFT output. The2.0/Nscaling normalizes the amplitude, and[0:N//2]considers only the positive frequency components (the FFT output is symmetric for real-valued signals).Power Spectrum: This shows the power of each frequency component. It’s often calculated as the square of the amplitude spectrum, sometimes with an additional normalization (
(1.0/N * np.abs(yf[0:N//2]))**2). Power spectra are useful because power is additive and directly related to the energy content at each frequency. This is often preferred when comparing the “strength” of different frequency bands.We can also compute the phase spectrum, which we omit here.
By plotting these spectra, we can visually identify dominant frequencies, broadband noise, or other spectral features in the signal. We can further plot the frequency (x-axis) and/or amplitude/power (y-axis) on a logarithmic scale to better visualize components across a wider range.
# Perform FFT on the raw signal
N = len(raw_ap) # Number of sample points
T = 1.0 / ap_fs # Sample spacing (1/sampling frequency)
yf = fft(raw_ap)
xf = fftfreq(N, T)[:N//2]
amp_spectrum = 2.0/N * np.abs(yf[0:N//2])
power_spectrum = (1.0/N * np.abs(yf[0:N//2]))**2
def combined_interactive_plot(spectrum_type, x_scale, y_scale):
fig, axes = plt.subplots(2, 1, figsize=(13, 10), sharex=False) # sharex=False as x-axes are different (time vs freq)
# Top Panel: Raw Signal Over Time
axes[0].plot(t_ap, raw_ap, color='C4', lw=0.6)
axes[0].set_title('Raw Broadband Signal Over Time')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Amplitude (µV)')
axes[0].grid(True, linestyle='--', alpha=0.7)
# Bottom Panel: Spectrum Plot
plot_xf = xf
plot_y_data = None
y_label = ''
if spectrum_type == 'Amplitude Spectrum':
plot_y_data = amp_spectrum
y_label = 'Amplitude (µV)'
else: # Power Spectrum
# Using the convention from previous notebook cells for consistency
plot_y_data = power_spectrum
y_label = 'Power (µV$^2$)'
# Exclude DC component (0 Hz) if x-axis is logarithmic to avoid issues with log(0)
if x_scale == 'log':
# Find the first non-zero frequency to plot
first_nonzero_freq_idx = np.where(plot_xf > 0)[0]
if len(first_nonzero_freq_idx) > 0:
start_idx = first_nonzero_freq_idx[0]
plot_xf_log = plot_xf[start_idx:]
plot_y_data_log = plot_y_data[start_idx:]
else:
# Fallback for extremely unusual cases, or if xf is all zeros (shouldn't happen with fftfreq)
plot_xf_log = np.array([1e-6]) # A small positive value to allow log plot
plot_y_data_log = np.array([0])
axes[1].plot(plot_xf_log, plot_y_data_log, color='C1')
else:
axes[1].plot(plot_xf, plot_y_data, color='C1')
axes[1].set_title(f'{spectrum_type} of Raw Broadband Signal')
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel(y_label)
axes[1].set_xscale(x_scale)
axes[1].set_yscale(y_scale)
axes[1].grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
# Create widgets for user control
spectrum_type_widget = widgets.RadioButtons(
options=['Amplitude Spectrum', 'Power Spectrum'],
value='Amplitude Spectrum',
description='Spectrum Type:',
disabled=False,
layout=widgets.Layout(width='auto')
)
x_scale_widget = widgets.RadioButtons(
options=['linear', 'log'],
value='linear',
description='X-axis Scale: ',
disabled=False,
layout=widgets.Layout(width='auto')
)
y_scale_widget = widgets.RadioButtons(
options=['linear', 'log'],
value='linear',
description='Y-axis Scale: ',
disabled=False,
layout=widgets.Layout(width='auto')
)
# The interactive_output itself is a widget that captures the plot
output_plot = widgets.interactive_output(combined_interactive_plot, {
'spectrum_type': spectrum_type_widget,
'x_scale': x_scale_widget,
'y_scale': y_scale_widget
})
def combined_interactive_plot(spectrum_type, x_scale, y_scale):
fig, axes = plt.subplots(2, 1, figsize=(13, 10), sharex=False) # sharex=False as x-axes are different (time vs freq)
# Top Panel: Raw Signal Over Time
axes[0].plot(t_ap, raw_ap, color='C4', lw=0.6)
axes[0].set_title('Raw Broadband Signal Over Time')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Amplitude (µV)')
axes[0].grid(True, linestyle='--', alpha=0.7)
# Bottom Panel: Spectrum Plot
plot_xf = xf
plot_y_data = None
y_label = ''
if spectrum_type == 'Amplitude Spectrum':
plot_y_data = amp_spectrum
y_label = 'Amplitude (µV)'
else: # Power Spectrum
# Using the convention from previous notebook cells for consistency
plot_y_data = power_spectrum
y_label = 'Power (µV$^2$)'
# Exclude DC component (0 Hz) if x-axis is logarithmic to avoid issues with log(0)
if x_scale == 'log':
# Find the first non-zero frequency to plot
first_nonzero_freq_idx = np.where(plot_xf > 0)[0]
if len(first_nonzero_freq_idx) > 0:
start_idx = first_nonzero_freq_idx[0]
plot_xf_log = plot_xf[start_idx:]
plot_y_data_log = plot_y_data[start_idx:]
else:
# Fallback for extremely unusual cases, or if xf is all zeros (shouldn't happen with fftfreq)
plot_xf_log = np.array([1e-6]) # A small positive value to allow log plot
plot_y_data_log = np.array([0])
axes[1].plot(plot_xf_log, plot_y_data_log, color='C1')
else:
axes[1].plot(plot_xf, plot_y_data, color='C1')
axes[1].set_title(f'{spectrum_type} of Raw Broadband Signal')
axes[1].set_xlabel('Frequency (Hz)')
axes[1].set_ylabel(y_label)
axes[1].set_xscale(x_scale)
axes[1].set_yscale(y_scale)
axes[1].grid(True, linestyle='--', alpha=0.7)
plt.tight_layout()
plt.show()
# Create widgets for user control
spectrum_type_widget = widgets.RadioButtons(
options=['Amplitude Spectrum', 'Power Spectrum'],
value='Amplitude Spectrum',
description='Spectrum Type:',
disabled=False,
layout=widgets.Layout(width='auto')
)
x_scale_widget = widgets.RadioButtons(
options=['linear', 'log'],
value='linear',
description='X-axis Scale: ',
disabled=False,
layout=widgets.Layout(width='auto')
)
y_scale_widget = widgets.RadioButtons(
options=['linear', 'log'],
value='linear',
description='Y-axis Scale: ',
disabled=False,
layout=widgets.Layout(width='auto')
)
# The interactive_output itself is a widget that captures the plot
output_plot = widgets.interactive_output(combined_interactive_plot, {
'spectrum_type': spectrum_type_widget,
'x_scale': x_scale_widget,
'y_scale': y_scale_widget
})
# @title Arrange the widgets and the plot output in a VBox
display(widgets.VBox([
output_plot, # The actual plot output from interactive_output
widgets.HBox([spectrum_type_widget, x_scale_widget, y_scale_widget]) # Controls second
]))
Think! 2: What do you see in the spectrum?#
Notice some prominent features in the neural signal’s spectrum:
there are some very sharp peaks at specific frequencies,
and a broad hump from ~200 to 4000Hz.
What do you think these features correspond to?
Summary#
Estimated timing of tutorial: 45 min
In this tutorial you built up the Fourier picture of signals from a couple of sine waves. You saw that:
any signals can be decomposed into a weighted sum of sines
the magnitude spectrum is just those weights indexed by frequency
the phase of each component is crucial for determining the time-domain shape of the signal
and that we can find the magnitude and phase coefficients through discrete Fourier transform (DFT), which you implemented yourself.
With these tools, you can decompose any signal in a computer into its frequency components, analyze its spectral content for contributions from different frequencies, and understand how phase affects the time-domain representation.
Almost any neural signal you work with — be it spikes, LFPs, EEG, or fMRI — will have been preprocessed or analyzed using Fourier techniques.
But be careful!#
Just because you can decompose a signal into sinusoidal components, it does not mean that the signal is actually made of sine waves in a physical or mechanistic sense! This is particularly true for putative neural oscillations.
Submit your feedback on this tutorial#
Show code cell source
# @title Submit your feedback on this tutorial
content_review(f"{feedback_prefix}_Overall")
Bonus: Coordinate System and Change of Basis#
In Demo 1 above, we essentially changed the representation of our mystery signal from a time-domain basis, i.e., the signal at every point as a function of time (top) to a frequency-domain basis, i.e., the signal as a function of frequency, specifically the sum of three sine waves at frequencies 0Hz, 3Hz, and 7Hz (bottom right).
This is an example of a change of basis. We will now build up to the “Fourier basis”, starting from a review of the concept of basis in 2D, which you learned about for Principal Component Analysis (PCA).
Basis in 2D#
Definition:
A set of vectors {\(\mathbf{e_1},\mathbf{e_2}\)} is a basis for \(\mathbb{R}^2\) if:
They are linearly independent, i.e., if the equation:
\[c_1 \mathbf{e_1} + c_2 \mathbf{e_2} = 0\]is only satisfied when both \(c_1\) and \(c_2\) are zero.
In \(\mathbb{R}^2\) this intuitively means that the two vectors that form a basis cannot have the same direction.
Every vector \(\mathbf{v}\) in \(\mathbb{R}^2\) can be written as their linear combination: \(\mathbf{v} = d_1\mathbf{e_1} + d_2\mathbf{e_2}\) for some \(d_1\), \(d_2\), which are called the coordinates in this basis.
Furthermore, if the dot product of basis vectors is zero (\(\mathbf{e_1} \cdot \mathbf{e_2} = 0\)) such a basis is called orthogonal.
Lastly, if the vectors \(\mathbf{e_1}\) and \(\mathbf{e_2}\) that form an orthogonal basis also have unit length (\(\mathbf{e_1} ⋅ \mathbf{e_1} = \mathbf{e_2} ⋅\mathbf{e_2} = 1\)), such a basis is then called orthonormal.
Cartesian (rectangular) coordinate system#
A coordinate system on \(\mathbb{R}^2\) consists of:
An origin - reference point (0,0)
A basis
Coordinates
So any point can be uniquely represented by its coordinate vector \(\mathbf{v}\): $\(\mathbf{v} = d_1 \mathbf{e_1} + d_2 \mathbf{e_2}\)$
The most common coordinate system is Cartesian (Standard):
Basis: \(\mathbf{e_1} = \begin{bmatrix} 1\\0 \end{bmatrix}\) - horizontal or “x” direction, and \(\mathbf{e_2} = \begin{bmatrix} 0\\1 \end{bmatrix}\) - vertical or “y” direction.
Any point can then be expressed by specifying numbers x and y such that:
But the Standard coordinate system is not the only one.
Polar coordinate system#
Points are specified by the radius \( r = \sqrt{x^2 + y^2}\) from the origin and angle \(ϕ = \arctan{\frac{y}{x}}\) from x-axis.
Basis: \(\mathbf{e_1} = \begin{bmatrix} \cos{ϕ}\\ \sin{ϕ} \end{bmatrix} = \begin{bmatrix} \frac{x}{r}\\ \frac{y}{r} \end{bmatrix}\), \(\mathbf{e_2} = \begin{bmatrix} -\sin{ϕ}\\ \cos{ϕ} \end{bmatrix} = \begin{bmatrix} -\frac{y}{r}\\ \frac{x}{r} \end{bmatrix}\)
Any point can then be expressed by specifying numbers \(r\) and \(ϕ\) such that \(\mathbf{v} = d_r \mathbf{e_1} + d_ϕ \mathbf{e_2}\)
Think! Choice of basis#
There are many different basis one could use to achieve one’s goal:
What are the advantages of orthogonal basis and orthonormal basis compared to the arbitrary one?
Are vectors \(\mathbf{e_1} = \begin{bmatrix} 1\\0 \end{bmatrix}\), \(\mathbf{e_2} = \begin{bmatrix} 2\\1 \end{bmatrix}\) forming a basis? Is it orthogonal? Orthonormal?
What about \(\mathbf{e_1} = \begin{bmatrix} sin(\pi/4)\\ cos(\pi/4) \end{bmatrix}\), \(\mathbf{e_2} = \begin{bmatrix} -sin(\pi/4)\\cos(\pi/4) \end{bmatrix}\)?
Functions as Basis#
Similar to 2D space, in higher dimensions, any vector can be represented as a linear combination of basis vectors. Thus if one finds an orthogonal basis {\(\mathbf{e_1},\mathbf{e_2}, …, \mathbf{e_N}\)} in \(N\) dimensional space, it would be possible to construct any vector \(\mathbf{v}\) as a linear combination of this basis:
So far we have only discussed vectors as ordered lists of numbers (arrays) and its geometrical representation as arrows but we could move to another level of abstraction.
It turns out that functions also form a vector space (function space), thus if one has an orthogonal set of functions {\({f_1},{f_2}, …, {f_N}\)}, it is possible to represent a function \(g\) from the same space as a linear combination of basis functions: $\({g} = a_1{f_1} + a_2{f_2} + a_3{f_3} + \dots + a_N{f_N}\)$
In function space, orthogonality is defined via the inner product \(\left<f_1,f_2\right>\). We will define the inner product as follows:
Similarly to the dot product, if the inner product of two different functions is zero, they are orthogonal.
Sinusoidal basis for signals#
Concretely, sine waves can act as a basis for signals:
In the previous demo, we already saw that the mystery signal can be represented as a linear combination of sine waves. This wasn’t surprising, since we told you that the mystery signal is composed of sine waves and gave you the frequencies.
But as it turns out, any signal can be represented as a sum of sine and cosine waves, and that certain sets of sine waves are orthogonal to each other (as we will see below), which makes them a convenient choice for a basis.
Bonus Interactive Demo: Orthogonal Sinusoids#
This demo shows that two sine waves with different frequencies are orthogonal to each other on the time interval \([0,T]\), where \(T = 1\).
Specifically, we look at “trigonometric” functions: \(\sin(n*ωt + ϕ)\) and \(\cos(n*ωt + ϕ)\), where \(ω = \frac{2π}{T}\) is the fundamental frequency, and \(m\) and \(n\) are integers that index different “vectors” in this basis. Note that the bigger the value of \(m\) and \(n\), the faster the oscillation.
In this figure, the plots are:
Top: two sinusoidal functions where you control their frequencies and phase shifts
Bottom-left: the element-wise multiplication of the two functions at each point in time, i.e., \(f_1(t)f_2(t)\)
Bottom-right: the cumulative integral of the product \(\int_0^t f_1(\tau)f_2(\tau)d\tau\), i.e., the running sum of the left plot, and the full integral value over the period is displayed in the legend box.
Start with \(n\) and \(m\): Find values for those parameters where the two sinusoidal functions are orthogonal to each other, i.e, the inner product or integral is zero.
What difference does the phase shift make?
Make sure you execute this cell to enable the widget!#
Show code cell source
# @title Make sure you execute this cell to enable the widget!
def get_function_and_label(mode, n, m, phi1, phi2, T):
"""Helper function to get the waves and the labels
Inputs:
mode: "sine-sine", "cosine-cosine" or "sine-cosine"
n,m: integer numbers for frequencies
phi1, phi2: phase shift
T: period
Returns:
f1, f2: functions with respect to chosen mode and arguments n,m,phi1,phi2
label1, label2: labels for the functions
"""
omega = 2 * np.pi / T
if mode == "sine-sine":
def f1(t):
return np.sin(omega * n * t + phi1 * np.pi)
def f2(t):
return np.sin(omega * m * t + phi2 * np.pi)
label1 = rf"$\sin({int(n)}\omega t + {phi1}\pi)$"
label2 = rf"$\sin({int(m)}\omega t + {phi2}\pi)$"
elif mode == "sine-cosine":
def f1(t):
return np.sin(omega * n * t + phi1 * np.pi)
def f2(t):
return np.cos(omega * m * t + phi2 * np.pi)
label1 = rf"$\sin({n}\omega t + {phi1}\pi)$"
label2 = rf"$\cos({m}\omega t + {phi2}\pi)$"
elif mode == "cosine-cosine":
def f1(t):
return np.cos(omega * n * t + phi1 * np.pi)
def f2(t):
return np.cos(omega * m * t + phi2 * np.pi)
label1 = rf"$\cos({n}\omega t + {phi1}\pi)$"
label2 = rf"$\cos({m}\omega t + {phi2}\pi)$"
return f1, f2, label1, label2
def plot_all(n=2, m=1, phi1=0, phi2=0, mode="sine-sine"):
T = 1
f1, f2, label1, label2 = get_function_and_label(mode, n, m, phi1, phi2, T)
t = np.linspace(0, T, 1000)
integrand = lambda t: f1(t) * f2(t)
integral_scalar = lambda tt: quad(integrand, t[0], tt)[0]
integral = np.vectorize(integral_scalar)
# arrange the plots
fig = plt.figure(figsize=(9, 6))
gs = fig.add_gridspec(2, 2, height_ratios=[1, 1])
ax1 = fig.add_subplot(gs[0, :]) # top row
ax2 = fig.add_subplot(gs[1, 0]) # bottom left
ax3 = fig.add_subplot(gs[1, 1]) # bottom right
# Row 1: sine/cosine waves
ax1.plot(t, f1(t), label=label1)
ax1.plot(t, f2(t), label=label2)
ax1.legend(loc="upper right")
ax1.set_xlabel("t")
ax1.set_ylabel("Value")
ax1.set_title("Waves")
ax1.axhline(0, color="black", linewidth=1)
ax1.grid(alpha=0.25)
ax1.set_xlim(t[0], t[-1])
# Row 2 left: integrand
label_integrand = label1 + r"$\cdot$" + label2
y = integrand(t)
ax2.plot(t, y)
ax2.fill_between(t, y, 0, where=(y >= 0), color="#99BA5A", alpha=0.8)
ax2.fill_between(t, y, 0, where=(y < 0), color="#C9B541", alpha=0.8)
ax2.set_xlabel("t")
ax2.set_ylabel("Value")
ax2.set_title(label_integrand, fontsize=12)
ax2.axhline(0, color="black", linewidth=1)
ax2.grid(alpha=0.25)
ax2.set_xlim(t[0], t[-1])
# Row 2 right: cumulative integral
label_integral = rf"$\int_{{0}}^{{t}}$" + label_integrand + r"$dt$"
ax3.plot(t, integral(t))
ax3.set_xlabel("t")
ax3.set_ylabel("Value")
ax3.set_title(label_integral, fontsize=12)
ax3.axhline(0, color="black", linewidth=1)
ax3.grid(alpha=0.25)
ax3.set_xlim(t[0], t[-1])
ax3.text(
0.05, 0.95,
f"Full integral: {integral(t[-1]):.3f}",
transform=ax3.transAxes,
color="black",
fontsize=11,
verticalalignment="top",
bbox=dict(facecolor="white", alpha=0.8, edgecolor="black")
)
fig.tight_layout()
plt.show()
### Widgets and Display ###
n_widget = widgets.FloatSlider(value=2, min=-10, max=10, step=1, description='n')
m_widget = widgets.FloatSlider(value=1, min=-10, max=10, step=1, description='m')
phi1_widget = widgets.FloatSlider(value=0, min=0, max=2, step=0.05, description='phi1')
phi2_widget = widgets.FloatSlider(value=0, min=0, max=2, step=0.05, description='phi2')
mode_widget = widgets.Dropdown(options=["sine-sine", "cosine-cosine", "sine-cosine"], value="sine-sine", description='mode')
controls = widgets.VBox([mode_widget, n_widget, m_widget, phi1_widget, phi2_widget])
output = widgets.interactive_output(
plot_all,
{
'n': n_widget,
'm': m_widget,
'phi1': phi1_widget,
'phi2': phi2_widget,
'mode': mode_widget
}
)
ui = widgets.VBox([controls, output])
display(ui)
Submit your feedback#
Show code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Demo_Bonus_FourierBasis")