Flying & singing zebra finches - Birdpark#

Data from Rüttimann et al. [2025] (paper, zenodo), a multimodal dataset of zebra finch groups with synchronized video, microphone arrays, and backpack-mounted vibration transducer (accelerometers).

The code below shows how one can convert that existing dataset into the Trials.nc format. The sampling rate of the vibration transducer is very high (24kHz), so the raw vibration trace loads slowly; plot the downsampled vibration_env envelope instead (built below, together with candidate changepoints on its bursts), or use the Downsample option next to Load.

Left: GUI screenshot, Right: Adapted from Rüttimann et al. [2025], Fig. 2C

import zipfile
from pathlib import Path

import h5py
import numpy as np
import pandas as pd
import requests
import xarray as xr
from audioio import write_audio

import ethograph as eto
from ethograph.features.changepoints import find_nearest_turning_points_binary
from ethograph.features.energy import bandpass_envelope, highpass_envelope
from ethograph.io.pairing import pair_media

Download dataset#

You can download the entire dataset from here or use the code below. I only tested the copExpBP08 recording. If problems arise, the ReadMe here is very helpful.

try:
    _here = Path(__vsc_ipynb_file__).parent
except NameError:
    _here = Path().resolve()

data_folder = _here.parent / "data" / "birdpark"
data_folder.mkdir(parents=True, exist_ok=True)

response = requests.get("https://zenodo.org/api/records/13144875")
data = response.json()

# Download dataset if not already present
if not (data_folder / "copExpBP08" / "BP_2021-05-25_08-12-51_655154_0380000.mp4").exists():
    for file in data["files"]:
        if file["checksum"] == "md5:32d1ae6049556c803f68b6d354c952ca":
            print(f"Checksum matches: {file['key']}")
            output_path = data_folder / file["key"]
            r = requests.get(file["links"]["self"], stream=True)
            with open(output_path, "wb") as f:
                for chunk in r.iter_content(chunk_size=8192):
                    f.write(chunk)
            if output_path.suffix == ".zip":
                with zipfile.ZipFile(output_path, "r") as zip_ref:
                    zip_ref.extractall(data_folder)
            break

Build NWB alignment and dataset#

# Recording: copExpBP08, session BP_2021-05-25_08-12-51_655154_0380000
fps = 47.6837158203125
audio_sr = 24414.0625  # audio and accelerometer sampling rate (from dataset metadata)

recording_folder = data_folder / "copExpBP08"
h5_path = recording_folder / "BP_2021-05-25_08-12-51_655154_0380000.h5"
video_path = recording_folder / "BP_2021-05-25_08-12-51_655154_0380000.mp4"
audio_path = video_path.with_suffix(".wav")
nc_path = video_path.with_suffix(".nc")

# Read H5 file
with h5py.File(h5_path, "r") as f1:
    radioSignals = f1["/radioSignals"][()]  # accelerometer (one row per channel)
    daqSignals = f1["/daqSignals"][()]  # microphone channels

# Create .wav file from microphone channels
write_audio(audio_path, daqSignals.T, audio_sr)

# ─── Build session table ───
# trial=1 matches the ID assigned when loading a plain Dataset as a single-trial TrialTree
# The table names each file by its full path; pair_media stores the basename per
# trial and the path per stream, so the GUI finds the media without being told
# the folder (and a folder setting still overrides it on another machine).
session_table = pd.DataFrame(
    {
        "trial": [1],
        "video_0": [str(video_path)],
        "audio_0": [str(audio_path)],
    }
)

nwb_path = recording_folder / ".ethograph" / "alignment.nwb"
pair_media(
    trial_table=session_table,
    stream_rates={"video": float(fps), "audio": float(audio_sr)},
    output_path=nwb_path,
    on_existing="replace"
)

# ─── Build xarray dataset ───
time_coords = np.arange(radioSignals.shape[1]) / audio_sr

ds = xr.Dataset(
    data_vars={
        "vibration": xr.DataArray(
            radioSignals.T,
            dims=["time", "individual"],
        ),
    },
    coords={
        "time": time_coords,
        "individual": ["male (red radio)", "female (yellow radio)"],  # specific to copExpBP08
    },
    attrs={
        "fps": fps,
        "audio_sr": audio_sr,
    },
)

# ─── Amplitude envelopes + changepoints ───
# Envelopes have their own clocks; any coord containing "time" is recognised by the GUI.
ds["vibration"] = ds.vibration.where(ds.vibration > -999_999)

envelopes, envelopes_high = [], []
for ind in ds.individual.values:
    raw = ds.vibration.sel(individual=ind)
    filled = raw.fillna(raw.median()).values
    env_time, env = bandpass_envelope(filled, audio_sr, band=(10.0, 5000.0), env_rate=200.0, cutoff=20.0)
    env_high_time, env_high = highpass_envelope(filled, audio_sr, env_rate=1000.0, cutoff=200.0)
    envelopes.append(env)
    envelopes_high.append(env_high)

ds["vibration_bandpass"] = xr.DataArray(
    np.stack(envelopes, axis=1), dims=["time_env", "individual"], coords={"time_env": env_time}
)
ds["vibration_highpass"] = xr.DataArray(
    np.stack(envelopes_high, axis=1), dims=["time_env_high", "individual"], coords={"time_env_high": env_high_time}
)

env_scale = float(np.nanpercentile(ds.vibration_highpass.values, 99))
ds = eto.add_changepoints_to_ds(
    ds=ds,
    target_feature="vibration_highpass",
    changepoint_name="turning_points",
    changepoint_func=find_nearest_turning_points_binary,
    threshold=0.02 * env_scale,
    max_value=0.1 * env_scale,
    prominence=0.2 * env_scale,
    distance=int(0.05 * 1000.0),
)

ds.to_netcdf(nc_path)
print(f"Saved to {nc_path}")

Saved to c:\Users\aksel\Documents\Code\ethograph\data\birdpark\copExpBP08\BP_2021-05-25_08-12-51_655154_0380000.nc

# ─── Microphone envelopes + sliding-window correlation with the vibration ───
# Each microphone gets the same highpass envelope as the vibration, so both live on
# the 1000 Hz `time_env_high` clock and can be correlated sample-for-sample.
corr_win_s = 0.25
env_rate_high = 1000.0

mic_names = [f"mic_{i}" for i in range(daqSignals.shape[0])]
audio_envs = []
for channel in daqSignals:
    _, mic_env = highpass_envelope(channel.astype(float), audio_sr, env_rate=env_rate_high, cutoff=200.0)
    audio_envs.append(mic_env)

vib_high = ds.vibration_highpass.values
n = min(vib_high.shape[0], min(len(e) for e in audio_envs))
audio_env = np.stack([e[:n] for e in audio_envs], axis=1)
vib_high = vib_high[:n]
env_high_time_c = ds.time_env_high.values[:n]


def sliding_corr(a: np.ndarray, b: np.ndarray, window: int) -> np.ndarray:
    """Pearson correlation of ``a`` and ``b`` over a centred sliding window."""
    kernel = np.ones(window)

    def roll(x):
        return np.convolve(x, kernel, mode="same")

    n_win = roll(np.ones_like(a))
    mean_a, mean_b = roll(a) / n_win, roll(b) / n_win
    cov = roll(a * b) / n_win - mean_a * mean_b
    var_a = np.maximum(roll(a * a) / n_win - mean_a**2, 0.0)
    var_b = np.maximum(roll(b * b) / n_win - mean_b**2, 0.0)
    denom = np.sqrt(var_a * var_b)
    return np.divide(cov, denom, out=np.zeros_like(cov), where=denom > 0)


window = int(round(corr_win_s * env_rate_high))
pair_names, pair_corrs = [], []
for i, ind in enumerate(ds.individual.values):
    for j, mic in enumerate(mic_names):
        pair_names.append(f"{ind} × {mic}")
        pair_corrs.append(sliding_corr(vib_high[:, i], audio_env[:, j], window))

ds["audio_env"] = xr.DataArray(
    audio_env, dims=["time_env_high", "microphone"], coords={"time_env_high": env_high_time_c, "microphone": mic_names}
)
ds["sliding_corr"] = xr.DataArray(
    np.stack(pair_corrs, axis=1),
    dims=["time_env_high", "pair"],
    coords={"time_env_high": env_high_time_c, "pair": pair_names},
)

ds.to_netcdf(nc_path)
print(f"Saved {list(ds.data_vars)} to {nc_path}")

Decent for segmentation#

voc.segment.meansquared(
    <data>,  # set by GUI
    <sr>,    # set by GUI
    threshold=15000,
    min_dur=0.003,
    min_silent_dur=0.0001,
    freq_cutoffs=(500, 10000),
    smooth_win=0.32,
    scale=True,
    scale_val=32768,
)