00 — Introduction to Dynamic Mode Decomposition¶

This notebook introduces Dynamic Mode Decomposition (DMD) and shows how it can decompose periodic motion into interpretable oscillatory modes.

We demonstrate the method on a synthetic example, then walk through the fields returned by run_dmd() in the DMDResult dataclass.

Mathematical background¶

DMD approximates a time series $\mathbf{x}(t)$ as a sum of modes:

$$ \mathbf{x}(t) \;=\; \sum_{k=1}^{r} b_k\,\boldsymbol{\varphi}_k\,e^{\omega_k\,t} $$

where $\boldsymbol{\varphi}_k$ are spatial modes, $\omega_k$ are complex eigenvalues encoding frequency and growth/decay, and $b_k$ are amplitudes.

In [1]:

%load_ext autoreload
%autoreload 2
%matplotlib widget
%config InlineBackend.figure_format='svg'

import matplotlib.pyplot as plt
import numpy as np

from birddmd import compute_rmse, reconstruct, run_dmd

%load_ext autoreload %autoreload 2 %matplotlib widget %config InlineBackend.figure_format='svg' import matplotlib.pyplot as plt import numpy as np from birddmd import compute_rmse, reconstruct, run_dmd

1 — Synthetic example: sum of sinusoids¶

We create a 2D signal that is a mixture of two frequencies, then use DMD to recover them.

In [2]:

# Create synthetic data: two frequencies in three spatial channels (1 marker)
dt = 0.01
t_synth = np.arange(0, 2, dt)
f1, f2 = 3.0, 7.5  # Hz

x1 = 2.0 * np.sin(2 * np.pi * f1 * t_synth) + 0.5 * np.cos(2 * np.pi * f2 * t_synth)
x2 = 1.0 * np.cos(2 * np.pi * f1 * t_synth) + 1.5 * np.sin(2 * np.pi * f2 * t_synth)
x3 = 0.8 * np.sin(2 * np.pi * f1 * t_synth) - 0.3 * np.cos(2 * np.pi * f2 * t_synth)

data_synth = np.column_stack([x1, x2, x3])  # shape (200, 3) -- 1 marker x 3 axes

fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(t_synth, x1, label="Channel 1")
ax.plot(t_synth, x2, label="Channel 2")
ax.plot(t_synth, x3, label="Channel 3")
ax.set_xlabel("Time (s)")
ax.set_ylabel("Amplitude")
ax.legend()
ax.set_title("Synthetic 3-channel signal")
plt.tight_layout()
plt.show()

# Create synthetic data: two frequencies in three spatial channels (1 marker) dt = 0.01 t_synth = np.arange(0, 2, dt) f1, f2 = 3.0, 7.5 # Hz x1 = 2.0 * np.sin(2 * np.pi * f1 * t_synth) + 0.5 * np.cos(2 * np.pi * f2 * t_synth) x2 = 1.0 * np.cos(2 * np.pi * f1 * t_synth) + 1.5 * np.sin(2 * np.pi * f2 * t_synth) x3 = 0.8 * np.sin(2 * np.pi * f1 * t_synth) - 0.3 * np.cos(2 * np.pi * f2 * t_synth) data_synth = np.column_stack([x1, x2, x3]) # shape (200, 3) -- 1 marker x 3 axes fig, ax = plt.subplots(figsize=(8, 3)) ax.plot(t_synth, x1, label="Channel 1") ax.plot(t_synth, x2, label="Channel 2") ax.plot(t_synth, x3, label="Channel 3") ax.set_xlabel("Time (s)") ax.set_ylabel("Amplitude") ax.legend() ax.set_title("Synthetic 3-channel signal") plt.tight_layout() plt.show()

Figure

In [3]:

# Run DMD on the synthetic data
synth_avg = np.zeros(3)  # zero mean for pre-centred synthetic data

result_synth = run_dmd(
    data=data_synth,  # (n_timesteps, 3)
    times=t_synth,
    n_modes=4,
    d=2,
    eig_constraints={"conjugate_pairs"},
    average_shape=synth_avg,
    n_markers=1,
    verbose=True,
)

# Recovered frequencies
freqs_hz = np.imag(result_synth.eigenvalues) / (2 * np.pi)
print(f"\nRecovered frequencies: {np.round(np.abs(freqs_hz), 2)} Hz")
print(f"True frequencies:      [{f1}, {f2}] Hz")

# Run DMD on the synthetic data synth_avg = np.zeros(3) # zero mean for pre-centred synthetic data result_synth = run_dmd( data=data_synth, # (n_timesteps, 3) times=t_synth, n_modes=4, d=2, eig_constraints={"conjugate_pairs"}, average_shape=synth_avg, n_markers=1, verbose=True, ) # Recovered frequencies freqs_hz = np.imag(result_synth.eigenvalues) / (2 * np.pi) print(f"\nRecovered frequencies: {np.round(np.abs(freqs_hz), 2)} Hz") print(f"True frequencies: [{f1}, {f2}] Hz")

Running DMD with 4 modes, delay d=2
Input shape: (200, 3)
Number of variables in DMD results: 6
Complex Eigenvalues (log): [ 0.+18.85j   0.-18.85j  -0.+47.124j -0.-47.124j]
Number of modes in Psi: 4

Recovered frequencies: [3.  3.  7.5 7.5] Hz
True frequencies:      [3.0, 7.5] Hz

/Users/lfrance/Library/CloudStorage/OneDrive-Personal/004 GitHub/BirdDMD/.venv/lib/python3.12/site-packages/pydmd/snapshots.py:73: UserWarning: Input data condition number 2037172144721110.0. Consider preprocessing data, passing in augmented data
matrix, or regularization methods.
  warnings.warn(

2 — Understanding DMDResult¶

run_dmd() returns a frozen DMDResult dataclass containing every output of the analysis. All fields are accessible as attributes.

In [4]:

print(f"DMDResult fields: {list(result_synth.__dataclass_fields__)}\n")
print(
    f"  eigenvalues:      {result_synth.eigenvalues.shape}"
    "  complex, lambda = sigma + i*omega"
)
print(f"  modes:            {result_synth.modes.shape}  spatial DMD modes phi")
print(f"  amplitudes:       {result_synth.amplitudes.shape}  mode weights b")
print(f"  frequencies_hz:   {result_synth.frequencies_hz.shape}  Im(lambda)/2pi")
print(f"  growth_rates:     {result_synth.growth_rates.shape}  Re(lambda)")
print(f"  mode_magnitudes:  {result_synth.mode_magnitudes.shape}  per-marker |phi|")
print(f"  phase_shifts:     {result_synth.phase_shifts.shape}  phase angles theta")
print(f"  conjugate_pairs:  {result_synth.conjugate_pairs}")
print(f"  reconstruction:   {result_synth.reconstruction.shape}  full reconstruction")

print(f"DMDResult fields: {list(result_synth.__dataclass_fields__)}\n") print( f" eigenvalues: {result_synth.eigenvalues.shape}" " complex, lambda = sigma + i*omega" ) print(f" modes: {result_synth.modes.shape} spatial DMD modes phi") print(f" amplitudes: {result_synth.amplitudes.shape} mode weights b") print(f" frequencies_hz: {result_synth.frequencies_hz.shape} Im(lambda)/2pi") print(f" growth_rates: {result_synth.growth_rates.shape} Re(lambda)") print(f" mode_magnitudes: {result_synth.mode_magnitudes.shape} per-marker |phi|") print(f" phase_shifts: {result_synth.phase_shifts.shape} phase angles theta") print(f" conjugate_pairs: {result_synth.conjugate_pairs}") print(f" reconstruction: {result_synth.reconstruction.shape} full reconstruction")

DMDResult fields: ['eigenvalues', 'modes', 'amplitudes', 'frequencies_hz', 'growth_rates', 'mode_magnitudes', 'phase_shifts', 'conjugate_pairs', 'reconstruction', 'average_shape', 'times', 'sort_indices', '_dmd_object']

  eigenvalues:      (4,)  complex, lambda = sigma + i*omega
  modes:            (3, 4)  spatial DMD modes phi
  amplitudes:       (4,)  mode weights b
  frequencies_hz:   (4,)  Im(lambda)/2pi
  growth_rates:     (4,)  Re(lambda)
  mode_magnitudes:  (1, 3, 4)  per-marker |phi|
  phase_shifts:     (3, 4)  phase angles theta
  conjugate_pairs:  [(0, 1), (2, 3)]
  reconstruction:   (199, 1, 3)  full reconstruction

Field reference¶

Field	Shape	Math	Description
eigenvalues	(n_modes,)	$\lambda = \sigma + i\omega$	Complex eigenvalues
modes	(n_coords, n_modes)	$\varphi$	Spatial DMD modes
amplitudes	(n_modes,)	$b$	Mode weights
frequencies_hz	(n_modes,)	$\text{Im}(\lambda)/2\pi$	Oscillation frequency
growth_rates	(n_modes,)	$\text{Re}(\lambda)$	Exponential growth/decay
mode_magnitudes	(n_markers, 3, n_modes)	$\lvert\varphi\rvert$	Per-marker magnitude
phase_shifts	(n_markers, 3, n_modes)	$\theta$	Phase angles
conjugate_pairs	list of (i, j)	—	Paired eigenvalue indices
reconstruction	(T-1, n_markers, 3)	$\tilde{\mathbf{x}}(t)$	Full reconstruction

Helper properties: n_modes, n_pairs, pair_frequency(i), pair_indices([0, 1]).

Annotated walkthrough¶

Using the synthetic example above, we step through the key fields and helper methods. Eigenvalues encode frequency and growth/decay; conjugate pairs group complex-conjugate eigenvalue pairs; and reconstruct() can target specific mode subsets.

In [5]:

# Eigenvalues encode frequency and growth/decay
for k, eig in enumerate(result_synth.eigenvalues):
    freq = np.abs(np.imag(eig)) / (2 * np.pi)
    growth = np.real(eig)
    print(f"Mode {k}: {freq:.1f} Hz, growth = {growth:.4f}")

# Conjugate pairs group complex-conjugate eigenvalues
print(f"\nConjugate pairs: {result_synth.conjugate_pairs}")
print(f"Number of pairs: {result_synth.n_pairs}")

# Pair frequency helper
for p in range(result_synth.n_pairs):
    print(f"  Pair {p}: {result_synth.pair_frequency(p):.1f} Hz")

# Reconstruct specific mode subsets
kp_pair0 = reconstruct(result_synth, pairs=[0])
print(f"\nPair 0 reconstruction: {kp_pair0.shape}")

# RMSE
rmse = compute_rmse(result_synth.reconstruction, data_synth[:-1].reshape(-1, 1, 3))
print(f"Full reconstruction RMSE: {np.mean(rmse):.6f}")

# Eigenvalues encode frequency and growth/decay for k, eig in enumerate(result_synth.eigenvalues): freq = np.abs(np.imag(eig)) / (2 * np.pi) growth = np.real(eig) print(f"Mode {k}: {freq:.1f} Hz, growth = {growth:.4f}") # Conjugate pairs group complex-conjugate eigenvalues print(f"\nConjugate pairs: {result_synth.conjugate_pairs}") print(f"Number of pairs: {result_synth.n_pairs}") # Pair frequency helper for p in range(result_synth.n_pairs): print(f" Pair {p}: {result_synth.pair_frequency(p):.1f} Hz") # Reconstruct specific mode subsets kp_pair0 = reconstruct(result_synth, pairs=[0]) print(f"\nPair 0 reconstruction: {kp_pair0.shape}") # RMSE rmse = compute_rmse(result_synth.reconstruction, data_synth[:-1].reshape(-1, 1, 3)) print(f"Full reconstruction RMSE: {np.mean(rmse):.6f}")

Mode 0: 3.0 Hz, growth = 0.0000
Mode 1: 3.0 Hz, growth = 0.0000
Mode 2: 7.5 Hz, growth = -0.0000
Mode 3: 7.5 Hz, growth = -0.0000

Conjugate pairs: [(0, 1), (2, 3)]
Number of pairs: 2
  Pair 0: 3.0 Hz
  Pair 1: 7.5 Hz

Pair 0 reconstruction: (199, 1, 3)
Full reconstruction RMSE: 0.000000

Next¶

The next notebook (01_flapping_modes) applies DMD to real hawk flight data — motion-capture recordings of a Harris's hawk in straight flapping flight — and decomposes the wingbeat into three interpretable conjugate mode pairs.