The Python lab section runs entirely in your browser using Pyodide(Python 3.11 compiled to WebAssembly). NumPy, SciPy and Matplotlib are available out of the box; you can modify the code, re-run it, and see plots immediately. The status indicator below shows when the runtime is ready — usually a few seconds on first load, then cached for subsequent visits.

Loading Pyodide + NumPy + SciPy + Matplotlib…

Experiment 1 — Acoustic Complexity Index (ACI)Signal Processing

Simulate a spectrogram with biotic and abiotic components, compute ACI, and visualise how index values respond to signal complexity. Modify n_species or noise_level to explore.

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import io, base64

np.random.seed(42)
n_species, noise_level = 8, 0.15
n_freqs, n_times = 64, 200

S = np.full((n_freqs, n_times), noise_level) + \
    np.random.exponential(noise_level * 0.5, (n_freqs, n_times))
for _ in range(n_species):
    f0 = np.random.randint(10, n_freqs - 10)
    bw = np.random.randint(2, 8)
    t0 = np.random.randint(0, n_times - 40)
    dur = np.random.randint(15, 50)
    amp = np.random.uniform(0.5, 1.5)
    freq_env = np.exp(-0.5 * ((np.arange(n_freqs) - f0) / bw)**2)
    time_env = np.zeros(n_times)
    time_env[t0:t0+dur] = amp * (0.5 + 0.5 * np.sin(
        np.linspace(0, np.pi * np.random.randint(3, 10), dur)))
    S += np.outer(freq_env, time_env)

def compute_ACI(S):
    total = 0.0
    for k in range(S.shape[0]):
        row = S[k]; diffs = np.abs(np.diff(row))
        total += np.sum(diffs) / (np.sum(row) + 1e-10)
    return total

ACI_val = compute_ACI(S)
S_null  = np.random.exponential(0.3, (n_freqs, n_times))
ACI_null = compute_ACI(S_null)

print(f"ACI (biotic simulation) : {ACI_val:.3f}")
print(f"ACI (white noise null)  : {ACI_null:.3f}")
print(f"Ratio (biotic/null)     : {ACI_val/ACI_null:.3f}")

fig, axes = plt.subplots(1, 3, figsize=(13, 4), facecolor='#0a1a0f',
    gridspec_kw={'width_ratios':[3,3,1]})
fig.suptitle('Acoustic Complexity Index - Simulation', color='#c8e0cc',
             fontsize=11, fontfamily='monospace')
for ax in axes:
    ax.set_facecolor('#040d07')
    for sp in ax.spines.values(): sp.set_color('#2d6a4f')
    ax.tick_params(colors='#6a9a7a', labelsize=8)
axes[0].imshow(S, aspect='auto', origin='lower', cmap='YlGn',
               extent=[0, n_times, 0, n_freqs])
axes[0].set_title('Biotic Spectrogram', color='#52b788', fontsize=9)
axes[1].imshow(S_null, aspect='auto', origin='lower', cmap='YlGn',
               extent=[0, n_times, 0, n_freqs])
axes[1].set_title('White Noise (null)', color='#52b788', fontsize=9)
axes[2].bar(['Biotic', 'Null'], [ACI_val, ACI_null],
            color=['#52b788', '#4a7a5a'], width=0.5)
axes[2].set_title('ACI', color='#52b788', fontsize=9)
plt.tight_layout()
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=130, bbox_inches='tight', facecolor='#0a1a0f')
plt.close(); buf.seek(0)
display_plot(base64.b64encode(buf.read()).decode())

Try: change n_species or noise_level

Experiment 2 — Dawn Acoustic Duct & Active SpaceAtmospheric Acoustics

Model sound propagation at dawn (temperature inversion → acoustic duct) vs midday. Compute and plot active space radius as a function of ambient noise.

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import io, base64

def c(T_C): return 331.3 * np.sqrt(1 + T_C / 273.15)

z = np.linspace(0, 50, 300)
c_dawn   = c(10 + 0.15 * z)
c_midday = c(25 - 0.006 * z)

def active_space(I0_dB, N_dB, beta):
    I0, N = 10**(I0_dB/10), 10**(N_dB/10)
    r = (I0 / N)**(1/beta); return r, np.pi*r**2

source_dB = 70
r_d, A_d = active_space(source_dB, 30, beta=1)
r_m, A_m = active_space(source_dB, 30, beta=2)
print(f"Dawn  (duct):    r = {r_d:.1f} m, A = {A_d/1e4:.2f} ha")
print(f"Midday (spher): r = {r_m:.1f} m, A = {A_m/1e4:.2f} ha")
print(f"Area ratio dawn/midday: {A_d/A_m:.1f}x")

noise = np.linspace(20, 60, 100)
r_d_arr = np.array([active_space(source_dB,n,1)[0] for n in noise])
r_m_arr = np.array([active_space(source_dB,n,2)[0] for n in noise])

fig, axes = plt.subplots(1, 3, figsize=(13, 4.5), facecolor='#0a1a0f')
fig.suptitle('Dawn Acoustic Duct - Active Space Amplification',
             color='#c8e0cc', fontsize=11, fontfamily='monospace')
for ax in axes:
    ax.set_facecolor('#040d07')
    for sp in ax.spines.values(): sp.set_color('#2d6a4f')
    ax.tick_params(colors='#6a9a7a', labelsize=8)
axes[0].plot(c_dawn, z, color='#f4a261', lw=2, label='Dawn')
axes[0].plot(c_midday, z, color='#a8dadc', lw=2, label='Midday')
axes[0].set_title('Atmospheric Profile', color='#52b788', fontsize=9)
axes[0].legend(fontsize=7, facecolor='#0a1a0f', edgecolor='#2d6a4f', labelcolor='#c8e0cc')
axes[1].plot(noise, r_d_arr, color='#f4a261', lw=2, label='Dawn (duct, β=1)')
axes[1].plot(noise, r_m_arr, color='#a8dadc', lw=2, label='Midday (β=2)')
axes[1].set_title('Active Space Radius', color='#52b788', fontsize=9)
axes[1].legend(fontsize=7, facecolor='#0a1a0f', edgecolor='#2d6a4f', labelcolor='#c8e0cc')
axes[2].plot(noise, (r_d_arr/r_m_arr)**2, color='#e9c46a', lw=2)
axes[2].set_title('Amplification Factor', color='#52b788', fontsize=9)
plt.tight_layout()
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=130, bbox_inches='tight', facecolor='#0a1a0f')
plt.close(); buf.seek(0)
display_plot(base64.b64encode(buf.read()).decode())

Try: change source_dB or the inversion gradient

Experiment 3 — Dawn Chorus Species Sequence & Niche OverlapCommunity Ecology

Model staggered species onset times, their spectral niches, and compute pairwise acoustic niche overlap O_ij. Visualise community structure.

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import io, base64

species = [
    ("Robin", 90, 3.5, 1.2, '#9ecae1'),
    ("Blackbird", 75, 2.8, 1.8, '#f4a261'),
    ("Song Thrush", 60, 4.2, 1.5, '#52b788'),
    ("Wren", 45, 5.5, 2.0, '#a8dadc'),
    ("Great Tit", 30, 4.8, 1.0, '#e9c46a'),
    ("Chaffinch", 15, 3.2, 0.8, '#e76f51'),
]
T = np.linspace(0, 120, 500); F = np.linspace(1, 10, 300); sunrise = 90
S_total = np.zeros((len(F), len(T))); S_ind = []
for name, onset, fc, bw, col in species:
    t0 = sunrise - onset
    t_env = np.zeros(len(T))
    active = (T >= t0) & (T <= sunrise + 20)
    if active.any():
        t_env[active] = np.exp(-0.5 * ((T[active] - (t0 + 15)) / 12)**2)
    f_env = np.exp(-0.5 * ((F - fc) / (bw/2))**2)
    layer = np.outer(f_env, t_env)
    S_ind.append(layer); S_total += layer

n = len(species); O = np.zeros((n,n))
for i in range(n):
    for j in range(n):
        num = np.sum(S_ind[i] * S_ind[j])
        den = np.sqrt(np.sum(S_ind[i]**2) * np.sum(S_ind[j]**2)) + 1e-12
        O[i,j] = num/den
print("Niche-overlap matrix O_ij:")
names = [s[0][:8] for s in species]
print("         " + "  ".join(f"{n:8s}" for n in names))
for i,row in enumerate(O):
    print(f"{names[i]:8s} " + "  ".join(f"{v:8.4f}" for v in row))
print(f"Mean off-diagonal overlap: {(O.sum()-np.trace(O))/(n*(n-1)):.4f}")

fig = plt.figure(figsize=(14,5), facecolor='#0a1a0f')
gs = fig.add_gridspec(1,3, width_ratios=[4,3,2.5], wspace=0.35)
fig.suptitle('Dawn Chorus - Spectrotemporal Niche Structure',
             color='#c8e0cc', fontsize=11, fontfamily='monospace')
ax1 = fig.add_subplot(gs[0]); ax1.set_facecolor('#040d07')
ax1.imshow(S_total, aspect='auto', origin='lower', cmap='YlGn',
           extent=[0,120,1,10], alpha=0.85)
ax1.axvline(sunrise, color='#f4a261', lw=1.5, ls='--', alpha=0.7)
ax1.set_xlabel('Min before/after sunrise', color='#6a9a7a', fontsize=8)
ax1.set_ylabel('Frequency (kHz)', color='#6a9a7a', fontsize=8)
ax1.set_title('Spectrotemporal Activity', color='#52b788', fontsize=9)
for sp in ax1.spines.values(): sp.set_color('#2d6a4f')
ax1.tick_params(colors='#6a9a7a', labelsize=8)
ax2 = fig.add_subplot(gs[1]); ax2.set_facecolor('#040d07')
y_pos = np.arange(len(species))
ax2.barh(y_pos, [s[1] for s in species], color=[s[4] for s in species], alpha=0.8, height=0.6)
ax2.set_yticks(y_pos); ax2.set_yticklabels([s[0] for s in species], fontsize=8, color='#c8e0cc')
ax2.set_title('Onset Sequence', color='#52b788', fontsize=9); ax2.invert_xaxis()
for sp in ax2.spines.values(): sp.set_color('#2d6a4f')
ax2.tick_params(colors='#6a9a7a', labelsize=8)
ax3 = fig.add_subplot(gs[2]); ax3.set_facecolor('#040d07')
im = ax3.imshow(O, cmap='Greens', vmin=0, vmax=1, aspect='auto')
ax3.set_xticks(range(n)); ax3.set_yticks(range(n))
ax3.set_xticklabels([s[0][:4] for s in species], rotation=45, ha='right', fontsize=7, color='#c8e0cc')
ax3.set_yticklabels([s[0][:4] for s in species], fontsize=7, color='#c8e0cc')
ax3.set_title('Niche Overlap O_ij', color='#52b788', fontsize=9)
for sp in ax3.spines.values(): sp.set_color('#2d6a4f')
plt.tight_layout()
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=130, bbox_inches='tight', facecolor='#0a1a0f')
plt.close(); buf.seek(0)
display_plot(base64.b64encode(buf.read()).decode())

Try: shift onset times or f_center values

Experiment 4 — Optimal Transport & Acoustic Restoration (W₂)Mathematical Physics

Compute the Wasserstein-2 distance between degraded and reference acoustic spectral distributions. Model the restoration trajectory and entropy production rate.

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.stats import wasserstein_distance
import io, base64

np.random.seed(7); F = np.linspace(0.2, 12, 200)

def healthy(F):
    s = 0.05*np.ones_like(F)
    for fc,bw,amp in [(2.5,0.6,1.0),(3.8,0.8,1.3),(5.2,0.5,0.8),(4.5,1.0,0.6),(7.0,0.4,0.5),(1.5,0.4,0.4)]:
        s += amp*np.exp(-0.5*((F-fc)/bw)**2)
    return s/s.sum()

def degraded(F):
    s = 0.3*np.exp(-0.5*((F-0.8)/0.4)**2)
    s += 0.15*np.exp(-0.5*((F-0.4)/0.2)**2) + 0.02*np.ones_like(F)
    return s/s.sum()

mu_ref, mu_deg = healthy(F), degraded(F)
print(f"W1(mu_deg, mu_ref) = {wasserstein_distance(F,F,mu_deg,mu_ref):.4f} kHz")

n_steps = 7; alphas = np.linspace(0,1,n_steps)
def interp(mu0,mu1,alpha,F):
    c0,c1 = np.cumsum(mu0)/mu0.sum(), np.cumsum(mu1)/mu1.sum()
    qg = np.linspace(0,1,len(F))
    Q = (1-alpha)*np.interp(qg,c0,F) + alpha*np.interp(qg,c1,F)
    mu = np.zeros_like(F)
    for q in Q: mu[np.argmin(np.abs(F-q))] += 1/len(F)
    return mu/(mu.sum()+1e-12)

spectra, W_arr, H_arr = [], [], []
for a in alphas:
    mu = interp(mu_deg, mu_ref, a, F); spectra.append(mu)
    W_arr.append(wasserstein_distance(F,F,mu,mu_ref))
    pp = mu[mu>0]; H_arr.append(-np.sum(pp*np.log(pp)))

print(f"\n{'alpha':>6} {'W1':>8} {'H':>8}")
for a,w,h in zip(alphas,W_arr,H_arr): print(f"{a:6.2f} {w:8.4f} {h:8.4f}")

fig, axes = plt.subplots(1,3, figsize=(14,4.5), facecolor='#0a1a0f')
fig.suptitle('Optimal Transport - Acoustic Restoration Trajectory (W2)',
             color='#c8e0cc', fontsize=11, fontfamily='monospace')
cmap_s = plt.cm.YlGn(np.linspace(0.2,1.0,n_steps))
for ax in axes:
    ax.set_facecolor('#040d07')
    for sp in ax.spines.values(): sp.set_color('#2d6a4f')
    ax.tick_params(colors='#6a9a7a', labelsize=8)
for i,(mu,a) in enumerate(zip(spectra,alphas)):
    lw = 2.5 if i in (0,n_steps-1) else 0.9
    axes[0].plot(F, mu, color=cmap_s[i], lw=lw)
axes[0].plot(F, mu_ref, color='#f4a261', lw=2, ls='--', label='Reference')
axes[0].plot(F, mu_deg, color='#e76f51', lw=2, ls='--', label='Degraded')
axes[0].set_title('Wasserstein Geodesic', color='#52b788', fontsize=9)
axes[0].legend(fontsize=7, facecolor='#0a1a0f', edgecolor='#2d6a4f', labelcolor='#c8e0cc')
axes[1].plot(alphas, W_arr, 'o-', color='#52b788', lw=2, ms=6)
axes[1].set_title('Distance to Reference', color='#52b788', fontsize=9)
axes[2].plot(alphas, H_arr, 's-', color='#e9c46a', lw=2, ms=6)
axes[2].set_title('Acoustic Entropy Recovery', color='#52b788', fontsize=9)
plt.tight_layout()
buf = io.BytesIO()
plt.savefig(buf, format='png', dpi=130, bbox_inches='tight', facecolor='#0a1a0f')
plt.close(); buf.seek(0)
display_plot(base64.b64encode(buf.read()).decode())

Direct connection to Villani T₂ / HWI framework

The four labs at a glance

Lab 1 — Spectrogram of a synthesised dawn chorus

Build a multi-species chorus by adding sinusoids of different frequencies and time envelopes, compute its short-time Fourier transform, and visualise the resulting spectrogram. Compare it to a real bird recording (provided as a small embedded WAV). The lab introduces SciPy's scipy.signal.spectrogramand the trade-off between time and frequency resolution governed by the window length.

Lab 2 — Acoustic indices from scratch

Implement the ACI (Acoustic Complexity Index), BI (Bioacoustic Index) and NDSI (Normalised Difference Soundscape Index) directly from their definitions in NumPy. Apply them to a set of test spectrograms with known soundscape composition and verify that the indices respond as expected. This lab is a good entry point to the technical literature on soundscape ecology.

Lab 3 — Niche-partitioning model

Simulate N species as bivariate Gaussian distributions in frequency-time space; compute pairwise acoustic-niche overlaps O_ij; and investigate how community-level overlap evolves under a Lotka-Volterra-like competition dynamic where high-overlap pairs penalise each other. Plot the resulting partition and verify Krause's prediction of evolved non-overlap.

Lab 4 — Optimal-transport restoration

Take a degraded soundscape and a healthy reference; compute the W₂ Wasserstein distance between their spectral measures; and find the optimal transport plan that maps one onto the other. This is the mathematical content of Module 11 made concrete with the POT (Python Optimal Transport) library. The lab illustrates how acoustic restoration can be quantified, not just described.

What you need

A modern browser (Chrome, Firefox, Safari, Edge — Pyodide works on all four).
~100 MB of free memory for the Pyodide runtime.
~30 s of patience for Pyodide's first download — the runtime is cached for subsequent visits.

Tip

Labs are independent: you can jump straight to Lab 4 if you have the prerequisites, or work through them in order to build up the workflow. All four labs persist their cell state in browser local storage, so you can leave and come back without losing work.

Share:X Reddit LinkedIn