API reference

Sessions and I/O

Session discovery and BWF metadata for AmbiX recordings.

A session is a folder of WAV files from one recording occasion. Files whose BWF timestamps chain end-to-start (recorder 2 GB splits) are treated as one continuous take; otherwise they are separate takes on a common absolute timeline (seconds since the session's first midnight).

Take dataclass

Source code in src/ambiscape/io.py

@dataclass
class Take:
    path: Path
    start: float          # seconds since session day 0 midnight
    duration: float
    frames: int
    samplerate: int
    channels: int
    date: str
    clock: str
    order: str = "ambix"  # 'ambix' (W,Y,Z,X) or 'fuma' (W,X,Y,Z)

    @property
    def end(self) -> float:
        return self.start + self.duration

    @property
    def wyzx(self) -> tuple[int, int, int, int]:
        """Column indices of (W, Y, Z, X) for this take's convention."""
        return (0, 2, 3, 1) if self.order == "fuma" else (0, 1, 2, 3)

wyzx property

Column indices of (W, Y, Z, X) for this take's convention.

Session dataclass

Source code in src/ambiscape/io.py

@dataclass
class Session:
    folder: Path
    takes: list[Take] = field(default_factory=list)
    day0: _dt.date | None = None

    @property
    def duration(self) -> float:
        return sum(t.duration for t in self.takes)

    @property
    def name(self) -> str:
        return getattr(self, "_name", None) or self.folder.name

    def clock(self, t: float) -> str:
        """Absolute seconds -> 'DD Mon HH:MM:SS' string."""
        base = _dt.datetime.combine(self.day0, _dt.time())
        return (base + _dt.timedelta(seconds=t)).strftime("%d %b %H:%M:%S")

clock(t)

Absolute seconds -> 'DD Mon HH:MM:SS' string.

Source code in src/ambiscape/io.py

def clock(self, t: float) -> str:
    """Absolute seconds -> 'DD Mon HH:MM:SS' string."""
    base = _dt.datetime.combine(self.day0, _dt.time())
    return (base + _dt.timedelta(seconds=t)).strftime("%d %b %H:%M:%S")

read_bext(path)

Parse the BWF 'bext' chunk (pure python RIFF walk).

Source code in src/ambiscape/io.py

def read_bext(path: str | Path) -> dict:
    """Parse the BWF 'bext' chunk (pure python RIFF walk)."""
    out = {}
    with open(path, "rb") as f:
        riff, _size, wave = struct.unpack("<4sI4s", f.read(12))
        if riff != b"RIFF" or wave != b"WAVE":
            raise ValueError(f"{path}: not a RIFF/WAVE file")
        while True:
            hdr = f.read(8)
            if len(hdr) < 8:
                break
            cid, csize = struct.unpack("<4sI", hdr)
            if cid == b"bext":
                data = f.read(min(csize, 604))
                out["description"] = data[0:256].split(b"\0")[0].decode("ascii", "replace")
                out["originator"] = data[256:288].split(b"\0")[0].decode("ascii", "replace")
                out["date"] = data[320:330].decode("ascii", "replace").strip("\0 ")
                out["time"] = data[330:338].decode("ascii", "replace").strip("\0 ")
                out["time_reference"] = struct.unpack("<Q", data[338:346])[0]
                break
            f.seek(csize + (csize & 1), 1)
    return out

channel_order(bext_description)

Detect B-format convention from the H3-VR's zTRK tags in the bext description: 'ambix' (W,Y,Z,X) or 'fuma' (W,X,Y,Z). Defaults to 'ambix' when no tags are present.

Source code in src/ambiscape/io.py

def channel_order(bext_description: str) -> str:
    """Detect B-format convention from the H3-VR's zTRK tags in the bext
    description: 'ambix' (W,Y,Z,X) or 'fuma' (W,X,Y,Z). Defaults to 'ambix'
    when no tags are present."""
    trk = {}
    for line in bext_description.replace("\r", "\n").split("\n"):
        if line.startswith("zTRK") and "=" in line:
            k, v = line.split("=", 1)
            trk[int(k[4:])] = v.strip().upper()
    seq = [trk.get(i) for i in (1, 2, 3, 4)]
    if seq == ["W", "X", "Y", "Z"]:
        return "fuma"
    return "ambix"

open_session(folder)

Scan a session folder.

If calibration.json contains clock_offset_s, that many seconds are added to every take's start time — the fix for a recorder whose clock was found to be off (positive offset = clock was slow). All clock-labeled outputs (figures, annotations, reports) then agree on corrected time.

Source code in src/ambiscape/io.py

def open_session(folder: str | Path) -> Session:
    """Scan a session folder.

    If ``calibration.json`` contains ``clock_offset_s``, that many seconds are
    added to every take's start time — the fix for a recorder whose clock was
    found to be off (positive offset = clock was slow). All clock-labeled
    outputs (figures, annotations, reports) then agree on corrected time.
    """
    folder = Path(folder)
    paths = sorted(p for p in folder.iterdir()
                   if p.suffix.lower() == ".wav" and p.is_file())
    if not paths:
        raise FileNotFoundError(f"no WAV files in {folder}")
    clock_offset = 0.0
    cal = folder / "calibration.json"
    if cal.exists():
        import json
        clock_offset = float(json.loads(cal.read_text())
                             .get("clock_offset_s", 0.0))
    sess = Session(folder=folder)
    for p in paths:
        info = sf.info(str(p))
        bx = read_bext(p)
        date = _dt.date.fromisoformat(bx["date"].replace(":", "-"))
        if sess.day0 is None:
            sess.day0 = date
        hh, mm, ss = (int(x) for x in bx["time"].split(":"))
        start = ((date - sess.day0).days * 86400 + hh * 3600 + mm * 60 + ss
                 + clock_offset)
        sess.takes.append(Take(
            path=p, start=float(start), duration=info.frames / info.samplerate,
            frames=info.frames, samplerate=info.samplerate,
            channels=info.channels, date=bx["date"], clock=bx["time"],
            order=channel_order(bx.get("description", "")),
        ))
    sess.takes.sort(key=lambda t: t.start)
    return sess

open_recording(path)

Open a single WAV file as a one-take session ("scene").

The folder-as-session model of :func:open_session assumes every WAV in a folder belongs to one recording occasion on a shared clock. A contributed corpus is often the opposite: one folder per recordist, each holding many independent one-off scenes from different places and dates. This opens exactly one file as its own session (day0 = that file's BWF date), so each scene can go through the full pipeline on its own. The session name is the file stem.

Source code in src/ambiscape/io.py

def open_recording(path: str | Path) -> Session:
    """Open a single WAV file as a one-take session ("scene").

    The folder-as-session model of :func:`open_session` assumes every WAV in
    a folder belongs to one recording occasion on a shared clock. A
    contributed corpus is often the opposite: one folder per recordist, each
    holding many independent one-off scenes from different places and dates.
    This opens exactly one file as its own session (day0 = that file's BWF
    date), so each scene can go through the full pipeline on its own. The
    session name is the file stem.
    """
    path = Path(path)
    if not path.is_file():
        raise FileNotFoundError(path)
    info = sf.info(str(path))
    bx = read_bext(path)
    date = _dt.date.fromisoformat(bx["date"].replace(":", "-"))
    hh, mm, ss = (int(x) for x in bx["time"].split(":"))
    sess = Session(folder=path.parent, day0=date)
    sess.takes.append(Take(
        path=path, start=float(hh * 3600 + mm * 60 + ss),
        duration=info.frames / info.samplerate, frames=info.frames,
        samplerate=info.samplerate, channels=info.channels,
        date=bx["date"], clock=bx["time"],
        order=channel_order(bx.get("description", "")),
    ))
    sess._name = path.stem
    return sess

read_span(sess, t0, dur, dtype='float32')

Read [t0, t0+dur) seconds (session time) from whichever take covers it.

Source code in src/ambiscape/io.py

def read_span(sess: Session, t0: float, dur: float, dtype="float32"):
    """Read [t0, t0+dur) seconds (session time) from whichever take covers it."""
    for tk in sess.takes:
        if tk.start <= t0 < tk.end:
            fs = tk.samplerate
            off = int((t0 - tk.start) * fs)
            n = min(int(dur * fs), tk.frames - off)
            with sf.SoundFile(str(tk.path)) as f:
                f.seek(off)
                return f.read(n, dtype=dtype, always_2d=True), fs
    raise ValueError(f"t={t0} not covered by session {sess.name}")

export_segment(sess, t0, dur, out_path)

Bit-exact 4-channel excerpt [t0, t0+dur) to a WAV.

Samples are copied in the source's own PCM subtype (no float round trip), so the excerpt is archival: the representative segments of a report stay citable against the raw takes. The span must lie within one take (recorder 2 GB splits chain seamlessly only in read_span's float path). Returns the output path.

Source code in src/ambiscape/io.py

def export_segment(sess: Session, t0: float, dur: float,
                   out_path: str | Path) -> Path:
    """Bit-exact 4-channel excerpt [t0, t0+dur) to a WAV.

    Samples are copied in the source's own PCM subtype (no float round
    trip), so the excerpt is archival: the representative segments of a
    report stay citable against the raw takes. The span must lie within one
    take (recorder 2 GB splits chain seamlessly only in ``read_span``'s
    float path). Returns the output path.
    """
    out_path = Path(out_path)
    for tk in sess.takes:
        if tk.start <= t0 < tk.end:
            fs = tk.samplerate
            off = int((t0 - tk.start) * fs)
            n = min(int(dur * fs), tk.frames - off)
            with sf.SoundFile(str(tk.path)) as f:
                subtype = f.subtype
                dtype = "int16" if subtype == "PCM_16" else "int32"
                f.seek(off)
                data = f.read(n, dtype=dtype, always_2d=True)
            sf.write(str(out_path), data, fs, subtype=subtype)
            return out_path
    raise ValueError(f"t={t0} not covered by session {sess.name}")

stereo_preview(x, wyzx=(0, 1, 2, 3), az_deg=90.0)

Side-facing cardioid stereo decode of an AmbiX block, for previews.

Left/right cardioids at ±az_deg in the horizontal plane: 0.5 * (W ± sin(az) * Y) (SN3D). Returns an (n, 2) float array — write it with soundfile for a listenable preview of an exported segment.

Source code in src/ambiscape/io.py

def stereo_preview(x, wyzx=(0, 1, 2, 3), az_deg: float = 90.0):
    """Side-facing cardioid stereo decode of an AmbiX block, for previews.

    Left/right cardioids at ±``az_deg`` in the horizontal plane:
    ``0.5 * (W ± sin(az) * Y)`` (SN3D). Returns an (n, 2) float array —
    write it with ``soundfile`` for a listenable preview of an exported
    segment.
    """
    import numpy as np
    W, Y = x[:, wyzx[0]], x[:, wyzx[1]]
    g = float(np.sin(np.radians(az_deg)))
    return np.stack([0.5 * (W + g * Y), 0.5 * (W - g * Y)], axis=1)

Feature extraction

Streaming per-second feature extraction from AmbiX B-format WAV files.

Designed for arbitrarily long recordings: files are read in 60-s blocks and never held in memory. Per second: broadband and A-weighted fast levels (125 ms), octave-band powers, spectral centroid/flatness, a 96-band log-frequency spectrogram row, per-octave pseudo-intensity vectors, broadband DOA (azimuth, elevation) and diffuseness. Per minute: full-resolution mean PSD (for narrowband hum tracking and fingerprinting).

Channel order is AmbiX ACN (W, Y, Z, X) as written by the Zoom H3-VR.

a_weighting_sos(fs)

IEC 61672 A-weighting as SOS (bilinear transform of the analog filter).

Source code in src/ambiscape/features.py

def a_weighting_sos(fs: int):
    """IEC 61672 A-weighting as SOS (bilinear transform of the analog filter)."""
    f1, f2, f3, f4 = 20.598997, 107.65265, 737.86223, 12194.217
    a1000 = 1.9997
    nums = [(2 * np.pi * f4) ** 2 * 10 ** (a1000 / 20), 0, 0, 0, 0]
    dens = np.polymul(
        np.polymul([1, 4 * np.pi * f4, (2 * np.pi * f4) ** 2],
                   [1, 4 * np.pi * f1, (2 * np.pi * f1) ** 2]),
        np.polymul([1, 2 * np.pi * f3], [1, 2 * np.pi * f2]),
    )
    b, a = signal.bilinear(nums, dens, fs)
    return signal.tf2sos(b, a)

extract_take(take, verbose=False)

Run the streaming extractor over one file; returns feature arrays.

Source code in src/ambiscape/features.py

def extract_take(take: Take, verbose: bool = False) -> dict:
    """Run the streaming extractor over one file; returns feature arrays."""
    fs = take.samplerate
    hop = int(HOP * fs / 48000)
    nfft = NFFT
    win = np.hanning(nfft).astype(np.float32)
    wsum2 = float((win ** 2).sum())
    freqs = np.fft.rfftfreq(nfft, 1 / fs)
    oct_idx = [np.where((freqs >= c / np.sqrt(2)) & (freqs < c * np.sqrt(2)))[0]
               for c in OCT_CENTERS]
    logf = np.geomspace(*LOGF_RANGE, N_LOGBANDS + 1)
    log_idx = np.clip(np.searchsorted(logf, freqs) - 1, -1, N_LOGBANDS - 1)
    doa_mask = (freqs >= DOA_BAND[0]) & (freqs <= DOA_BAND[1])
    spec_mask = (freqs >= 50) & (freqs <= 16000)
    a_sos = a_weighting_sos(fs)
    a_state = np.zeros((a_sos.shape[0], 2), dtype=np.float64)

    nsec = int(take.frames // fs)
    nfast = int(take.frames // int(FAST * fs))
    ffs = int(FAST * fs)
    hfs = int(HI_ENV * fs)
    nhi = int(take.frames // hfs)
    F = {
        "fast_db": np.zeros(nfast, np.float32),
        "fast_dba": np.zeros(nfast, np.float32),
        "env_hi": np.zeros(nhi, np.float32),
        "rms_w": np.zeros(nsec, np.float32),
        "peak": np.zeros(nsec, np.float32),
        "oct_pow": np.zeros((nsec, len(OCT_CENTERS)), np.float32),
        "centroid": np.zeros(nsec, np.float32),
        "flatness": np.zeros(nsec, np.float32),
        "logspec": np.zeros((nsec, N_LOGBANDS), np.float32),
        "I_band": np.zeros((nsec, len(OCT_CENTERS), 3), np.float32),
        "az": np.zeros(nsec, np.float32),
        "el": np.zeros(nsec, np.float32),
        "diffuse": np.zeros(nsec, np.float32),
    }
    nmin = -(-nsec // 60) if nsec else 0
    minspec = np.zeros((nmin, len(freqs)), np.float64)
    mincnt = np.zeros(nmin, np.int64)
    eps = 1e-20

    carry = np.zeros((0, 4), np.float32)
    sec_base = 0
    fast_base = 0
    hi_base = 0
    with sf.SoundFile(str(take.path)) as f:
        while True:
            block = f.read(60 * fs, dtype="float32", always_2d=True)
            if block.shape[0] == 0:
                break
            data = np.concatenate([carry, block]) if carry.shape[0] else block
            navail = data.shape[0]
            nsec_blk = min(navail // fs, nsec - sec_base)
            nwin = (navail - nfft) // hop + 1 if navail >= nfft else 0
            if nsec_blk <= 0:
                break

            # fast levels on W (contiguous 125 ms frames)
            nfast_blk = min((nsec_blk * fs) // ffs, nfast - fast_base)
            wseg = data[: nfast_blk * ffs, 0].reshape(nfast_blk, ffs)
            F["fast_db"][fast_base:fast_base + nfast_blk] = 10 * np.log10(
                (wseg.astype(np.float64) ** 2).mean(1) + eps)
            wa, a_state = signal.sosfilt(a_sos, data[: nfast_blk * ffs, 0]
                                         .astype(np.float64), zi=a_state)
            F["fast_dba"][fast_base:fast_base + nfast_blk] = 10 * np.log10(
                (wa.reshape(nfast_blk, ffs) ** 2).mean(1) + eps)
            fast_base += nfast_blk

            # 20 ms broadband envelope on W (linear power, for modulation)
            nhi_blk = min((nsec_blk * fs) // hfs, nhi - hi_base)
            hseg = data[: nhi_blk * hfs, 0].reshape(nhi_blk, hfs)
            F["env_hi"][hi_base:hi_base + nhi_blk] = \
                (hseg.astype(np.float64) ** 2).mean(1)
            hi_base += nhi_blk

            if nwin > 0:
                iW, iY, iZ, iX = take.wyzx
                idx = np.arange(nfft)[None, :] + hop * np.arange(nwin)[:, None]
                Wf = np.fft.rfft(data[:, iW][idx] * win)
                Yf = np.fft.rfft(data[:, iY][idx] * win)
                Zf = np.fft.rfft(data[:, iZ][idx] * win)
                Xf = np.fft.rfft(data[:, iX][idx] * win)
                Pw = (Wf.real ** 2 + Wf.imag ** 2) / wsum2
                IX = (Wf.conj() * Xf).real / wsum2
                IY = (Wf.conj() * Yf).real / wsum2
                IZ = (Wf.conj() * Zf).real / wsum2
                Ev = (Xf.real ** 2 + Xf.imag ** 2 + Yf.real ** 2 + Yf.imag ** 2
                      + Zf.real ** 2 + Zf.imag ** 2) / wsum2
                centers = (idx[:, 0] + nfft // 2) / fs

            for s in range(nsec_blk):
                g = sec_base + s
                seg = data[s * fs:(s + 1) * fs]
                F["rms_w"][g] = np.sqrt((seg[:, 0].astype(np.float64) ** 2).mean())
                F["peak"][g] = float(np.abs(seg).max())
                if nwin == 0:
                    continue
                sel = np.where((centers >= s) & (centers < s + 1))[0]
                if len(sel) == 0:
                    continue
                pw = Pw[sel].mean(0)
                ix, iy, iz = (IX[sel].mean(0), IY[sel].mean(0), IZ[sel].mean(0))
                ev = Ev[sel].mean(0)
                for b, bi in enumerate(oct_idx):
                    F["oct_pow"][g, b] = pw[bi].sum()
                    F["I_band"][g, b] = (ix[bi].sum(), iy[bi].sum(), iz[bi].sum())
                p = pw[spec_mask]
                F["centroid"][g] = float((freqs[spec_mask] * p).sum() / (p.sum() + eps))
                F["flatness"][g] = float(np.exp(np.log(p + eps).mean()) / (p.mean() + eps))
                np.add.at(F["logspec"][g], log_idx[log_idx >= 0], pw[log_idx >= 0])
                Ix, Iy, Iz = (ix[doa_mask].sum(), iy[doa_mask].sum(), iz[doa_mask].sum())
                F["az"][g] = np.degrees(np.arctan2(Iy, Ix))
                F["el"][g] = np.degrees(np.arctan2(Iz, np.hypot(Ix, Iy)))
                etot = (pw[doa_mask].sum() + ev[doa_mask].sum()) / 2
                inorm = float(np.sqrt(Ix ** 2 + Iy ** 2 + Iz ** 2))
                F["diffuse"][g] = 1.0 - min(1.0, inorm / (etot + eps))
                minspec[g // 60] += pw
                mincnt[g // 60] += 1

            carry = data[nsec_blk * fs:].copy()
            sec_base += nsec_blk
            if sec_base >= nsec:
                break

    minspec[mincnt > 0] /= mincnt[mincnt > 0, None]
    F["minspec"] = minspec.astype(np.float32)
    F["freqs"] = freqs.astype(np.float32)
    F["logf"] = logf.astype(np.float32)
    F["start"] = np.float64(take.start)
    F["fs"] = np.int64(fs)
    F["fast_dt"] = np.float64(FAST)
    F["hi_dt"] = np.float64(HI_ENV)
    return F

extract_session(sess, out_dir, verbose=True)

Extract features for every take; save one npz per take. Returns paths.

Source code in src/ambiscape/features.py

def extract_session(sess: Session, out_dir: str | Path, verbose=True) -> list[Path]:
    """Extract features for every take; save one npz per take. Returns paths."""
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    paths = []
    for tk in sess.takes:
        out = out_dir / (tk.path.stem + ".npz")
        if not out.exists():
            F = extract_take(tk)
            np.savez_compressed(out, **F)
            if verbose:
                print(f"  extracted {tk.path.name} ({tk.duration:.0f}s)", flush=True)
        paths.append(out)
    return paths

load_features(npz_paths)

Concatenate per-take feature files onto one absolute time axis.

Source code in src/ambiscape/features.py

def load_features(npz_paths: list[str | Path]) -> dict:
    """Concatenate per-take feature files onto one absolute time axis."""
    parts = [np.load(str(p)) for p in npz_paths]
    parts.sort(key=lambda p: float(p["start"]))
    out = {}
    out["t"] = np.concatenate([p["start"] + np.arange(len(p["rms_w"]))
                               for p in parts])
    fd = float(parts[0]["fast_dt"])
    out["t_fast"] = np.concatenate([p["start"] + fd * np.arange(len(p["fast_db"]))
                                    for p in parts])
    for k in ("fast_db", "fast_dba", "rms_w", "peak", "oct_pow", "centroid",
              "flatness", "logspec", "I_band", "az", "el", "diffuse"):
        out[k] = np.concatenate([p[k] for p in parts])
    if all("env_hi" in p for p in parts):    # absent in pre-0.2 caches
        hd = float(parts[0]["hi_dt"])
        out["hi_dt"] = hd
        out["t_hi"] = np.concatenate(
            [p["start"] + hd * np.arange(len(p["env_hi"])) for p in parts])
        out["env_hi"] = np.concatenate([p["env_hi"] for p in parts])
    out["min_t"] = np.concatenate([p["start"] + 60 * np.arange(p["minspec"].shape[0])
                                   for p in parts])
    out["minspec"] = np.concatenate([p["minspec"] for p in parts])
    out["freqs"] = parts[0]["freqs"]
    out["logf"] = parts[0]["logf"]
    return out

Descriptors, events, reverberation

Session-level descriptors, event detection, and reverberation estimation.

Descriptor conventions follow the Intercontinental-database report (2026-07-10): fast level = 125 ms RMS on W; events = fast level exceeding a running background (10th percentile in a sliding 60 s window) by >= 8 dB for

= 0.25 s; diffuseness/DOA from per-second pseudo-intensity vectors.

detect_events(fast_db, fast_dt, thresh_db=8.0, min_dur=0.25)

Return list of dicts (onset index, length, peak index, exceedance).

Source code in src/ambiscape/analysis.py

def detect_events(fast_db, fast_dt, thresh_db=8.0, min_dur=0.25):
    """Return list of dicts (onset index, length, peak index, exceedance)."""
    bg = running_background(fast_db, fast_dt)
    above = fast_db > bg + thresh_db
    events = []
    i, n = 0, len(above)
    min_len = max(1, int(round(min_dur / fast_dt)))
    while i < n:
        if above[i]:
            j = i
            while j + 1 < n and above[j + 1]:
                j += 1
            if j - i + 1 >= min_len:
                k = i + int(np.argmax(fast_db[i:j + 1]))
                events.append(dict(i0=i, i1=j, ipk=k,
                                   exceed=float(fast_db[k] - bg[k])))
            i = j + 1
        else:
            i += 1
    return events, bg

intermittency_ratio(level_db, dt, k_db=3.0)

Intermittency ratio IR (Wunderli et al. 2016), in percent.

The share of total sound energy carried by "events": frames whose level exceeds the whole-period Leq by k_db (3 dB per the original definition, there on 1 s LAeq frames — here on the fast frames, which is equivalent for events longer than the frame). IR ≈ 0 for steady scenes (drones, dense traffic), high for scenes whose energy arrives in distinct events (rail, church bells, sparse traffic).

Source code in src/ambiscape/analysis.py

def intermittency_ratio(level_db: np.ndarray, dt: float,
                        k_db: float = 3.0) -> float:
    """Intermittency ratio IR (Wunderli et al. 2016), in percent.

    The share of total sound energy carried by "events": frames whose
    level exceeds the whole-period Leq by ``k_db`` (3 dB per the original
    definition, there on 1 s LAeq frames — here on the fast frames, which
    is equivalent for events longer than the frame). IR ≈ 0 for steady
    scenes (drones, dense traffic), high for scenes whose energy arrives
    in distinct events (rail, church bells, sparse traffic).
    """
    p = 10 ** (np.asarray(level_db, np.float64) / 10)
    leq = db(p.mean())
    mask = level_db >= leq + k_db
    return float(100.0 * p[mask].sum() / (p.sum() + EPS))

decay_metrics(x, fs, bands=((250, 500), (500, 1000), (1000, 2000), (2000, 4000), (4000, 8000)))

T60, EDT, C50, C80 (dB) and D50 per octave band from an impulse.

Same truncated-Schroeder machinery as :func:decay_time (which is kept unchanged — its output feeds frozen corpus reports), plus the standard companions: EDT from the 0…−10 dB fit (perceived reverberance), clarity C50/C80 = 10·log10 of the early/late energy ratio at 50/80 ms, and definition D50 = early fraction at 50 ms. Returns {band: {"T60", "EDT", "C50", "C80", "D50", "dr_db"}}.

Source code in src/ambiscape/analysis.py

def decay_metrics(x: np.ndarray, fs: int, bands=((250, 500), (500, 1000),
                  (1000, 2000), (2000, 4000), (4000, 8000))) -> dict:
    """T60, EDT, C50, C80 (dB) and D50 per octave band from an impulse.

    Same truncated-Schroeder machinery as :func:`decay_time` (which is
    kept unchanged — its output feeds frozen corpus reports), plus the
    standard companions: EDT from the 0…−10 dB fit (perceived
    reverberance), clarity C50/C80 = 10·log10 of the early/late energy
    ratio at 50/80 ms, and definition D50 = early fraction at 50 ms.
    Returns ``{band: {"T60", "EDT", "C50", "C80", "D50", "dr_db"}}``.
    """
    from scipy import signal as sg
    pk_i = int(np.abs(x).argmax())
    env_bb = sg.convolve(x ** 2, np.ones(480) / 480, "same")
    tail = 10 * np.log10(env_bb[pk_i:pk_i + 3 * fs] + 1e-15)
    run_min = np.minimum.accumulate(tail)
    re = np.flatnonzero((tail - run_min > 8) & (np.arange(len(tail)) > fs // 10))
    cut = int(re[0]) if len(re) else 2 * fs
    out = {}
    for lo, hi in bands:
        sos = sg.butter(4, [lo, hi], "bandpass", fs=fs, output="sos")
        y = sg.sosfilt(sos, x)
        env = sg.convolve(y ** 2, np.ones(240) / 240, "same")
        pk = int(env[max(0, pk_i - 2400):pk_i + 2400].argmax()) \
            + max(0, pk_i - 2400)
        if pk < fs // 4:
            continue
        noise = float(np.median(env[:pk - fs // 8]))
        dr = 10 * np.log10(env[pk] / (noise + EPS))
        if dr < 20:
            continue
        seg = np.maximum(y[pk:pk + cut] ** 2 - noise, 0)
        sch = np.cumsum(seg[::-1])[::-1]
        sch_db = 10 * np.log10(sch / (sch[0] + EPS) + 1e-15)
        tax = np.arange(len(sch_db)) / fs
        res = {"dr_db": round(float(dr), 0)}
        for key, hi_db, lo_db in (("T60", -5.0, max(-35.0, -dr + 8)),
                                  ("EDT", 0.0, -10.0)):
            m = (sch_db <= hi_db) & (sch_db >= lo_db)
            if m.sum() < 150:
                continue
            A = np.vstack([tax[m], np.ones(int(m.sum()))]).T
            slope, _ = np.linalg.lstsq(A, sch_db[m], rcond=None)[0]
            if slope < 0:
                res[key] = round(-60.0 / slope, 2)
        for key, ms in (("C50", 50), ("C80", 80)):
            i = int(ms * fs / 1000)
            if i < len(sch) and sch[i] > 0:
                res[key] = round(float(10 * np.log10(
                    (sch[0] - sch[i]) / (sch[i] + EPS) + EPS)), 1)
        i50 = int(0.05 * fs)
        if i50 < len(sch):
            res["D50"] = round(float((sch[0] - sch[i50]) / (sch[0] + EPS)), 2)
        if "T60" in res:
            out[f"{lo}-{hi}"] = res
    return out

circular_stats(az_deg, weights=None)

Energy-weighted circular mean (deg) and resultant length R.

Source code in src/ambiscape/analysis.py

def circular_stats(az_deg, weights=None):
    """Energy-weighted circular mean (deg) and resultant length R."""
    from .circstats import mean_resultant
    mu, R = mean_resultant(np.radians(np.asarray(az_deg, float)), weights)
    return float(np.degrees(mu)), R

summarize(F)

Session descriptor dict from concatenated features (see features.load_features).

Source code in src/ambiscape/analysis.py

def summarize(F: dict) -> dict:
    """Session descriptor dict from concatenated features (see features.load_features)."""
    fast, fasta = F["fast_db"], F["fast_dba"]
    dt = float(np.median(np.diff(F["t_fast"]))) if len(F["t_fast"]) > 1 else 0.125
    leq = db(np.mean(10 ** (fast.astype(np.float64) / 10)))
    laeq = db(np.mean(10 ** (fasta.astype(np.float64) / 10)))
    l10, l50, l90 = (float(np.percentile(fast, q)) for q in (90, 50, 10))
    events, bg = detect_events(fast, dt)
    dur = float(len(F["t"]))  # 1 s per feature frame; robust across take gaps

    p = F["rms_w"].astype(np.float64) ** 2
    e_fg = p >= np.percentile(p, 75)
    e_bg = p <= np.percentile(p, 25)
    az_mean, R = circular_stats(F["az"], weights=p)
    az_fg, R_fg = circular_stats(F["az"][e_fg], weights=p[e_fg])
    el_fg = float(np.median(F["el"][e_fg]))
    psi = F["diffuse"]

    return {
        "duration_min": round(dur / 60, 1),
        "leq_dbfs": round(float(leq), 1),
        "laeq_dbfs": round(float(laeq), 1),
        "leq_minus_laeq_db": round(float(leq - laeq), 1),
        "L10": round(l10, 1), "L50": round(l50, 1), "L90": round(l90, 1),
        "dynamics_L10_L90": round(l10 - l90, 1),
        "events_per_min": round(len(events) / max(dur / 60, 1e-9), 1),
        "event_median_dur_s": round(float(np.median(
            [(e["i1"] - e["i0"] + 1) * dt for e in events])), 2) if events else None,
        "centroid_median_hz": int(np.median(F["centroid"])),
        "flatness_median": round(float(np.median(F["flatness"])), 3),
        "diffuseness_median": round(float(np.median(psi)), 2),
        "diffuseness_iqr": round(float(np.percentile(psi, 75)
                                       - np.percentile(psi, 25)), 2),
        "azimuth_mean_deg": round(az_mean, 0),
        "azimuth_R": round(R, 2),
        "azimuth_fg_deg": round(az_fg, 0),
        "elevation_fg_median_deg": round(el_fg, 0),
        "n_events": len(events),
        "emergence_db": round(float(laeq - np.percentile(fasta, 10)), 1),
        "intermittency_ratio_pct": round(intermittency_ratio(fasta, dt), 1),
    }

decay_time(x, fs, bands=((250, 500), (500, 1000), (1000, 2000), (2000, 4000), (4000, 8000)))

T60 estimates from an impulse via truncated Schroeder integration.

The decay is truncated at the first re-attack (envelope rising >= 8 dB above its running minimum) and at the noise floor; a linear fit of -5 dB .. max(-35 dB, floor + 8 dB) is extrapolated to 60 dB. Returns {band: (T60, dynamic_range_db)}.

Source code in src/ambiscape/analysis.py

def decay_time(x: np.ndarray, fs: int, bands=((250, 500), (500, 1000),
               (1000, 2000), (2000, 4000), (4000, 8000))) -> dict:
    """T60 estimates from an impulse via truncated Schroeder integration.

    The decay is truncated at the first re-attack (envelope rising >= 8 dB
    above its running minimum) and at the noise floor; a linear fit of
    -5 dB .. max(-35 dB, floor + 8 dB) is extrapolated to 60 dB.
    Returns {band: (T60, dynamic_range_db)}.
    """
    from scipy import signal as sg
    pk_i = int(np.abs(x).argmax())
    env_bb = sg.convolve(x ** 2, np.ones(480) / 480, "same")
    tail = 10 * np.log10(env_bb[pk_i:pk_i + 3 * fs] + 1e-15)
    run_min = np.minimum.accumulate(tail)
    re = np.flatnonzero((tail - run_min > 8) & (np.arange(len(tail)) > fs // 10))
    cut = int(re[0]) if len(re) else 2 * fs
    out = {}
    for lo, hi in bands:
        sos = sg.butter(4, [lo, hi], "bandpass", fs=fs, output="sos")
        y = sg.sosfilt(sos, x)
        env = sg.convolve(y ** 2, np.ones(240) / 240, "same")
        pk = int(env[max(0, pk_i - 2400):pk_i + 2400].argmax()) + max(0, pk_i - 2400)
        if pk < fs // 4:
            continue
        noise = float(np.median(env[:pk - fs // 8]))
        dr = 10 * np.log10(env[pk] / (noise + EPS))
        if dr < 20:
            continue
        seg = np.maximum(y[pk:pk + cut] ** 2 - noise, 0)
        sch = np.cumsum(seg[::-1])[::-1]
        sch_db = 10 * np.log10(sch / (sch[0] + EPS) + 1e-15)
        tax = np.arange(len(sch_db)) / fs
        lo_db = max(-35.0, -dr + 8)
        m = (sch_db <= -5) & (sch_db >= lo_db)
        if m.sum() < 150:
            continue
        A = np.vstack([tax[m], np.ones(int(m.sum()))]).T
        slope, _ = np.linalg.lstsq(A, sch_db[m], rcond=None)[0]
        if slope < 0:
            out[f"{lo}-{hi}"] = (round(-60.0 / slope, 2), round(float(dr), 0))
    return out

pick_segments(F, n=4, seg_s=600.0)

Suggest representative windows: quietest, most active, median-typical, and (if present) the strongest state transition.

Source code in src/ambiscape/analysis.py

def pick_segments(F: dict, n=4, seg_s=600.0) -> list[dict]:
    """Suggest representative windows: quietest, most active, median-typical,
    and (if present) the strongest state transition."""
    t, fast = F["t_fast"], F["fast_db"]
    dt = float(np.median(np.diff(t)))
    win = max(1, int(seg_s / dt))
    if len(fast) < win:
        return [dict(kind="whole", t0=float(t[0]), dur=float(t[-1] - t[0]))]
    k = np.ones(win) / win
    m_lvl = np.convolve(10 ** (fast.astype(np.float64) / 10), k, "valid")
    var = np.convolve((fast - fast.mean()) ** 2, k, "valid")
    picks = []
    for kind, idx in (("quietest", int(np.argmin(m_lvl))),
                      ("most_active", int(np.argmax(var))),
                      ("typical", int(np.argmin(np.abs(db(m_lvl) - np.median(db(m_lvl))))))):
        picks.append(dict(kind=kind, t0=float(t[idx]), dur=seg_s))
    smooth = median_filter(fast, size=max(3, int(30 / dt)) | 1)
    jump = np.abs(np.diff(smooth))
    if jump.max() > 6:
        picks.append(dict(kind="transition",
                          t0=float(max(t[0], t[int(np.argmax(jump))] - seg_s / 2)),
                          dur=seg_s))
    return picks[:n]

Figures

Session figures.

Names follow ambiviz conventions where the plots correspond (https://github.com/fisheggg/ambiviz): the azimuth-vs-time panel is an anglegram and the polar energy histogram a directogram, computed here from streaming per-second pseudo-intensity features rather than a full AEM, so they scale to many-hour recordings. For rich spherical maps (AEM) of short excerpts, export a segment and use ambiviz directly.

overview(F, out_path, title='', clock=None)

4-row overview: fast level + background; log spectrogram; anglegram (energy-weighted azimuth x time); diffuseness. Takes separated by more than 10 minutes get their own column (width ~ duration).

Source code in src/ambiscape/figures.py

def overview(F, out_path, title="", clock=None):
    """4-row overview: fast level + background; log spectrogram; anglegram
    (energy-weighted azimuth x time); diffuseness. Takes separated by more
    than 10 minutes get their own column (width ~ duration)."""
    with plt.rc_context(RC):
        segs = _gap_split(F["t"])
        widths = [F["t"][i1 - 1] - F["t"][i0] + 1 for i0, i1 in segs]
        fig, axes = plt.subplots(
            4, len(segs), figsize=(12, 9.5), dpi=130, sharey="row",
            gridspec_kw={"width_ratios": widths, "wspace": 0.04}, squeeze=False)

        tf_all, fast_all = F["t_fast"], F["fast_db"]
        dt = float(np.median(np.diff(tf_all))) if len(tf_all) > 1 else 0.125
        n_events_total = 0
        nonempty = F["logspec"].sum(0) > 0
        Sall = db(F["logspec"][:, nonempty])
        vmax = np.percentile(Sall, 99.5)
        fc = np.sqrt(F["logf"][:-1] * F["logf"][1:])[nonempty]
        p_all = F["rms_w"].astype(np.float64) ** 2

        for col, (i0, i1) in enumerate(segs):
            t = F["t"][i0:i1]
            fm = (tf_all >= t[0]) & (tf_all < t[-1] + 1)
            tf, fast = tf_all[fm], fast_all[fm]
            events, bg = detect_events(fast, dt)
            n_events_total += len(events)

            ax = axes[0][col]
            ax.plot(tf, fast, color=BLUE, lw=0.4, alpha=0.7)
            ax.plot(tf, bg, color=YELLOW, lw=1.4)
            for q, ls in ((90, ":"), (50, "-"), (10, ":")):
                ax.axhline(np.percentile(fast_all, q), color=MUT, lw=0.7, ls=ls)
            ax.set_xlim(t[0], t[-1])

            ax = axes[1][col]
            ax.pcolormesh(t, fc, Sall[i0:i1].T, cmap=SEQ, vmin=vmax - 65,
                          vmax=vmax, shading="auto", rasterized=True)
            ax.set_yscale("log")
            ax.set_ylim(25, 16000)
            ax.grid(False)

            ax = axes[2][col]
            nbins = 36
            nb_t = max(int((t[-1] - t[0]) / 30.0), 1)
            tb = np.linspace(t[0], t[-1] + 1, nb_t + 1)
            H = np.zeros((nbins, nb_t))
            azb = np.linspace(-180, 180, nbins + 1)
            ti = np.clip(np.searchsorted(tb, t) - 1, 0, nb_t - 1)
            ai = np.clip(np.searchsorted(azb, F["az"][i0:i1]) - 1, 0, nbins - 1)
            np.add.at(H, (ai, ti), p_all[i0:i1])
            ax.pcolormesh(tb[:-1], azb[:-1], db(H), cmap=SEQ, shading="auto",
                          rasterized=True, vmin=db(H).max() - 40,
                          vmax=db(H).max())
            ax.set_yticks([-180, -90, 0, 90, 180])
            ax.grid(False)

            ax = axes[3][col]
            d = F["diffuse"][i0:i1]
            ax.plot(t, d, color=GREEN, lw=0.5, alpha=0.55)
            k = min(121, max(3, len(t) // 50) | 1)
            ax.plot(t, np.convolve(d, np.ones(k) / k, "same"),
                    color=GREEN, lw=1.6)
            ax.set_ylim(0, 1)
            for row in range(4):
                _time_axis(axes[row][col], t, clock)
                if row < 3:
                    axes[row][col].tick_params(labelbottom=False)
            if clock and len(segs) > 1:
                axes[0][col].set_title(clock(t[0])[:6], loc="left",
                                       fontsize=8.5, color=SEC)

        axes[0][0].set_ylabel("fast level (dBFS)")
        axes[1][0].set_ylabel("frequency (Hz)")
        axes[2][0].set_ylabel("azimuth (°)\nanglegram")
        axes[3][0].set_ylabel("diffuseness ψ")
        fig.suptitle(f"{title} — level (blue), running background (yellow), "
                     f"L10/L50/L90 (grey); {n_events_total} events",
                     x=0.01, ha="left", fontsize=10, color=INK)
        fig.tight_layout(rect=(0, 0, 1, 0.985))
        fig.savefig(out_path, bbox_inches="tight")
        plt.close(fig)

ltas_percentiles(F, out_path, title='')

10/50/90th percentile long-term spectra (background vs foreground).

Source code in src/ambiscape/figures.py

def ltas_percentiles(F, out_path, title=""):
    """10/50/90th percentile long-term spectra (background vs foreground)."""
    with plt.rc_context(RC):
        nonempty = F["logspec"].sum(0) > 0
        S = db(F["logspec"][:, nonempty])
        fc = np.sqrt(F["logf"][:-1] * F["logf"][1:])[nonempty]
        fig, ax = plt.subplots(figsize=(8, 4), dpi=130)
        p10, p50, p90 = (np.percentile(S, q, axis=0) for q in (10, 50, 90))
        ax.fill_between(fc, p10, p90, color=BLUE, alpha=0.18, lw=0)
        ax.plot(fc, p50, color=BLUE, lw=1.5)
        ax.plot(fc, p10, color=MUT, lw=0.8)
        ax.plot(fc, p90, color=MAGENTA, lw=1.0)
        ax.annotate("90th pct (foreground)", (fc[-30], p90[-30]), color=MAGENTA,
                    fontsize=8, xytext=(0, 6), textcoords="offset points")
        ax.annotate("10th pct (background)", (fc[-30], p10[-30]), color=MUT,
                    fontsize=8, xytext=(0, -12), textcoords="offset points")
        ax.set_xscale("log")
        ax.set_xlabel("frequency (Hz)")
        ax.set_ylabel("PSD (dB)")
        ax.set_title(f"{title} — percentile LTAS", loc="left", fontsize=10)
        fig.tight_layout()
        fig.savefig(out_path, bbox_inches="tight")
        plt.close(fig)

directogram(F, out_path, title='')

Polar azimuth histograms: foreground (loudest 25 %) vs background (quietest 25 %) energy-weighted. ambiviz-style directogram.

Source code in src/ambiscape/figures.py

def directogram(F, out_path, title=""):
    """Polar azimuth histograms: foreground (loudest 25 %) vs background
    (quietest 25 %) energy-weighted. ambiviz-style directogram."""
    with plt.rc_context(RC):
        p = F["rms_w"].astype(np.float64) ** 2
        fg, bgm = p >= np.percentile(p, 75), p <= np.percentile(p, 25)
        fig, ax = plt.subplots(figsize=(5, 5), dpi=130,
                               subplot_kw={"projection": "polar"})
        bins = np.linspace(-180, 180, 37)
        th = np.radians((bins[:-1] + bins[1:]) / 2)
        for m, c, lab in ((bgm, MUT, "background (quietest 25%)"),
                          (fg, BLUE, "foreground (loudest 25%)")):
            h, _ = np.histogram(F["az"][m], bins=bins, weights=p[m])
            ax.bar(th, h / (h.max() + 1e-20), width=np.radians(10) * 0.92,
                   color=c, alpha=0.65, label=lab)
        ax.set_theta_zero_location("N")
        ax.set_theta_direction(1)
        ax.set_thetagrids([0, 90, 180, 270],
                          ["front", "left", "rear", "right"], fontsize=8.5)
        ax.set_rticks([])
        ax.set_title(f"{title} — directogram", fontsize=10, pad=14)
        ax.legend(loc="lower left", bbox_to_anchor=(-0.1, -0.12),
                  frameon=False, fontsize=8)
        fig.tight_layout()
        fig.savefig(out_path, bbox_inches="tight")
        plt.close(fig)

Spectral foreground

Per-band running background and spectral foreground decomposition.

The broadband event detector in :mod:analysis misses band-limited events riding on a loud bed in other bands (distant bells over traffic move their octave a few dB while the broadband level barely changes). This module works on the cached 1 Hz log-band spectrogram (F["logspec"], 96 bands):

• band_background — running low-percentile background per band;
• foreground — dB exceedance and the per-second foreground fraction (share of total power sitting above the spectral background);
• spectral_events — connected spectro-temporal regions of exceedance (time x band blobs), each with onset, duration, band span, and peak rise;
• summarize_foreground — session descriptors appended to the analyze summary and README.

All functions are pure array transforms on cached features — no audio pass.

band_background(logspec, win_s=300.0, pct=10.0)

Running pct-percentile background per log band.

logspec is the (nsec, nband) power array from the cached features; the window is in seconds (= rows). Returns the same shape.

Source code in src/ambiscape/background.py

def band_background(logspec: np.ndarray, win_s: float = 300.0,
                    pct: float = 10.0) -> np.ndarray:
    """Running ``pct``-percentile background per log band.

    ``logspec`` is the (nsec, nband) power array from the cached features;
    the window is in seconds (= rows). Returns the same shape.
    """
    n = max(3, int(round(win_s)) | 1)
    return ndimage.percentile_filter(logspec, pct, size=(n, 1),
                                     mode="nearest")

foreground(logspec, bg)

dB rise above the spectral background and per-second foreground fraction (share of total power more than 3 dB above background).

Source code in src/ambiscape/background.py

def foreground(logspec: np.ndarray, bg: np.ndarray):
    """dB rise above the spectral background and per-second foreground
    fraction (share of total power more than 3 dB above background)."""
    rise_db = 10 * np.log10((logspec + EPS) / (bg + EPS))
    fg_mask = rise_db > 3.0
    frac = (logspec * fg_mask).sum(1) / (logspec.sum(1) + EPS)
    return rise_db, frac

spectral_events(rise_db, logf, thresh_db=6.0, min_dur_s=2.0, min_bands=2)

Connected regions of band-wise exceedance as event dicts.

A spectral event is a blob in the (time x band) plane where the rise exceeds thresh_db, lasting >= min_dur_s and spanning >= min_bands bands. Returns onset/duration (s), band span (Hz), and peak rise (dB), sorted by onset.

Source code in src/ambiscape/background.py

def spectral_events(rise_db: np.ndarray, logf: np.ndarray,
                    thresh_db: float = 6.0, min_dur_s: float = 2.0,
                    min_bands: int = 2) -> list[dict]:
    """Connected regions of band-wise exceedance as event dicts.

    A spectral event is a blob in the (time x band) plane where the rise
    exceeds ``thresh_db``, lasting >= ``min_dur_s`` and spanning >=
    ``min_bands`` bands. Returns onset/duration (s), band span (Hz), and
    peak rise (dB), sorted by onset.
    """
    lab, nlab = ndimage.label(rise_db > thresh_db)
    out = []
    for sl_t, sl_b in ndimage.find_objects(lab):
        dur = sl_t.stop - sl_t.start
        nb = sl_b.stop - sl_b.start
        if dur < min_dur_s or nb < min_bands:
            continue
        blob = rise_db[sl_t, sl_b]
        out.append({
            "t0_s": int(sl_t.start),
            "dur_s": int(dur),
            "f_lo_hz": round(float(logf[sl_b.start]), 1),
            "f_hi_hz": round(float(logf[min(sl_b.stop, len(logf) - 1)]), 1),
            "peak_rise_db": round(float(blob.max()), 1),
        })
    return sorted(out, key=lambda e: e["t0_s"])

masking_index(F, active, quiet)

How much a dominant source hides the rest of the field — the "lo-fi" claim as a number.

active/quiet are boolean second-masks (source on / off). Per log band, the floor elevation is the rise of the active-state median level above the quiet-state median: ambient sounds in that band must now exceed the elevated typical floor to be audible. Returns the median and maximum elevation over 250 Hz–8 kHz, the fraction of bands elevated by more than 6 dB, and the per-band curve.

Source code in src/ambiscape/background.py

def masking_index(F: dict, active: np.ndarray, quiet: np.ndarray) -> dict:
    """How much a dominant source hides the rest of the field — the "lo-fi"
    claim as a number.

    ``active``/``quiet`` are boolean second-masks (source on / off). Per log
    band, the floor elevation is the rise of the active-state median level
    above the quiet-state median: ambient sounds in that band must now
    exceed the elevated typical floor to be audible. Returns the median and
    maximum elevation over 250 Hz–8 kHz, the fraction of bands elevated by
    more than 6 dB, and the per-band curve.
    """
    ls = F["logspec"]
    logf = F["logf"]
    el = 10 * np.log10(np.median(ls[active], axis=0) + EPS) \
        - 10 * np.log10(np.median(ls[quiet], axis=0) + EPS)
    band = (logf[:-1] >= 250) & (logf[:-1] <= 8000)
    return {
        "floor_elevation_median_db": round(float(np.median(el[band])), 1),
        "floor_elevation_max_db": round(float(el[band].max()), 1),
        "bands_masked_gt6db_fraction": round(float((el[band] > 6).mean()), 2),
        "elevation_db_per_band": [round(float(v), 1) for v in el],
    }

source_fingerprint(F, active, quiet, fmin=25.0, fmax=16000.0, min_prom_db=6.0, max_peaks=20)

Spectral fingerprint of a source: active-minus-quiet mean PSDs.

active/quiet are boolean masks over the minutes of F["minspec"] (source clearly on / clearly off, e.g. from :func:ambiscape.states.state_segments). The rise curve is the dB difference of the two mean spectra — the source's own spectrum with the room ambience subtracted. Narrowband peaks of the rise are extracted and passed through the harmonic sieve, so a blade-pass or compressor comb reports its base frequency.

Returns dict: freqs/rise_db (the full curve), rise_max_db/ rise_max_hz (the turbulence hump), peaks (list of {f_hz, rise_db}), and comb ({f0_hz, harmonicity} of the peak set, f0_hz None when there are no peaks).

Source code in src/ambiscape/background.py

def source_fingerprint(F: dict, active: np.ndarray, quiet: np.ndarray,
                       fmin: float = 25.0, fmax: float = 16000.0,
                       min_prom_db: float = 6.0, max_peaks: int = 20) -> dict:
    """Spectral fingerprint of a source: active-minus-quiet mean PSDs.

    ``active``/``quiet`` are boolean masks over the *minutes* of
    ``F["minspec"]`` (source clearly on / clearly off, e.g. from
    :func:`ambiscape.states.state_segments`). The rise curve is the dB
    difference of the two mean spectra — the source's own spectrum with the
    room ambience subtracted. Narrowband peaks of the rise are extracted and
    passed through the harmonic sieve, so a blade-pass or compressor comb
    reports its base frequency.

    Returns dict: ``freqs``/``rise_db`` (the full curve), ``rise_max_db``/
    ``rise_max_hz`` (the turbulence hump), ``peaks`` (list of
    ``{f_hz, rise_db}``), and ``comb`` (``{f0_hz, harmonicity}`` of the peak
    set, ``f0_hz`` None when there are no peaks).
    """
    from scipy.ndimage import median_filter as _medf
    from scipy.signal import find_peaks as _find_peaks
    from .tonality import harmonic_sieve

    freqs = np.asarray(F["freqs"], float)
    S_a = F["minspec"][np.asarray(active, bool)].mean(0)
    S_q = F["minspec"][np.asarray(quiet, bool)].mean(0)
    m = (freqs >= fmin) & (freqs <= fmax)
    rise = 10 * np.log10((S_a[m] + EPS) / (S_q[m] + EPS))
    fsel = freqs[m]

    # hump: broad maximum of the smoothed rise
    smooth = _medf(rise, size=51, mode="nearest")
    i_max = int(np.argmax(smooth))
    # peaks: narrowband lines above the smoothed curve
    line = rise - smooth
    pk, props = _find_peaks(line, height=min_prom_db, distance=5)
    order = np.argsort(props["peak_heights"])[::-1][:max_peaks]
    keep = np.sort(pk[order])
    peaks = [{"f_hz": round(float(fsel[i]), 1),
              "rise_db": round(float(rise[i]), 1)} for i in keep]
    f0, h = harmonic_sieve(fsel[keep], 10 ** (rise[keep] / 10)) \
        if len(keep) else (None, 0.0)
    return {
        "freqs": fsel, "rise_db": rise,
        "rise_max_db": round(float(smooth[i_max]), 1),
        "rise_max_hz": round(float(fsel[i_max]), 1),
        "peaks": peaks,
        "comb": {"f0_hz": round(f0, 1) if f0 else None, "harmonicity": h},
    }

summarize_foreground(F, win_s=300.0)

Foreground descriptors for the analyze summary.

Source code in src/ambiscape/background.py

def summarize_foreground(F: dict, win_s: float = 300.0) -> dict:
    """Foreground descriptors for the analyze summary."""
    bg = band_background(F["logspec"], win_s=win_s)
    rise_db, frac = foreground(F["logspec"], bg)
    ev = spectral_events(rise_db, F["logf"])
    dur_min = max(len(frac) / 60.0, 1e-9)
    return {
        "fg_fraction_median": round(float(np.median(frac)), 2),
        "fg_fraction_p90": round(float(np.percentile(frac, 90)), 2),
        "spectral_events_per_min": round(len(ev) / dur_min, 1),
        "spectral_event_median_dur_s": (
            round(float(np.median([e["dur_s"] for e in ev])), 1)
            if ev else None),
    }

Machine states

Machine states: on/off segmentation, switch points, and duty cycles.

Domestic and mechanical sources (ventilation, fridges, pumps, HVAC) show up in a soundscape as a state — a band-limited floor that is either present or absent — rather than as events. This module segments a band-level timeline into those states from the cached features, no audio pass:

• band_level — per-second dB level in a frequency band from the cached log-band spectrogram (the "machine band" of a source, e.g. 250–1000 Hz for a ventilation unit);
• state_segments — two-state (on/off) segmentation of that level with an automatic bimodal threshold, hysteresis, and a minimum duration, each segment carrying its median level and within-state stability (SD);
• switch_points — the transitions between segments (the 07:53:55 switch-off moments);
• duty_cycle — cycle statistics of a cycling machine (a fridge's ~24 min period at ~50 % duty): period, duty fraction, cycle count.

Typical use: segs = state_segments(band_level(F, (250, 1000))) and mask other analyses (fingerprints, masking, taxonomy states) by segment.

band_level(F, band=(250.0, 1000.0))

Per-second dB level in band (Hz) from the cached logspec.

Source code in src/ambiscape/states.py

def band_level(F: dict, band=(250.0, 1000.0)) -> np.ndarray:
    """Per-second dB level in ``band`` (Hz) from the cached ``logspec``."""
    logf = np.asarray(F["logf"], float)
    fc = np.sqrt(logf[:-1] * logf[1:])
    m = (fc >= band[0]) & (fc <= band[1])
    return 10 * np.log10(F["logspec"][:, m].sum(1) + EPS)

bimodal_threshold(level_db)

Otsu's threshold on the level histogram: the split that best separates the two modes of an on/off timeline.

Source code in src/ambiscape/states.py

def bimodal_threshold(level_db: np.ndarray) -> float:
    """Otsu's threshold on the level histogram: the split that best
    separates the two modes of an on/off timeline."""
    lo, hi = np.percentile(level_db, (0.5, 99.5))
    hist, edges = np.histogram(level_db, bins=128, range=(lo, hi))
    p = hist.astype(float) / max(hist.sum(), 1)
    centers = (edges[:-1] + edges[1:]) / 2
    w0 = np.cumsum(p)
    mu = np.cumsum(p * centers)
    mu_t = mu[-1]
    var = (mu_t * w0 - mu) ** 2 / (w0 * (1 - w0) + EPS)
    k = int(np.argmax(var))
    # between-class variance is flat across an empty inter-mode gap; the
    # midpoint of the two class means splits the gap centrally
    mu0 = mu[k] / (w0[k] + EPS)
    mu1 = (mu_t - mu[k]) / (1 - w0[k] + EPS)
    return float((mu0 + mu1) / 2)

state_segments(level_db, thresh_db=None, smooth_s=11.0, hysteresis_db=1.0, min_dur_s=30.0)

Two-state segmentation of a 1 Hz band-level timeline.

The level is median-smoothed over smooth_s; the threshold defaults to the bimodal (Otsu) split of the histogram — pass thresh_db when the timeline is not clearly bimodal. Hysteresis of hysteresis_db around the threshold suppresses chatter, and segments shorter than min_dur_s are merged into their neighbors. Returns segments in time order as dicts: state ('on'/'off'), t0_s/dur_s (seconds into the timeline), median_db, and sd_db (within-state stability of the raw level — a running machine is steady, ambience is not).

Source code in src/ambiscape/states.py

def state_segments(level_db: np.ndarray, thresh_db: float | None = None,
                   smooth_s: float = 11.0, hysteresis_db: float = 1.0,
                   min_dur_s: float = 30.0) -> list[dict]:
    """Two-state segmentation of a 1 Hz band-level timeline.

    The level is median-smoothed over ``smooth_s``; the threshold defaults
    to the bimodal (Otsu) split of the histogram — pass ``thresh_db`` when
    the timeline is not clearly bimodal. Hysteresis of ``hysteresis_db``
    around the threshold suppresses chatter, and segments shorter than
    ``min_dur_s`` are merged into their neighbors. Returns segments in time
    order as dicts: state ('on'/'off'), t0_s/dur_s (seconds into the
    timeline), median_db, and sd_db (within-state stability of the raw
    level — a running machine is *steady*, ambience is not).
    """
    x = np.asarray(level_db, float)
    k = max(3, int(round(smooth_s)) | 1)
    sm = median_filter(x, size=k, mode="nearest")
    th = bimodal_threshold(sm) if thresh_db is None else float(thresh_db)
    on = np.zeros(len(sm), bool)
    cur = sm[0] > th
    for i, v in enumerate(sm):
        if cur and v < th - hysteresis_db / 2:
            cur = False
        elif not cur and v > th + hysteresis_db / 2:
            cur = True
        on[i] = cur

    def bounds(mask):
        edges = np.flatnonzero(np.diff(mask.astype(int))) + 1
        return [0, *edges.tolist(), len(mask)]

    # merge runs shorter than min_dur_s into the surrounding state
    b = bounds(on)
    for i0, i1 in zip(b[:-1], b[1:]):
        if i1 - i0 < min_dur_s and i0 > 0 and i1 < len(on):
            on[i0:i1] = on[i0 - 1]
    b = bounds(on)
    segs = []
    for i0, i1 in zip(b[:-1], b[1:]):
        seg = x[i0:i1]
        segs.append({
            "state": "on" if on[i0] else "off",
            "t0_s": float(i0), "dur_s": float(i1 - i0),
            "median_db": round(float(np.median(seg)), 1),
            "sd_db": round(float(seg.std()), 2),
        })
    return segs

switch_points(segments)

Transitions between consecutive segments: time and direction ('on' = machine starts, 'off' = machine stops).

Source code in src/ambiscape/states.py

def switch_points(segments: list[dict]) -> list[dict]:
    """Transitions between consecutive segments: time and direction
    ('on' = machine starts, 'off' = machine stops)."""
    out = []
    for a, b in zip(segments[:-1], segments[1:]):
        out.append({"t_s": float(b["t0_s"]),
                    "direction": b["state"],
                    "step_db": round(b["median_db"] - a["median_db"], 1)})
    return out

duty_cycle(segments)

Cycle statistics of a cycling machine from its state segments: median period (consecutive on-starts), duty fraction (median on-time over period), and the number of complete cycles observed.

Source code in src/ambiscape/states.py

def duty_cycle(segments: list[dict]) -> dict:
    """Cycle statistics of a cycling machine from its state segments:
    median period (consecutive on-starts), duty fraction (median on-time
    over period), and the number of complete cycles observed."""
    on_starts = np.array([s["t0_s"] for s in segments if s["state"] == "on"])
    on_durs = np.array([s["dur_s"] for s in segments if s["state"] == "on"])
    if len(on_starts) < 2:
        return {"period_s": None, "duty": None,
                "n_cycles": int(len(on_starts))}
    period = float(np.median(np.diff(on_starts)))
    return {"period_s": round(period, 1),
            "duty": round(float(np.median(on_durs)) / period, 3),
            "n_cycles": int(len(on_starts))}

Electric network frequency (ENF)

Electric network frequency (ENF) traces from mains hum.

Buildings hum at the mains frequency and its harmonics (50 Hz nominal in Europe; magnetostriction is strongest at 100 Hz), and the grid's actual frequency wanders by tens of millihertz as load and generation balance. A long indoor recording therefore carries a continuous, involuntary log of the grid — usable as a session descriptor (how electrified is this room?), as a source separator (a "50 Hz" line that does not follow the grid is a rotor, not electricity), and forensically: matched against published grid-frequency archives, an ENF trace timestamps a recording independently of the recorder clock.

• hum_peak — sub-millihertz line frequency in one mono window (zero-padded FFT + parabolic interpolation) with its rise over the local spectral floor;
• enf_track — the trace: windows every step_s across a whole session, one or more harmonics, all scaled to the fundamental;
• enf_summary — mean/SD/max deviation, coverage, and cross-harmonic agreement — the latter is the authenticity check (independent acoustic lines reporting the same electrical frequency).

Needs raw audio (one streaming pass over the W channel); the cached per-minute spectra are far too coarse (5.9 Hz bins) for millihertz work.

hum_peak(w, fs, nominal=50.0, search_hz=0.2, nfft_mult=4)

Frequency and floor-rise of the strongest line near nominal.

Zero-padded FFT of the Hann-windowed mono signal, parabolic interpolation of the log-power peak within nominal ± search_hz. Returns (freq_hz, rise_db); rise is measured against the median power in a ±1.5 Hz-widened neighborhood, so a genuine line scores high even on a rumble shoulder.

Source code in src/ambiscape/enf.py

def hum_peak(w: np.ndarray, fs: int, nominal: float = 50.0,
             search_hz: float = 0.2, nfft_mult: int = 4):
    """Frequency and floor-rise of the strongest line near ``nominal``.

    Zero-padded FFT of the Hann-windowed mono signal, parabolic
    interpolation of the log-power peak within ``nominal ± search_hz``.
    Returns ``(freq_hz, rise_db)``; rise is measured against the median
    power in a ±1.5 Hz-widened neighborhood, so a genuine line scores
    high even on a rumble shoulder.
    """
    n = len(w)
    W = np.fft.rfft(w * np.hanning(n), n * nfft_mult)
    f = np.fft.rfftfreq(n * nfft_mult, 1 / fs)
    P = W.real ** 2 + W.imag ** 2
    m = (f >= nominal - search_hz) & (f <= nominal + search_hz)
    j = int(np.flatnonzero(m)[0] + np.argmax(P[m]))
    a, b, c = (np.log(P[j - 1] + EPS), np.log(P[j] + EPS),
               np.log(P[j + 1] + EPS))
    d = 0.5 * (a - c) / (a - 2 * b + c + EPS)
    floor = np.median(P[(f >= nominal - search_hz - 1.5)
                        & (f <= nominal + search_hz + 1.5)])
    return float(f[j] + d * (f[1] - f[0])), \
        float(10 * np.log10(P[j] / (floor + EPS)))

enf_track(sess, step_s=300.0, win_s=60.0, nominal=50.0, search_hz=0.2, harmonics=(1, 2), channel=0)

Track the mains hum across a whole session.

One window of win_s every step_s, per take (windows start 1 s into each take and reads shorter than 90 % of the window are skipped — recorder 2 GB splits overlap by a fraction of a second, so a read at an exact take start can return a sliver of the previous file). Each harmonic k is searched at k*nominal ± k*search_hz and reported scaled to the fundamental.

Returns {"t": absolute seconds, "f": {k: freq_hz/k}, "rise": {k: rise_db}}.

Source code in src/ambiscape/enf.py

def enf_track(sess: Session, step_s: float = 300.0, win_s: float = 60.0,
              nominal: float = 50.0, search_hz: float = 0.2,
              harmonics=(1, 2), channel: int = 0) -> dict:
    """Track the mains hum across a whole session.

    One window of ``win_s`` every ``step_s``, per take (windows start 1 s
    into each take and reads shorter than 90 % of the window are skipped —
    recorder 2 GB splits overlap by a fraction of a second, so a read at an
    exact take start can return a sliver of the previous file). Each
    harmonic ``k`` is searched at ``k*nominal ± k*search_hz`` and reported
    scaled to the fundamental.

    Returns ``{"t": absolute seconds, "f": {k: freq_hz/k}, "rise":
    {k: rise_db}}``.
    """
    ts, fk, rk = [], {k: [] for k in harmonics}, {k: [] for k in harmonics}
    for tk in sess.takes:
        t = tk.start + 1.0
        while t + win_s <= tk.end:
            x, fs = read_span(sess, t, win_s)
            if x.shape[0] >= 0.9 * win_s * fs:
                w = x[:, channel].astype(np.float64)
                for k in harmonics:
                    f, r = hum_peak(w, fs, nominal=k * nominal,
                                    search_hz=k * search_hz)
                    fk[k].append(f / k)
                    rk[k].append(r)
                ts.append(t)
            t += step_s
    return {"t": np.array(ts),
            "f": {k: np.array(v) for k, v in fk.items()},
            "rise": {k: np.array(v) for k, v in rk.items()}}

enf_summary(track, nominal=50.0, min_rise_db=6.0)

Descriptors of an ENF trace.

Statistics use only windows where the first harmonic rises min_rise_db above the floor; coverage is the fraction of windows that qualify. harmonic_agreement_mhz is the median absolute difference between the first two tracked harmonics (fundamental-scaled) where both are detected — millihertz-level agreement authenticates the line as electrical.

Source code in src/ambiscape/enf.py

def enf_summary(track: dict, nominal: float = 50.0,
                min_rise_db: float = 6.0) -> dict:
    """Descriptors of an ENF trace.

    Statistics use only windows where the first harmonic rises
    ``min_rise_db`` above the floor; ``coverage`` is the fraction of
    windows that qualify. ``harmonic_agreement_mhz`` is the median absolute
    difference between the first two tracked harmonics (fundamental-scaled)
    where both are detected — millihertz-level agreement authenticates the
    line as electrical.
    """
    ks = sorted(track["f"])
    k0 = ks[0]
    good = track["rise"][k0] >= min_rise_db
    f = track["f"][k0][good]
    out = {
        "n_windows": int(len(track["t"])),
        "coverage": round(float(good.mean()), 2) if len(good) else 0.0,
        "mean_hz": round(float(f.mean()), 4) if len(f) else None,
        "sd_mhz": round(float(f.std() * 1000), 1) if len(f) else None,
        "max_dev_mhz": round(float(np.abs(f - nominal).max() * 1000), 1)
        if len(f) else None,
        "median_rise_db": round(float(np.median(track["rise"][k0])), 1)
        if len(good) else None,
    }
    if len(ks) > 1:
        k1 = ks[1]
        both = good & (track["rise"][k1] >= min_rise_db)
        if both.any():
            d = np.abs(track["f"][k0][both] - track["f"][k1][both])
            out["harmonic_agreement_mhz"] = round(float(np.median(d) * 1000),
                                                  2)
    return out

Ecoacoustic indices

Ecoacoustic indices from the cached log-band spectrogram.

The standard soundscape-ecology battery, so sessions are reportable in the idiom global acoustic-monitoring corpora expect. All computed from the cached 1 Hz features — no audio pass:

• ACI (acoustic complexity, Pieretti et al. 2011): per-band temporal variation |ΔP|/ΣP summed over bands, averaged over 5-min chunks — sensitive to biophonic modulation, blind to steady drones;
• ADI / AEI (diversity / evenness, Villanueva-Rivera et al. 2011): Shannon entropy / Gini coefficient of the occupancy of 1 kHz bins (fraction of cells above a threshold re the band maximum);
• NDSI (Kasten et al. 2012): (biophony − anthropophony) / (biophony + anthropophony) with the conventional bands 2–8 kHz vs 1–2 kHz, in [−1, 1];
• BI (bioacoustic index, Boelman et al. 2007): area of the mean 2–8 kHz dB spectrum above its minimum;
• acoustic entropy H (Sueur et al. 2008): spectral entropy × temporal (envelope) entropy, in [0, 1].

Caveats for an ambisonic indoor corpus: these indices were designed for outdoor terrestrial monitoring; report them for comparability, but read NDSI/BI as "energy in the bird band", not as proof of birds — a 4 kHz ventilation hiss scores as "biophony". Combine with the taxonomy layer before interpreting.

aci(F, chunk_s=300.0)

Acoustic complexity index, mean over chunk_s chunks.

Source code in src/ambiscape/ecology.py

def aci(F: dict, chunk_s: float = 300.0) -> float:
    """Acoustic complexity index, mean over ``chunk_s`` chunks."""
    S = np.asarray(F["logspec"], float)
    n = max(2, int(chunk_s))
    vals = []
    for i0 in range(0, S.shape[0] - n + 1, n):
        c = S[i0:i0 + n]
        vals.append(float((np.abs(np.diff(c, axis=0)).sum(0)
                           / (c.sum(0) + EPS)).sum()))
    return float(np.mean(vals)) if vals else 0.0

adi_aei(F, **kw)

Acoustic diversity (Shannon, normalized) and evenness (Gini).

Source code in src/ambiscape/ecology.py

def adi_aei(F: dict, **kw):
    """Acoustic diversity (Shannon, normalized) and evenness (Gini)."""
    occ = _occupancy(F, **kw)
    p = occ / (occ.sum() + EPS)
    adi = float(-(p * np.log(p + EPS)).sum() / np.log(len(p) + EPS))
    x = np.sort(occ)
    n = len(x)
    gini = float((2 * np.arange(1, n + 1) - n - 1).dot(x)
                 / (n * x.sum() + EPS))
    return adi, gini

ndsi(F, anthro=(1000.0, 2000.0), bio=(2000.0, 8000.0))

Normalized difference soundscape index in [−1, 1].

Source code in src/ambiscape/ecology.py

def ndsi(F: dict, anthro=(1000.0, 2000.0), bio=(2000.0, 8000.0)) -> float:
    """Normalized difference soundscape index in [−1, 1]."""
    fc = _band_centers(F["logf"])
    S = np.asarray(F["logspec"], float).mean(0)
    a = S[(fc >= anthro[0]) & (fc < anthro[1])].sum()
    b = S[(fc >= bio[0]) & (fc < bio[1])].sum()
    return float((b - a) / (b + a + EPS))

bioacoustic_index(F, band=(2000.0, 8000.0))

Boelman BI: area of the mean band dB spectrum above its minimum.

Source code in src/ambiscape/ecology.py

def bioacoustic_index(F: dict, band=(2000.0, 8000.0)) -> float:
    """Boelman BI: area of the mean band dB spectrum above its minimum."""
    fc = _band_centers(F["logf"])
    m = (fc >= band[0]) & (fc <= band[1])
    s = 10 * np.log10(np.asarray(F["logspec"], float).mean(0)[m] + EPS)
    return float((s - s.min()).sum())

acoustic_entropy(F)

Sueur H = spectral entropy × temporal entropy, in [0, 1].

Source code in src/ambiscape/ecology.py

def acoustic_entropy(F: dict) -> float:
    """Sueur H = spectral entropy × temporal entropy, in [0, 1]."""
    S = np.asarray(F["logspec"], float)
    ps = S.mean(0)
    ps = ps / (ps.sum() + EPS)
    hf = float(-(ps * np.log(ps + EPS)).sum() / np.log(len(ps)))
    env = np.asarray(F["rms_w"], float)
    pe = env / (env.sum() + EPS)
    ht = float(-(pe * np.log(pe + EPS)).sum() / np.log(len(pe)))
    return hf * ht

indices(F)

The full battery as one dict.

Source code in src/ambiscape/ecology.py

def indices(F: dict) -> dict:
    """The full battery as one dict."""
    adi_, aei_ = adi_aei(F)
    return {
        "aci": round(aci(F), 1),
        "adi": round(adi_, 3),
        "aei": round(aei_, 3),
        "ndsi": round(ndsi(F), 3),
        "bi": round(bioacoustic_index(F), 1),
        "acoustic_entropy": round(acoustic_entropy(F), 3),
    }

summarize_ecology(F)

Alias of :func:indices for the analyze-summary pipeline.

Source code in src/ambiscape/ecology.py

def summarize_ecology(F: dict) -> dict:
    """Alias of :func:`indices` for the analyze-summary pipeline."""
    return indices(F)

Biophony

Biophony measures: capturing nature and animal sounds by structure.

The ecoacoustic battery in :mod:ambiscape.ecology (NDSI, BI, ADI) reads energy in a band; it cannot tell a dawn chorus from a ventilation hiss. Biophony is distinguished by how it is structured — narrowband, tonal, bursty in time, and (in an ambisonic recording) arriving from many elevated bearings at once. This module measures that structure from the cached features, no audio pass:

• narrowband_activity — persistent narrow spectral peaks in the bird band per minute (from the per-minute high-resolution PSD): birdsong is narrowband and tonal, wind and machines are broadband;
• band_temporal_entropy — Sueur temporal entropy of the bird-band envelope: structured vocalization concentrates energy in time (low Ht), a steady noise floor spreads it evenly (Ht → 1);
• band_activity — fraction of seconds and event rate where the bird band rises above its own running background (Towsey-style acoustic activity), restricted to the biophony band;
• spatial_dispersion — the ambisonic layer no other corpus tool has: the directional entropy and above-horizon energy fraction of the bird-band foreground — a chorus of many birds from many elevated directions is unmistakable, and it cross-checks a suspicious NDSI.

summarize_biophony returns the descriptor set for the analyze summary.

Caveats: these are acoustic-structure proxies, not detections. A tonal alarm, a whistling kettle, or a squealing fan belt can mimic biophonic structure; confirm species with the BirdNET hook (:func:ambiscape.ml.birdnet_session, [ml] extra) on the hi-fi windows. The default band (2–11 kHz) targets temperate birdsong; widen it (insects reach 8–16 kHz, many mammals sit below 2 kHz) per habitat.

narrowband_activity(F, band=BIRD_BAND, min_prom_db=6.0, min_peaks=2)

Per-minute count of narrowband tonal peaks in band.

Uses :func:ambiscape.tonality.tonal_peaks on each row of the cached per-minute PSD. Returns median peaks/min, the per-minute counts, and the fraction of minutes with at least min_peaks (an "active" biophonic minute).

On its own this is a weak biophony discriminator and must not be read as a bird count: (a) minute-averaging smears frequency-swept birdsong, so a busy dawn chorus shows only a modest per-minute peak count while its max and active_minute_fraction spike; (b) steady machine harmonics are narrowband too and score just as high. What separates birds from machines is the combination with the temporal (:func:band_temporal_entropy) and spatial (:func:spatial_dispersion) measures — biophony is narrowband and bursty and spread across elevated bearings; a machine tone is narrowband but steady, low, and directional.

Source code in src/ambiscape/biophony.py

def narrowband_activity(F: dict, band=BIRD_BAND, min_prom_db: float = 6.0,
                        min_peaks: int = 2) -> dict:
    """Per-minute count of narrowband tonal peaks in ``band``.

    Uses :func:`ambiscape.tonality.tonal_peaks` on each row of the cached
    per-minute PSD. Returns median peaks/min, the per-minute counts, and
    the fraction of minutes with at least ``min_peaks`` (an "active"
    biophonic minute).

    On its own this is a *weak* biophony discriminator and must not be read
    as a bird count: (a) minute-averaging smears frequency-swept birdsong,
    so a busy dawn chorus shows only a modest per-minute peak count while
    its ``max`` and ``active_minute_fraction`` spike; (b) steady machine
    harmonics are narrowband too and score just as high. What separates
    birds from machines is the *combination* with the temporal
    (:func:`band_temporal_entropy`) and spatial
    (:func:`spatial_dispersion`) measures — biophony is narrowband **and**
    bursty **and** spread across elevated bearings; a machine tone is
    narrowband but steady, low, and directional.
    """
    from .tonality import tonal_peaks
    ms, freqs = np.asarray(F["minspec"], float), np.asarray(F["freqs"], float)
    counts = np.array([
        len(tonal_peaks(ms[i], freqs, fmin=band[0], fmax=band[1],
                        min_prom_db=min_prom_db)[0])
        for i in range(ms.shape[0])])
    return {
        "median_peaks_per_min": float(np.median(counts)) if len(counts) else 0.0,
        "max_peaks_per_min": int(counts.max()) if len(counts) else 0,
        "active_minute_fraction": round(float((counts >= min_peaks).mean()), 2)
        if len(counts) else 0.0,
        "per_min": counts,
    }

band_temporal_entropy(F, band=BIRD_BAND)

Sueur temporal entropy Ht of the bird-band envelope, in [0, 1].

Low = energy concentrated in time (structured vocalization); near 1 = even over time (steady band, no biophonic events).

Source code in src/ambiscape/biophony.py

def band_temporal_entropy(F: dict, band=BIRD_BAND) -> float:
    """Sueur temporal entropy Ht of the bird-band envelope, in [0, 1].

    Low = energy concentrated in time (structured vocalization); near 1 =
    even over time (steady band, no biophonic events).
    """
    env = _band_envelope(F, band)
    p = env / (env.sum() + EPS)
    return round(float(-(p * np.log(p + EPS)).sum() / np.log(len(p))), 3) \
        if len(p) > 1 else 0.0

band_activity(F, band=BIRD_BAND, k_db=3.0, bg_win_s=300.0, min_dur_s=1)

Acoustic activity of the bird band above its running background.

The band envelope (dB) is compared to a running 10th-percentile background over bg_win_s; seconds exceeding it by k_db are active. Returns the active-second fraction, event rate per minute (runs of >= min_dur_s active seconds), and median event duration.

Source code in src/ambiscape/biophony.py

def band_activity(F: dict, band=BIRD_BAND, k_db: float = 3.0,
                  bg_win_s: float = 300.0, min_dur_s: int = 1) -> dict:
    """Acoustic activity of the bird band above its running background.

    The band envelope (dB) is compared to a running 10th-percentile
    background over ``bg_win_s``; seconds exceeding it by ``k_db`` are
    active. Returns the active-second fraction, event rate per minute
    (runs of >= ``min_dur_s`` active seconds), and median event duration.
    """
    env_db = 10 * np.log10(_band_envelope(F, band) + EPS)
    n = max(3, int(round(bg_win_s)) | 1)
    bg = percentile_filter(env_db, 10, size=n, mode="nearest")
    active = env_db > bg + k_db
    lab, nlab = label(active)
    durs = [int((lab == i).sum()) for i in range(1, nlab + 1)]
    durs = [d for d in durs if d >= min_dur_s]
    dur_min = max(len(env_db) / 60.0, 1e-9)
    return {
        "active_fraction": round(float(active.mean()), 3),
        "event_rate_per_min": round(len(durs) / dur_min, 1),
        "event_median_dur_s": round(float(np.median(durs)), 1) if durs else None,
    }

spatial_dispersion(F, band=BIRD_BAND, nbins=36, limit_deg=10.0)

Directional spread and elevation of the bird-band foreground.

The azimuth histogram is weighted by the per-second bird-band foreground energy (band level above its running background), so only seconds carrying biophonic energy contribute. Returns the normalized directional entropy (0 = one bearing, 1 = all around) and the fraction of that foreground energy arriving from above limit_deg elevation (birds aloft) versus at/below the horizon.

Source code in src/ambiscape/biophony.py

def spatial_dispersion(F: dict, band=BIRD_BAND, nbins: int = 36,
                       limit_deg: float = 10.0) -> dict:
    """Directional spread and elevation of the bird-band foreground.

    The azimuth histogram is weighted by the per-second bird-band
    foreground energy (band level above its running background), so only
    seconds carrying biophonic energy contribute. Returns the normalized
    directional entropy (0 = one bearing, 1 = all around) and the fraction
    of that foreground energy arriving from above ``limit_deg`` elevation
    (birds aloft) versus at/below the horizon.
    """
    w = _fg_weights(F, band)
    az, el = np.asarray(F["az"], float), np.asarray(F["el"], float)
    tot = w.sum() + EPS
    h, _ = np.histogram(az, bins=nbins, range=(-180, 180), weights=w)
    q = h / (h.sum() + EPS)
    ent = float(-(q * np.log(q + EPS)).sum() / np.log(nbins))
    return {
        "directional_entropy": round(ent, 3),
        "above_horizon_fraction": round(float(w[el > limit_deg].sum() / tot), 2),
        "below_horizon_fraction": round(float(w[el < -limit_deg].sum() / tot), 2),
    }

summarize_biophony(F, band=BIRD_BAND, min_active_fraction=0.02)

Biophony descriptors for the analyze summary.

The spatial biophony descriptors are only meaningful when the bird band actually carries foreground energy: in a quiet, birdless room a trickle of high-frequency energy that happens to arrive from above would otherwise read as above_horizon_fraction = 1.0, a false positive. When the band-active fraction is below min_active_fraction the directional and horizon descriptors are therefore reported as None rather than as spurious numbers.

Source code in src/ambiscape/biophony.py

def summarize_biophony(F: dict, band=BIRD_BAND,
                       min_active_fraction: float = 0.02) -> dict:
    """Biophony descriptors for the analyze summary.

    The *spatial* biophony descriptors are only meaningful when the bird
    band actually carries foreground energy: in a quiet, birdless room a
    trickle of high-frequency energy that happens to arrive from above would
    otherwise read as ``above_horizon_fraction`` = 1.0, a false positive.
    When the band-active fraction is below ``min_active_fraction`` the
    directional and horizon descriptors are therefore reported as ``None``
    rather than as spurious numbers.
    """
    na = narrowband_activity(F, band)
    act = band_activity(F, band)
    active = act["active_fraction"] >= min_active_fraction
    sd = spatial_dispersion(F, band) if active else {}
    return {
        "bird_peaks_per_min": na["median_peaks_per_min"],
        "bird_active_minute_fraction": na["active_minute_fraction"],
        "bird_band_activity_pct": round(act["active_fraction"] * 100, 1),
        "bird_event_rate_per_min": act["event_rate_per_min"],
        "bird_temporal_entropy": band_temporal_entropy(F, band),
        "bird_directional_entropy": sd.get("directional_entropy"),
        "bird_above_horizon_fraction": sd.get("above_horizon_fraction"),
    }

State-resolved descriptors

State-resolved descriptors: summarize each state of a session separately.

A single descriptor row for a multi-state session is a duration-weighted average of things that never coexisted — the Haarlem loft's row is dominated by the 9-hour air-pump night and barely reflects the hi-fi afternoon. This module slices the cached features by time and runs the full summary pipeline on each state, so "vent on" and "vent off" (or day / night, or any supplied intervals) get their own complete descriptor set.

• slice_features — a sub-F restricted to a set of time intervals, valid across every feature axis (1 s, 125 ms fast, 20 ms envelope, per-minute PSD) for all the summarize_* functions;
• full_summary — the merged descriptor dict (level/event, spectral foreground, ecoacoustic, spatial, biophony) — the same set analyze writes, computed on any F;
• resolve — {state: full_summary} for a dict of named intervals;
• machine_states / diel_states — auto-discover the states from a machine band (via :func:ambiscape.states.state_segments) or the wall clock (day / night).

Event detection and percentiles run per state, so an interval need not be contiguous; a state shorter than a few frames is skipped.

intervals_from_mask(t, mask)

Contiguous [start, stop) intervals (in t units) of a boolean mask over the 1 s frames.

Source code in src/ambiscape/resolve.py

def intervals_from_mask(t: np.ndarray, mask: np.ndarray) -> list:
    """Contiguous ``[start, stop)`` intervals (in ``t`` units) of a boolean
    mask over the 1 s frames."""
    t = np.asarray(t, float)
    m = np.asarray(mask, bool)
    if not m.any():
        return []
    edges = np.flatnonzero(np.diff(m.astype(int)))
    starts = [0] if m[0] else []
    stops = []
    for e in edges:
        (stops if m[e] else starts).append(e + 1)
    if m[-1]:
        stops.append(len(m))
    return [(float(t[a]), float(t[b - 1]) + 1.0) for a, b in zip(starts, stops)]

slice_features(F, intervals)

Restrict cached features to intervals (absolute seconds).

Returns a sub-F with every time-indexed array masked to the intervals and every scalar/axis array copied through — accepted by all the summarize_* functions.

Source code in src/ambiscape/resolve.py

def slice_features(F: dict, intervals: list) -> dict:
    """Restrict cached features to ``intervals`` (absolute seconds).

    Returns a sub-``F`` with every time-indexed array masked to the
    intervals and every scalar/axis array copied through — accepted by all
    the ``summarize_*`` functions.
    """
    out = {}
    m1 = _axis_mask(F["t"], intervals)
    out["t"] = F["t"][m1]
    for k in _SEC_KEYS:
        if k in F:
            out[k] = F[k][m1]
    if "t_fast" in F:
        mf = _axis_mask(F["t_fast"], intervals)
        out["t_fast"] = F["t_fast"][mf]
        for k in _FAST_KEYS:
            if k in F:
                out[k] = F[k][mf]
    if "t_hi" in F and "env_hi" in F:
        mh = _axis_mask(F["t_hi"], intervals)
        out["t_hi"] = F["t_hi"][mh]
        out["env_hi"] = F["env_hi"][mh]
    if "min_t" in F:
        mm = _axis_mask(F["min_t"], intervals)
        out["min_t"] = F["min_t"][mm]
        out["minspec"] = F["minspec"][mm]
    for k in _SCALAR_KEYS:
        if k in F:
            out[k] = F[k]
    return out

full_summary(F)

The complete analyze descriptor set for any F (no calibration).

Source code in src/ambiscape/resolve.py

def full_summary(F: dict) -> dict:
    """The complete analyze descriptor set for any ``F`` (no calibration)."""
    from . import analysis, background, biophony, ecology, spatial
    s = analysis.summarize(F)
    s.update(background.summarize_foreground(F))
    s.update(ecology.summarize_ecology(F))
    s.update(spatial.summarize_spatial(F))
    s.update(biophony.summarize_biophony(F))
    return s

resolve(F, states, min_frames=30)

{state: full_summary} for a dict of {label: intervals}.

States whose sliced 1 s length is below min_frames are skipped (too short for stable percentiles/events).

Source code in src/ambiscape/resolve.py

def resolve(F: dict, states: dict, min_frames: int = 30) -> dict:
    """``{state: full_summary}`` for a dict of ``{label: intervals}``.

    States whose sliced 1 s length is below ``min_frames`` are skipped
    (too short for stable percentiles/events).
    """
    out = {}
    for label, intervals in states.items():
        sub = slice_features(F, intervals)
        if len(sub["t"]) >= min_frames:
            out[label] = full_summary(sub)
    return out

machine_states(F, band=(250.0, 1000.0), min_dur_s=120.0, labels=('machine_on', 'machine_off'))

Auto-discover on/off states from a machine band.

Segments the band level with :func:ambiscape.states.state_segments and returns {labels[0]: on-intervals, labels[1]: off-intervals} in absolute seconds (empty sides dropped).

Source code in src/ambiscape/resolve.py

def machine_states(F: dict, band=(250.0, 1000.0), min_dur_s: float = 120.0,
                   labels=("machine_on", "machine_off")) -> dict:
    """Auto-discover on/off states from a machine band.

    Segments the band level with :func:`ambiscape.states.state_segments`
    and returns ``{labels[0]: on-intervals, labels[1]: off-intervals}``
    in absolute seconds (empty sides dropped).
    """
    from .states import band_level, state_segments
    segs = state_segments(band_level(F, band), min_dur_s=min_dur_s)
    t = F["t"]
    n = len(t)
    states = {labels[0]: [], labels[1]: []}
    for s in segs:
        i0 = int(s["t0_s"])
        i1 = int(min(s["t0_s"] + s["dur_s"], n))
        if i1 <= i0:
            continue
        key = labels[0] if s["state"] == "on" else labels[1]
        states[key].append((float(t[i0]), float(t[i1 - 1]) + 1.0))
    return {k: v for k, v in states.items() if v}

auto_states(F, band=(250.0, 1000.0), min_dur_s=180.0, min_step_db=4.0, min_dur_frac=0.05)

Machine on/off states only if the session is genuinely two-state.

The gate for automatic use in analyze: returns the :func:machine_states split when both states last at least min_dur_s (and min_dur_frac of the session) and the band-level step between them is at least min_step_db; otherwise {} (a single-state session gets no state rows). Prevents spurious splitting of a flat, steady soundscape.

Source code in src/ambiscape/resolve.py

def auto_states(F: dict, band=(250.0, 1000.0), min_dur_s: float = 180.0,
                min_step_db: float = 4.0, min_dur_frac: float = 0.05) -> dict:
    """Machine on/off states *only if* the session is genuinely two-state.

    The gate for automatic use in ``analyze``: returns the
    :func:`machine_states` split when both states last at least
    ``min_dur_s`` (and ``min_dur_frac`` of the session) and the band-level
    step between them is at least ``min_step_db``; otherwise ``{}`` (a
    single-state session gets no state rows). Prevents spurious splitting of
    a flat, steady soundscape.
    """
    from .states import band_level
    st = machine_states(F, band=band, min_dur_s=min_dur_s)
    if len(st) < 2:
        return {}
    total = len(F["t"])
    floor = max(min_dur_s, min_dur_frac * total)
    if any(sum(b - a for a, b in v) < floor for v in st.values()):
        return {}
    lvl = band_level(F, band)
    meds = {}
    for label, intervals in st.items():
        m = _axis_mask(F["t"], intervals)
        meds[label] = float(np.median(lvl[m])) if m.any() else 0.0
    if max(meds.values()) - min(meds.values()) < min_step_db:
        return {}
    return st

diel_states(F, sess, night=(22, 6), labels=('night', 'day'))

Split a session into night / day by the wall clock.

night is (start_hour, end_hour) wrapping midnight; uses the session's day0 to turn absolute seconds into hour-of-day.

Source code in src/ambiscape/resolve.py

def diel_states(F: dict, sess, night=(22, 6),
                labels=("night", "day")) -> dict:
    """Split a session into night / day by the wall clock.

    ``night`` is ``(start_hour, end_hour)`` wrapping midnight; uses the
    session's ``day0`` to turn absolute seconds into hour-of-day.
    """
    import datetime as _dt
    base = _dt.datetime.combine(sess.day0, _dt.time())
    hours = np.array([(base + _dt.timedelta(seconds=float(x))).hour
                      for x in F["t"]])
    lo, hi = night
    if lo < hi:
        night_mask = (hours >= lo) & (hours < hi)
    else:
        night_mask = (hours >= lo) | (hours < hi)
    out = {labels[0]: intervals_from_mask(F["t"], night_mask),
           labels[1]: intervals_from_mask(F["t"], ~night_mask)}
    return {k: v for k, v in out.items() if v}

Corpus catalog

Corpus aggregation: one cross-session table from cached summaries.

Every ambiscape analyze writes <session>/analysis/summary.json. This module collects them across a corpus folder into one table — CSV for analysis, a transposed Markdown table (descriptor rows, session columns) for a consolidated report — plus simple ranking and outlier helpers. No audio, no features: it reads only the cached summaries, so a whole corpus aggregates in milliseconds.

The Markdown layout follows the Intercontinental-database CONSOLIDATED.md convention. Sessions with differing descriptor sets (older caches, optional modules) are handled by taking the union of keys and leaving blanks where a session lacks one.

collect(corpus_dir, pattern='*/analysis/summary.json', include_states=False)

Map session name → summary dict for every summary under corpus_dir.

The session name is the top-level folder (the parent of analysis/). Unreadable or malformed files are skipped. With include_states, each session's analysis/states.json (if present) contributes extra "<session>::<state>" rows right after the pooled session row — the state-resolved corpus view.

Source code in src/ambiscape/catalog.py

def collect(corpus_dir: str | Path,
            pattern: str = "*/analysis/summary.json",
            include_states: bool = False) -> dict:
    """Map session name → summary dict for every summary under ``corpus_dir``.

    The session name is the top-level folder (the parent of ``analysis/``).
    Unreadable or malformed files are skipped. With ``include_states``, each
    session's ``analysis/states.json`` (if present) contributes extra
    ``"<session>::<state>"`` rows right after the pooled session row — the
    state-resolved corpus view.
    """
    corpus_dir = Path(corpus_dir)
    out = {}
    for p in sorted(corpus_dir.glob(pattern)):
        name = p.parent.parent.name
        try:
            out[name] = json.loads(p.read_text())
        except (json.JSONDecodeError, OSError):
            continue
        if include_states:
            sp = p.parent / "states.json"
            if sp.exists():
                try:
                    states = json.loads(sp.read_text()).get("states", {})
                except (json.JSONDecodeError, OSError):
                    states = {}
                for label, summ in states.items():
                    out[f"{name}::{label}"] = {
                        k: v for k, v in summ.items() if k != "intervals_s"}
    return out

to_csv(collected, path, keys=None)

Write a session-per-row CSV (union of keys, blanks for missing).

Uses the standard :mod:csv writer, so any field containing a comma, quote, or newline is quoted correctly (a value like "interior, morning" round-trips) — consistent with every other CSV the toolkit emits.

Source code in src/ambiscape/catalog.py

def to_csv(collected: dict, path: str | Path, keys=None) -> Path:
    """Write a session-per-row CSV (union of keys, blanks for missing).

    Uses the standard :mod:`csv` writer, so any field containing a comma,
    quote, or newline is quoted correctly (a value like ``"interior,
    morning"`` round-trips) — consistent with every other CSV the toolkit
    emits.
    """
    import csv
    path = Path(path)
    cols = _all_keys(collected, keys)
    with open(path, "w", newline="") as f:
        w = csv.writer(f)
        w.writerow(["session", *cols])
        for name, summary in collected.items():
            w.writerow([name, *("" if summary.get(k) is None else summary.get(k, "")
                                for k in cols)])
    return path

to_markdown(collected, keys=None, labels=None)

Transposed Markdown table: one row per descriptor, one column per session (the consolidated-report layout). labels optionally maps descriptor keys to human labels.

Source code in src/ambiscape/catalog.py

def to_markdown(collected: dict, keys=None, labels: dict | None = None) -> str:
    """Transposed Markdown table: one row per descriptor, one column per
    session (the consolidated-report layout). ``labels`` optionally maps
    descriptor keys to human labels."""
    names = list(collected)
    cols = _all_keys(collected, keys)
    labels = labels or {}
    head = "| Descriptor | " + " | ".join(names) + " |"
    sep = "|" + "---|" * (len(names) + 1)
    rows = [head, sep]
    for k in cols:
        cells = []
        for n in names:
            v = collected[n].get(k)
            cells.append("" if v is None else str(v))
        rows.append(f"| {labels.get(k, k)} | " + " | ".join(cells) + " |")
    return "\n".join(rows)

rank(collected, key, descending=True)

(session, value) pairs sorted by key (numeric sessions only).

Source code in src/ambiscape/catalog.py

def rank(collected: dict, key: str, descending: bool = True) -> list:
    """(session, value) pairs sorted by ``key`` (numeric sessions only)."""
    vals = _values(collected, key)
    return sorted(vals.items(), key=lambda kv: kv[1], reverse=descending)

outliers(collected, key, z=1.5)

(session, z-score) for sessions more than z SDs from the mean on key, most extreme first — the cheap "what stands out" query.

Source code in src/ambiscape/catalog.py

def outliers(collected: dict, key: str, z: float = 1.5) -> list:
    """(session, z-score) for sessions more than ``z`` SDs from the mean
    on ``key``, most extreme first — the cheap "what stands out" query."""
    vals = _values(collected, key)
    if len(vals) < 3:
        return []
    x = np.array(list(vals.values()))
    mu, sd = float(x.mean()), float(x.std())
    if sd == 0:
        return []
    scored = [(n, round((v - mu) / sd, 2)) for n, v in vals.items()]
    return sorted([s for s in scored if abs(s[1]) >= z],
                  key=lambda s: -abs(s[1]))

Longitudinal analysis

Longitudinal analysis: how a place sounds across weeks, months, a year.

The unit here is the dated session summary, not the audio. A year-long study is best run as many short sessions---one a day, say, as in the StillStanding archive---each analyzed to a small summary.json; a year is then 365 tiny rows, so the longitudinal analysis is inherently out-of-core however large the underlying audio was. (A single continuous multi-month recording exceeds what the in-memory feature pipeline can hold; the supported path for year-scale work is to segment it into per-day sessions first.)

• collect_series --- read every session's summary under a corpus, ordered by date (from a date field, else parsed from the folder name), into per-descriptor time series;
• decompose --- additive split of one descriptor into a slow trend (day-windowed rolling median), a repeating seasonal component (monthly climatology of the detrended series), and the residual;
• seasonal_climatology / trend_slope --- the two components on their own (per-month means; long-term change per year);
• summarize_longitudinal --- trend per year, seasonal amplitude, peak and trough months, span;
• render --- a figure: the descriptor over time with its trend, plus the monthly climatology.

Everything is numpy-only. The motivating example is already in the StillStanding data: bird mentions peak in July and fall to zero in winter--- not because the birds leave, but because the windows close.

parse_date(name)

Parse a leading YYYY-MM-DD or YYYYMMDD date from a string.

Source code in src/ambiscape/longitudinal.py

def parse_date(name: str) -> _dt.date | None:
    """Parse a leading ``YYYY-MM-DD`` or ``YYYYMMDD`` date from a string."""
    for rx in (_ISO, _COMPACT):
        m = rx.search(name)
        if m:
            try:
                return _dt.date(int(m[1]), int(m[2]), int(m[3]))
            except ValueError:
                continue
    return None

collect_series(corpus_dir, keys=None, pattern='*/analysis/summary.json')

Dated, date-ordered descriptor time series from a corpus of sessions.

Each session's date comes from the summary's date field if present, otherwise from a date parsed out of the session folder name; sessions with no resolvable date are skipped. Returns {"dates": [date, ...], "sessions": [name, ...], "series": {key: np.array}} with all arrays in date order. keys limits the descriptors (default: the union across sessions).

Source code in src/ambiscape/longitudinal.py

def collect_series(corpus_dir: str | Path, keys=None,
                   pattern: str = "*/analysis/summary.json") -> dict:
    """Dated, date-ordered descriptor time series from a corpus of sessions.

    Each session's date comes from the summary's ``date`` field if present,
    otherwise from a date parsed out of the session folder name; sessions
    with no resolvable date are skipped. Returns ``{"dates": [date, ...],
    "sessions": [name, ...], "series": {key: np.array}}`` with all arrays in
    date order. ``keys`` limits the descriptors (default: the union across
    sessions).
    """
    import json
    corpus_dir = Path(corpus_dir)
    rows = []
    for p in sorted(corpus_dir.glob(pattern)):
        name = p.parent.parent.name
        try:
            summ = json.loads(p.read_text())
        except (json.JSONDecodeError, OSError):
            continue
        d = None
        if isinstance(summ.get("date"), str):
            d = parse_date(summ["date"])
        if d is None:
            d = parse_date(name)
        if d is not None:
            rows.append((d, name, summ))
    rows.sort(key=lambda r: r[0])
    if keys is None:
        keys, seen = [], set()
        for _d, _n, s in rows:
            for k, v in s.items():
                if k not in seen and isinstance(v, (int, float)):
                    seen.add(k)
                    keys.append(k)
    series = {}
    for k in keys:
        series[k] = np.array([float(s[k]) if isinstance(s.get(k), (int, float))
                              else np.nan for _d, _n, s in rows])
    return {"dates": [r[0] for r in rows],
            "sessions": [r[1] for r in rows], "series": series}

rolling_trend(dates, values, window_days=365.0)

Day-windowed rolling median: a slow trend robust to spikes and to irregular sampling (each point is the median of all points within window_days of it).

Source code in src/ambiscape/longitudinal.py

def rolling_trend(dates, values, window_days: float = 365.0) -> np.ndarray:
    """Day-windowed rolling median: a slow trend robust to spikes and to
    irregular sampling (each point is the median of all points within
    ``window_days`` of it)."""
    t = _ordinals(dates)
    y = np.asarray(values, float)
    out = np.empty(len(y))
    for i in range(len(y)):
        m = np.abs(t - t[i]) <= window_days / 2
        out[i] = np.median(y[m])
    return out

seasonal_climatology(dates, values)

Per-calendar-month mean and count of a series (month 1..12).

Source code in src/ambiscape/longitudinal.py

def seasonal_climatology(dates, values):
    """Per-calendar-month mean and count of a series (month 1..12)."""
    ds, y = _finite(dates, values)
    months = np.array([d.month for d in ds])
    clim, counts = {}, {}
    for m in range(1, 13):
        sel = months == m
        counts[m] = int(sel.sum())
        clim[m] = float(y[sel].mean()) if sel.any() else np.nan
    return clim, counts

decompose(dates, values, window_days=365.0)

Additive decomposition: trend + seasonal (monthly) + residual.

Returns date-ordered arrays plus the monthly climatology (mean of the detrended series per calendar month, mean-centered so the seasonal component sums to ~0).

Source code in src/ambiscape/longitudinal.py

def decompose(dates, values, window_days: float = 365.0) -> dict:
    """Additive decomposition: trend + seasonal (monthly) + residual.

    Returns date-ordered arrays plus the monthly ``climatology`` (mean of the
    detrended series per calendar month, mean-centered so the seasonal
    component sums to ~0).
    """
    ds, y = _finite(dates, values)
    order = np.argsort(_ordinals(ds))
    ds = [ds[i] for i in order]
    y = y[order]
    trend = rolling_trend(ds, y, window_days)
    detr = y - trend
    clim, _counts = seasonal_climatology(ds, detr)
    vals = [v for v in clim.values() if np.isfinite(v)]
    center = float(np.mean(vals)) if vals else 0.0
    clim = {m: (v - center if np.isfinite(v) else 0.0) for m, v in clim.items()}
    months = np.array([d.month for d in ds])
    seasonal = np.array([clim[m] for m in months])
    return {"dates": ds, "trend": trend, "seasonal": seasonal,
            "residual": y - trend - seasonal, "values": y,
            "climatology": clim}

trend_slope(dates, values)

Long-term linear change per year (least-squares slope × 365.25).

Source code in src/ambiscape/longitudinal.py

def trend_slope(dates, values) -> float:
    """Long-term linear change per year (least-squares slope × 365.25)."""
    ds, y = _finite(dates, values)
    if len(y) < 2:
        return 0.0
    t = _ordinals(ds)
    A = np.vstack([t - t[0], np.ones(len(t))]).T
    slope = np.linalg.lstsq(A, y, rcond=None)[0][0]
    return float(slope * 365.25)

summarize_longitudinal(dates, values, window_days=365.0)

Trend/seasonal descriptors of one dated series.

Source code in src/ambiscape/longitudinal.py

def summarize_longitudinal(dates, values, window_days: float = 365.0) -> dict:
    """Trend/seasonal descriptors of one dated series."""
    ds, y = _finite(dates, values)
    if len(y) == 0:
        return {"n": 0}
    dec = decompose(ds, y, window_days)
    clim = {m: v for m, v in dec["climatology"].items()}
    finite_clim = {m: v for m, v in clim.items()
                   if seasonal_climatology(ds, y)[1][m] > 0}
    t = _ordinals(dec["dates"])
    return {
        "n": int(len(y)),
        "span_days": int(t[-1] - t[0]),
        "trend_per_year": round(trend_slope(ds, y), 3),
        "seasonal_amplitude": round(
            max(finite_clim.values()) - min(finite_clim.values()), 3)
        if finite_clim else 0.0,
        "peak_month": max(finite_clim, key=finite_clim.get)
        if finite_clim else None,
        "trough_month": min(finite_clim, key=finite_clim.get)
        if finite_clim else None,
        "residual_sd": round(float(np.std(dec["residual"])), 3),
    }

render(dates, values, out_path, key='', window_days=365.0)

Two-panel figure: series + rolling trend over time, and the monthly climatology.

Source code in src/ambiscape/longitudinal.py

def render(dates, values, out_path, key: str = "", window_days: float = 365.0):
    """Two-panel figure: series + rolling trend over time, and the monthly
    climatology."""
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from .figures import RC, BLUE, YELLOW, MUT

    dec = decompose(dates, values, window_days)
    ds, y, trend = dec["dates"], dec["values"], dec["trend"]
    t = _ordinals(ds)
    clim = dec["climatology"]
    with plt.rc_context(RC):
        fig, ax = plt.subplots(1, 2, figsize=(12, 4), dpi=130,
                               gridspec_kw=dict(width_ratios=[2.6, 1]))
        ax[0].plot(t, y, ".", ms=3, color=BLUE, alpha=0.5, label=key or "value")
        ax[0].plot(t, trend, color=YELLOW, lw=2,
                   label=f"trend ({window_days:.0f}-day median)")
        yr = np.arange(np.ceil(ds[0].year), ds[-1].year + 1)
        ax[0].set_xticks([_dt.date(int(y_), 1, 1).toordinal() for y_ in yr],
                         [str(int(y_)) for y_ in yr])
        ax[0].set(ylabel=key or "descriptor",
                  title=f"{key} over time")
        ax[0].legend(fontsize=8)
        months = list(range(1, 13))
        ax[1].bar(months, [clim[m] for m in months], color=BLUE, alpha=0.7)
        ax[1].axhline(0, color=MUT, lw=0.7)
        ax[1].set_xticks(months, ["J", "F", "M", "A", "M", "J", "J", "A",
                                  "S", "O", "N", "D"], fontsize=8)
        ax[1].set(title="seasonal (monthly, detrended)")
        fig.tight_layout()
        fig.savefig(out_path, bbox_inches="tight")
        plt.close(fig)
    return out_path

run_corpus(corpus_dir, out_dir, keys=None, window_days=365.0)

CLI driver: per-descriptor longitudinal summaries + a figure for each, writing longitudinal.json.

Source code in src/ambiscape/longitudinal.py

def run_corpus(corpus_dir, out_dir, keys=None, window_days: float = 365.0) -> dict:
    """CLI driver: per-descriptor longitudinal summaries + a figure for each,
    writing ``longitudinal.json``."""
    import json
    out_dir = Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    s = collect_series(corpus_dir, keys=keys)
    doc = {"n_sessions": len(s["dates"]),
           "date_range": [s["dates"][0].isoformat(), s["dates"][-1].isoformat()]
           if s["dates"] else None,
           "descriptors": {}}
    for k, y in s["series"].items():
        if np.isfinite(y).sum() >= 3:
            doc["descriptors"][k] = summarize_longitudinal(s["dates"], y,
                                                           window_days)
    (out_dir / "longitudinal.json").write_text(json.dumps(doc, indent=2))
    return doc

Visual features

Per-frame visual features: the light/vision analogue of the audio deposit.

The AMBIENT project treats a room as an audio-visual subject, and a room's look has a diurnal rhythm just as its sound does. This module extracts a compact descriptor from a single video frame so a camera can log visual behaviour rather than store imagery --- the same "features, not recordings" privacy stance as :mod:ambiscape.capture. It is numpy-only (no camera or OpenCV dependency); frame grabbing lives in the capture rig, the feature definitions live here so they are versioned and tested with ambiscape.

Per frame (:func:frame_features):

• brightness / brightness_sd --- mean and spread of Rec.709 luma (0 = dark, 1 = bright): the room's overall light level and its unevenness;
• r_frac / g_frac / b_frac, warm_cool_ratio --- colour balance and a warm/cool proxy (daylight vs incandescent vs the blue of a screen);
• saturation, colourfulness --- how colourful the scene is (Hasler--Susstrunk colourfulness), near zero for a grey room;
• spatial_entropy --- entropy of a coarse brightness grid: 1 when the room is evenly lit, low when light is concentrated (a lamp, a window) --- the visual analogue of acoustic diffuseness;
• bright_centroid_x / _y --- the luma-weighted centre of light in the frame (0..1): where the light comes from, a visual direction-of-arrival.

:func:frame_delta is a motion proxy (mean absolute luma change between two frames); :func:summarize_vision rolls a day of per-frame features into a vis_-prefixed summary that joins the audio summary.json.

luma(rgb)

Rec.709 luma of an RGB frame, in [0, 1].

Source code in src/ambiscape/vision.py

def luma(rgb) -> np.ndarray:
    """Rec.709 luma of an RGB frame, in [0, 1]."""
    x = _to_unit(rgb)
    return x @ np.array(_LUMA)

frame_features(rgb, grid=3)

Visual descriptor of one RGB frame (uint8 0..255 or float 0..1).

Source code in src/ambiscape/vision.py

def frame_features(rgb, grid: int = 3) -> dict:
    """Visual descriptor of one RGB frame (uint8 0..255 or float 0..1)."""
    x = _to_unit(rgb)
    R, G, B = x[..., 0], x[..., 1], x[..., 2]
    lum = x @ np.array(_LUMA)
    mr, mg, mb = float(R.mean()), float(G.mean()), float(B.mean())
    s = mr + mg + mb + EPS

    mx = x.max(-1)
    mn = x.min(-1)
    sat = float(np.where(mx > 0, (mx - mn) / (mx + EPS), 0.0).mean())
    rg = R - G
    yb = 0.5 * (R + G) - B
    colourfulness = float(np.sqrt(rg.std() ** 2 + yb.std() ** 2)
                          + 0.3 * np.sqrt(rg.mean() ** 2 + yb.mean() ** 2))

    gm = _grid_means(lum, grid)
    p = gm.ravel() / (gm.sum() + EPS)
    spatial_entropy = float(-(p * np.log(p + EPS)).sum() / np.log(grid * grid))

    h, w = lum.shape
    ys, xs = np.mgrid[0:h, 0:w]
    tot = lum.sum() + EPS
    cx = float((lum * xs).sum() / tot / max(w - 1, 1))
    cy = float((lum * ys).sum() / tot / max(h - 1, 1))

    return {
        "brightness": round(float(lum.mean()), 4),
        "brightness_sd": round(float(lum.std()), 4),
        "r_frac": round(mr / s, 4), "g_frac": round(mg / s, 4),
        "b_frac": round(mb / s, 4),
        "warm_cool_ratio": round(mr / (mb + EPS), 4),
        "saturation": round(sat, 4),
        "colourfulness": round(colourfulness, 4),
        "spatial_entropy": round(spatial_entropy, 4),
        "bright_centroid_x": round(cx, 4),
        "bright_centroid_y": round(cy, 4),
    }

frame_delta(prev_rgb, cur_rgb)

Mean absolute luma change between two frames (motion proxy, 0..1).

Source code in src/ambiscape/vision.py

def frame_delta(prev_rgb, cur_rgb) -> float:
    """Mean absolute luma change between two frames (motion proxy, 0..1)."""
    return round(float(np.abs(luma(cur_rgb) - luma(prev_rgb)).mean()), 5)

summarize_vision(features, motion=None)

Roll a day of per-frame feature dicts into a vis_ day summary.

Source code in src/ambiscape/vision.py

def summarize_vision(features, motion=None) -> dict:
    """Roll a day of per-frame feature dicts into a ``vis_`` day summary."""
    rows = list(features)
    if not rows:
        return {"n_frames": 0}
    def col(k):
        return np.array([r[k] for r in rows if k in r], float)
    b = col("brightness")
    out = {
        "n_frames": len(rows),
        "vis_brightness_median": round(float(np.median(b)), 4),
        "vis_brightness_range": round(float(np.percentile(b, 90)
                                            - np.percentile(b, 10)), 4),
        "vis_colourfulness_median": round(float(np.median(col("colourfulness"))), 4),
        "vis_saturation_median": round(float(np.median(col("saturation"))), 4),
        "vis_warm_cool_median": round(float(np.median(col("warm_cool_ratio"))), 4),
        "vis_spatial_entropy_median": round(
            float(np.median(col("spatial_entropy"))), 4),
    }
    if motion is not None and len(motion):
        m = np.asarray(motion, float)
        out["vis_motion_mean"] = round(float(m.mean()), 5)
        out["vis_motion_p90"] = round(float(np.percentile(m, 90)), 5)
    return out

Circular statistics

Circular statistics, shared by spatial (azimuth) and rhythm (phase) code.

Angles are radians throughout; the degree-facing wrapper lives in :func:ambiscape.analysis.circular_stats, and period-facing helpers here convert times to phases. The resultant length R in [0, 1] measures concentration; circular SD = sqrt(-2 ln R); the Rayleigh test gives the probability of R under uniformity (p ~ exp(-n R^2), adequate for n >= 10).

mean_resultant(angles, weights=None)

Weighted circular mean (rad) and resultant length R.

Source code in src/ambiscape/circstats.py

def mean_resultant(angles: np.ndarray, weights=None):
    """Weighted circular mean (rad) and resultant length R."""
    a = np.asarray(angles, float)
    w = np.ones_like(a) if weights is None else np.asarray(weights, float)
    z = (w * np.exp(1j * a)).sum() / (w.sum() + EPS)
    return float(np.angle(z)), float(np.abs(z))

circular_sd(R)

Circular standard deviation (rad) from a resultant length.

Source code in src/ambiscape/circstats.py

def circular_sd(R: float) -> float:
    """Circular standard deviation (rad) from a resultant length."""
    return float(np.sqrt(-2 * np.log(max(R, EPS))))

rayleigh_p(R, n)

Rayleigh-test p-value (uniformity null), small-p approximation.

Source code in src/ambiscape/circstats.py

def rayleigh_p(R: float, n: int) -> float:
    """Rayleigh-test p-value (uniformity null), small-p approximation."""
    z = n * R * R
    return float(np.clip(np.exp(-z) * (1 + (2 * z - z * z) / (4 * n)),
                         0.0, 1.0))

phase_stats(times, period)

Circular statistics of event times folded at period.

Returns mean phase (cycles), R, circular SD in seconds, and the Rayleigh p-value — the standard summary for one strike stream.

Source code in src/ambiscape/circstats.py

def phase_stats(times: np.ndarray, period: float) -> dict:
    """Circular statistics of event times folded at ``period``.

    Returns mean phase (cycles), R, circular SD in seconds, and the
    Rayleigh p-value — the standard summary for one strike stream.
    """
    ph = 2 * np.pi * (np.asarray(times, float) / period % 1.0)
    mu, R = mean_resultant(ph)
    return {
        "mean_phase": float((mu / (2 * np.pi)) % 1.0),
        "R": round(R, 4),
        "circ_sd_s": round(circular_sd(R) / (2 * np.pi) * period, 4),
        "rayleigh_p": rayleigh_p(R, len(ph)),
        "n": int(len(ph)),
    }

relative_phase(times, ref_times, period)

Phase of each event relative to the preceding reference event.

The per-event lock between two streams sharing one period: mean offset (cycles and seconds), R, and circular SD in seconds. R near 1 means the two streams are phase-locked at strike level.

Source code in src/ambiscape/circstats.py

def relative_phase(times: np.ndarray, ref_times: np.ndarray,
                   period: float) -> dict:
    """Phase of each event relative to the preceding reference event.

    The per-event lock between two streams sharing one period: mean offset
    (cycles and seconds), R, and circular SD in seconds. R near 1 means the
    two streams are phase-locked at strike level.
    """
    ref = np.sort(np.asarray(ref_times, float))
    t = np.asarray(times, float)
    i = np.clip(np.searchsorted(ref, t) - 1, 0, len(ref) - 1)
    d = 2 * np.pi * (((t - ref[i]) / period) % 1.0)
    mu, R = mean_resultant(d)
    off = (mu / (2 * np.pi)) % 1.0
    return {
        "mean_offset_cycles": round(off, 4),
        "mean_offset_s": round(off * period, 4),
        "R": round(R, 4),
        "circ_sd_s": round(circular_sd(R) / (2 * np.pi) * period, 4),
        "n": int(len(t)),
    }

Modulation profile

Environmental rhythm: multi-scale envelope modulation profile.

Soundscapes are rhythmic on very different time scales at once — strike patterns (micro), traffic waves and surf (meso), duty cycles of machines and human activity (macro). This module measures all three from cached envelopes, no audio pass:

• micro (0.5–20 Hz) from the 20 ms broadband envelope (env_hi, extractor ≥ 0.2 caches; older caches fall back to the 8 Hz fast level, which limits micro to < 4 Hz);
• meso (0.01–0.5 Hz) from the 125 ms fast level;
• macro (below 0.01 Hz, floor set by session length) from the 1 s RMS.

profile returns, per scale, a log-frequency modulation spectrum with the dominant modulation frequency, its prominence, and the band modulation depth. modulation_spectrogram computes the windowed version — the "rhythm spectrogram of the day" — and render writes the combined figure.

modulation_spectrum(env, dt, fmin, fmax, n_bins=48)

Welch modulation spectrum of a linear-power envelope, log-resampled.

The envelope is normalized to zero-mean unit-mean (x = env/mean − 1) so spectra are comparable across levels; returns (freqs, power density).

Source code in src/ambiscape/modulation.py

def modulation_spectrum(env: np.ndarray, dt: float, fmin: float, fmax: float,
                        n_bins: int = 48):
    """Welch modulation spectrum of a linear-power envelope, log-resampled.

    The envelope is normalized to zero-mean unit-mean (x = env/mean − 1) so
    spectra are comparable across levels; returns (freqs, power density).
    """
    x = env.astype(np.float64) / (env.mean() + EPS) - 1.0
    nper = int(min(len(x), max(64, round(8.0 / (fmin * dt)))))
    f, P = signal.welch(x, fs=1.0 / dt, nperseg=nper,
                        noverlap=nper // 2, detrend="linear")
    grid = np.geomspace(fmin, fmax, n_bins)
    idx = np.clip(np.searchsorted(f, grid), 1, len(f) - 1)
    return grid, P[idx]

profile(F)

Three-scale modulation profile from cached features.

Source code in src/ambiscape/modulation.py

def profile(F: dict) -> dict:
    """Three-scale modulation profile from cached features."""
    dur = float(len(F["t"]))
    out = {"scales": {}, "spectra": {}}
    sources = {
        "meso": (F["fast_db"], float(np.median(np.diff(F["t_fast"])))),
        "macro": (F["rms_w"] ** 2, 1.0),
    }
    if "env_hi" in F:
        sources["micro"] = (F["env_hi"], float(F["hi_dt"]))
    else:
        sources["micro"] = (F["fast_db"], sources["meso"][1])
        out["micro_limited"] = "no env_hi in cache; micro band tops out at 4 Hz"
    for scale in SCALES:
        env, dt = sources[scale]
        if scale in ("meso", "macro"):     # dB level -> linear power
            env = 10 ** (env.astype(np.float64) / 10) if scale == "meso" else env
        lo, hi = BANDS[scale]
        lo = max(lo or 4.0 / dur, 4.0 / dur)
        hi = min(hi, 0.45 / dt)
        if hi <= lo * 1.5:
            continue
        f, P = modulation_spectrum(env, dt, lo, hi)
        out["spectra"][scale] = {"freq_hz": [round(float(v), 5) for v in f],
                                 "power": [float(v) for v in P]}
        out["scales"][scale] = _scale_stats(f, P)
    return out

modulation_spectrogram(env, dt, win_s=600.0, step_s=120.0, fmin=0.02, fmax=20.0, n_bins=64)

Windowed modulation spectra: the rhythm spectrogram of the session.

Returns (t_centers, mod_freqs, S) with S in dB relative to each window's median (so rhythmic structure reads as ridges regardless of level).

Source code in src/ambiscape/modulation.py

def modulation_spectrogram(env: np.ndarray, dt: float, win_s: float = 600.0,
                           step_s: float = 120.0, fmin: float = 0.02,
                           fmax: float = 20.0, n_bins: int = 64):
    """Windowed modulation spectra: the rhythm spectrogram of the session.

    Returns (t_centers, mod_freqs, S) with S in dB relative to each window's
    median (so rhythmic structure reads as ridges regardless of level).
    """
    nwin = int(win_s / dt)
    nstep = int(step_s / dt)
    fmax = min(fmax, 0.45 / dt)
    grid = np.geomspace(fmin, fmax, n_bins)
    ts, rows = [], []
    for i0 in range(0, len(env) - nwin + 1, nstep):
        f, P = modulation_spectrum(env[i0:i0 + nwin], dt, fmin, fmax, n_bins)
        rows.append(10 * np.log10((P + EPS) / (np.median(P) + EPS)))
        ts.append((i0 + nwin / 2) * dt)
    return np.array(ts), grid, np.array(rows)

render(F, prof, out_path, title='', clock=None)

Combined figure: per-scale spectra + rhythm spectrogram.

Source code in src/ambiscape/modulation.py

def render(F: dict, prof: dict, out_path, title="", clock=None):
    """Combined figure: per-scale spectra + rhythm spectrogram."""
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    fig, (ax0, ax1) = plt.subplots(
        2, 1, figsize=(12.8, 7.2), dpi=130,
        gridspec_kw=dict(height_ratios=[1, 1.3]))
    colors = {"micro": "#2a78d6", "meso": "#d66a2a", "macro": "#3d9970"}
    for scale in SCALES:
        sp = prof["spectra"].get(scale)
        if not sp:
            continue
        f = np.array(sp["freq_hz"])
        P = np.array(sp["power"])
        ax0.plot(f, 10 * np.log10(P + EPS), color=colors[scale], lw=1.4,
                 label=f"{scale} (peak {prof['scales'][scale]['peak_period_s']} s)")
    ax0.set(xscale="log", xlabel="modulation frequency (Hz)",
            ylabel="power (dB)", title=f"{title} — envelope modulation "
            "spectra by scale")
    ax0.legend(fontsize=8)
    ax0.grid(alpha=0.2, which="both")

    env, dt = (F["env_hi"], float(F["hi_dt"])) if "env_hi" in F else \
        (10 ** (F["fast_db"].astype(np.float64) / 10),
         float(np.median(np.diff(F["t_fast"]))))
    t0 = float(F["t_hi"][0] if "t_hi" in F else F["t_fast"][0])
    ts, mf, S = modulation_spectrogram(env, dt)
    if len(ts):
        pc = ax1.pcolormesh(t0 + ts, mf, S.T, cmap="magma", shading="auto",
                            vmin=0, vmax=max(6, np.percentile(S, 99)))
        ax1.set(yscale="log", ylabel="modulation frequency (Hz)",
                title="rhythm spectrogram (10 min windows, dB re window median)")
        if clock is not None:
            xt = ax1.get_xticks()
            ax1.set_xticks(xt)
            ax1.set_xticklabels([clock(x)[7:] for x in xt], fontsize=8)
            ax1.set_xlim(t0 + ts[0], t0 + ts[-1])
        fig.colorbar(pc, ax=ax1, pad=0.01)
    fig.tight_layout()
    fig.savefig(out_path)
    return out_path

run_session(sess, out_dir)

CLI driver: profile + figure + modulation.json.

Source code in src/ambiscape/modulation.py

def run_session(sess, out_dir) -> dict:
    """CLI driver: profile + figure + modulation.json."""
    import json
    from .features import load_features
    out_dir = Path(out_dir)
    F = load_features(sorted((out_dir / "features").glob("*.npz")))
    prof = profile(F)
    render(F, prof, out_dir / "modulation_profile.png", title=sess.name,
           clock=sess.clock)
    (out_dir / "modulation.json").write_text(json.dumps(prof, indent=2))
    return prof

Tonality and harmonicity

Tonalness timeline, harmonic sieve, and inharmonicity.

Works entirely from the cached per-minute mean PSD (minspec):

• tonal_peaks — prominent narrowband components per minute (dB above a running spectral floor), the raw material for everything below;
• tonal_tracks — peaks linked across minutes into tracks (a hum, a bell partial, a beep), each with duration, median frequency, and cents drift: the tonalness timeline;
• harmonic_sieve — best f0 explaining a minute's peak set as a harmonic series; the unexplained remainder is the inharmonic tonal content. Voices, engines, and music score high harmonicity; bells score low (their partial series 1 : 2 : 2.4 : 3 : 4 is not harmonic);
• pitch_class_profile — tonal peak energy folded onto the 12 pitch classes: what "key" the soundscape hums in.

run_session writes tonality.json + tonality.png.

tonal_peaks(spec_row, freqs, fmin=40.0, fmax=8000.0, min_prom_db=8.0, max_n=40)

Prominent narrowband peaks in one mean spectrum.

The floor is a wide median filter on the log spectrum; peaks must rise min_prom_db above it. Returns (freq, prominence_db, power) arrays.

Source code in src/ambiscape/tonality.py

def tonal_peaks(spec_row: np.ndarray, freqs: np.ndarray, fmin=40.0,
                fmax=8000.0, min_prom_db=8.0, max_n=40):
    """Prominent narrowband peaks in one mean spectrum.

    The floor is a wide median filter on the log spectrum; peaks must rise
    ``min_prom_db`` above it. Returns (freq, prominence_db, power) arrays.
    """
    m = (freqs >= fmin) & (freqs <= fmax)
    ls = 10 * np.log10(spec_row[m] + EPS)
    floor = median_filter(ls, size=101, mode="nearest")
    rise = ls - floor
    pk, props = find_peaks(rise, height=min_prom_db, distance=3)
    order = np.argsort(props["peak_heights"])[::-1][:max_n]
    keep = np.sort(pk[order])
    return freqs[m][keep], rise[keep], spec_row[m][keep]

tonal_tracks(minspec, freqs, tol_cents=40.0, min_minutes=2, **peak_kw)

Link per-minute peaks into tracks. Returns a list of dicts sorted by duration (longest first): median freq, span of minutes, mean prominence, and total drift in cents.

Source code in src/ambiscape/tonality.py

def tonal_tracks(minspec: np.ndarray, freqs: np.ndarray, tol_cents=40.0,
                 min_minutes=2, **peak_kw):
    """Link per-minute peaks into tracks. Returns a list of dicts sorted by
    duration (longest first): median freq, span of minutes, mean prominence,
    and total drift in cents."""
    per_min = [tonal_peaks(minspec[i], freqs, **peak_kw)
               for i in range(minspec.shape[0])]
    open_tracks, done = [], []
    for mi, (fq, prom, _pw) in enumerate(per_min):
        used = np.zeros(len(fq), bool)
        for tr in list(open_tracks):
            cents = 1200 * np.abs(np.log2(fq / (tr["f"][-1] + EPS) + EPS))
            j = int(np.argmin(cents)) if len(cents) else -1
            if j >= 0 and cents[j] < tol_cents and not used[j]:
                tr["f"].append(float(fq[j]))
                tr["prom"].append(float(prom[j]))
                tr["m1"] = mi
                used[j] = True
            elif mi - tr["m1"] > 1:
                open_tracks.remove(tr)
                done.append(tr)
        for j in np.flatnonzero(~used):
            open_tracks.append(dict(f=[float(fq[j])], prom=[float(prom[j])],
                                    m0=mi, m1=mi))
    done += open_tracks
    out = []
    for tr in done:
        if tr["m1"] - tr["m0"] + 1 < min_minutes:
            continue
        f = np.array(tr["f"])
        out.append({
            "f_median_hz": round(float(np.median(f)), 1),
            "t0_min": tr["m0"], "t1_min": tr["m1"],
            "minutes": tr["m1"] - tr["m0"] + 1,
            "prominence_db": round(float(np.mean(tr["prom"])), 1),
            "drift_cents": round(float(1200 * np.log2(
                (f[-1] + EPS) / (f[0] + EPS))), 1),
        })
    return sorted(out, key=lambda t: -t["minutes"])

harmonic_sieve(fq, power, f0_min=60.0, f0_max=1200.0, tol_cents=35.0, max_harm=12)

Best f0 explaining the peak set as harmonics k*f0.

Candidate f0s are every peak divided by k = 1..6; the score is the power-weighted fraction of peaks within tol_cents of a harmonic. Returns (f0, harmonicity in [0,1]) — harmonicity is the explained power fraction; 1 − harmonicity is the inharmonicity index.

Source code in src/ambiscape/tonality.py

def harmonic_sieve(fq: np.ndarray, power: np.ndarray, f0_min=60.0,
                   f0_max=1200.0, tol_cents=35.0, max_harm=12):
    """Best f0 explaining the peak set as harmonics k*f0.

    Candidate f0s are every peak divided by k = 1..6; the score is the
    power-weighted fraction of peaks within ``tol_cents`` of a harmonic.
    Returns (f0, harmonicity in [0,1]) — harmonicity is the explained power
    fraction; 1 − harmonicity is the inharmonicity index.
    """
    if len(fq) == 0:
        return None, 0.0
    cands = np.concatenate([fq / k for k in range(1, 7)])
    cands = cands[(cands >= f0_min) & (cands <= f0_max)]
    best_f0, best = None, 0.0
    total = power.sum() + EPS
    for f0 in cands:
        k = np.clip(np.round(fq / f0), 1, max_harm)
        cents = 1200 * np.abs(np.log2(fq / (k * f0)))
        score = float(power[cents < tol_cents].sum() / total)
        if score > best:
            best_f0, best = float(f0), score
    return best_f0, round(best, 3)

pitch_class_profile(minspec, freqs, minutes=None, **peak_kw)

Tonal peak power folded onto 12 pitch classes (A4 = 440 Hz).

Source code in src/ambiscape/tonality.py

def pitch_class_profile(minspec: np.ndarray, freqs: np.ndarray,
                        minutes=None, **peak_kw):
    """Tonal peak power folded onto 12 pitch classes (A4 = 440 Hz)."""
    pcp = np.zeros(12)
    rows = range(minspec.shape[0]) if minutes is None else minutes
    for i in rows:
        fq, _prom, pw = tonal_peaks(minspec[i], freqs, **peak_kw)
        if len(fq):
            pc = np.mod(np.round(12 * np.log2(fq / 440.0) + 69), 12).astype(int)
            np.add.at(pcp, pc, pw)
    return pcp / (pcp.sum() + EPS)

run_session(sess, out_dir)

Tonality timeline + per-minute tonalness/harmonicity + PCP + figure.

Source code in src/ambiscape/tonality.py

def run_session(sess, out_dir) -> dict:
    """Tonality timeline + per-minute tonalness/harmonicity + PCP + figure."""
    import json
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from .features import load_features

    out_dir = Path(out_dir)
    F = load_features(sorted((out_dir / "features").glob("*.npz")))
    minspec, freqs = F["minspec"], F["freqs"]
    nmin = minspec.shape[0]

    tracks = tonal_tracks(minspec, freqs)
    tonalness, harmonicity = np.zeros(nmin), np.full(nmin, np.nan)
    for i in range(nmin):
        fq, _prom, pw = tonal_peaks(minspec[i], freqs)
        band = (freqs >= 40) & (freqs <= 8000)
        tonalness[i] = float(pw.sum() / (minspec[i][band].sum() + EPS))
        if len(fq) >= 3:
            _f0, h = harmonic_sieve(fq, pw)
            harmonicity[i] = h
    pcp = pitch_class_profile(minspec, freqs)

    doc = {
        "tracks": tracks[:40],
        "tonalness_median": round(float(np.median(tonalness)), 3),
        "harmonicity_median": round(float(np.nanmedian(harmonicity)), 3),
        "inharmonicity_median": round(1 - float(np.nanmedian(harmonicity)), 3),
        "pitch_class_profile": {NOTE[i]: round(float(pcp[i]), 3)
                                for i in range(12)},
        "top_pitch_classes": [NOTE[i] for i in np.argsort(pcp)[::-1][:3]],
    }
    (out_dir / "tonality.json").write_text(json.dumps(doc, indent=2))

    fig, ax = plt.subplots(1, 2, figsize=(12.8, 4.6), dpi=130,
                           gridspec_kw=dict(width_ratios=[2.4, 1]))
    for tr in tracks:
        ax[0].plot([tr["t0_min"], tr["t1_min"] + 1],
                   [tr["f_median_hz"]] * 2, lw=max(0.8, tr["prominence_db"] / 6),
                   color="#2a78d6", alpha=0.7, solid_capstyle="butt")
    ax[0].set(yscale="log", xlabel="minute of session", ylabel="Hz",
              title=f"{sess.name} — tonal tracks (width = prominence)")
    ax[0].grid(alpha=0.2, which="both")
    ax[1].bar(range(12), pcp, color="#2a78d6")
    ax[1].set_xticks(range(12), NOTE, fontsize=8)
    ax[1].set(title="pitch-class profile", ylabel="tonal power share")
    ax[1].grid(alpha=0.2, axis="y")
    fig.tight_layout()
    fig.savefig(out_dir / "tonality.png")
    plt.close(fig)
    return doc

Strike-level rhythm

Strike-level rhythm analysis of quasi-periodic pitched sources (bells, machines, signals) in long ambisonic recordings.

The 1 Hz features in :mod:features are too coarse for strike rhythm, so this module makes one extra streaming pass at ~20 ms resolution, restricted to the narrowband partials of the sources of interest:

1. detect_partials — narrowband peaks from the cached per-minute mean PSD, contrasting source-active against quiet minutes;
2. partial_pass — streaming STFT pass storing, per frame, the power envelope and pseudo-intensity (for DOA) at each partial, plus a broadband spectral-flux onset function;
3. cluster_partials — group partials into sources by correlating their half-wave-rectified log-envelope derivatives (strike-synchronous);
4. pick_strikes — adaptive, strongest-first onset picking per source;
5. rayleigh_period / period_track — point-process periodicity (resultant length over a period grid, with harmonics for multi-strike cycles);
6. cycle_grid — repetition-with-variation statistics: per-cycle timing residuals and amplitudes for each position of the repeating pattern;
7. rise_spectrum — strike-triggered post/pre spectral rise, the tonal content of one rhythmic position (also exposes cross-talk: a stream whose rise spectrum shows another source's partials is leakage, not a strike).

All statistics are on the W channel except pseudo-intensity, which follows the AmbiX ACN (W, Y, Z, X) convention used throughout ambiscape.

detect_partials(F, active, quiet, band=(350.0, 4500.0), min_rise_db=6.0, max_n=30)

Narrowband partials of the active-state source(s).

active/quiet are boolean masks over the minutes of F["minspec"] (source clearly present / clearly absent). Returns (freqs, rise_db) sorted by frequency.

Source code in src/ambiscape/rhythm.py

def detect_partials(F: dict, active: np.ndarray, quiet: np.ndarray,
                    band=(350.0, 4500.0), min_rise_db=6.0, max_n=30):
    """Narrowband partials of the active-state source(s).

    ``active``/``quiet`` are boolean masks over the minutes of ``F["minspec"]``
    (source clearly present / clearly absent). Returns (freqs, rise_db)
    sorted by frequency.
    """
    S_a = F["minspec"][active].mean(0)
    S_q = F["minspec"][quiet].mean(0)
    freqs = F["freqs"]
    m = (freqs >= band[0]) & (freqs <= band[1])
    ratio = 10 * np.log10((S_a[m] + EPS) / (S_q[m] + EPS))
    pk, props = find_peaks(ratio, height=min_rise_db, prominence=5.0,
                           distance=5)
    order = np.argsort(props["peak_heights"])[::-1][:max_n]
    keep = np.sort(pk[order])
    return freqs[m][keep], ratio[keep]

partial_pass(take, pfreq, nfft=4096, hop=960)

Streaming STFT pass; per ~20 ms frame: 3-bin power envelope and pseudo-intensity at each partial, plus a 400-4500 Hz log-flux onset function. Returns dict of arrays keyed t/env/ix/iy/iz/odf/pfreq.

Source code in src/ambiscape/rhythm.py

def partial_pass(take: Take, pfreq, nfft=4096, hop=960) -> dict:
    """Streaming STFT pass; per ~20 ms frame: 3-bin power envelope and
    pseudo-intensity at each partial, plus a 400-4500 Hz log-flux onset
    function. Returns dict of arrays keyed t/env/ix/iy/iz/odf/pfreq."""
    fs = take.samplerate
    hop = int(hop * fs / 48000)
    win = np.hanning(nfft).astype(np.float32)
    freqs = np.fft.rfftfreq(nfft, 1 / fs)
    pbin = np.array([int(round(f * nfft / fs)) for f in np.atleast_1d(pfreq)])
    psel = np.stack([pbin - 1, pbin, pbin + 1], 1)
    fmask = (freqs >= 400) & (freqs <= 4500)

    env, ix, iy, iz, odf = [], [], [], [], []
    prevL = None
    carry = np.zeros((0, take.channels), np.float32)
    iW, iY, iZ, iX = take.wyzx
    with sf.SoundFile(str(take.path)) as f:
        while True:
            block = f.read(60 * fs, dtype="float32", always_2d=True)
            if block.shape[0] == 0:
                break
            data = np.concatenate([carry, block]) if carry.shape[0] else block
            nwin = (data.shape[0] - nfft) // hop + 1
            if nwin <= 0:
                carry = data
                continue
            idx = np.arange(nfft)[None, :] + hop * np.arange(nwin)[:, None]
            W = np.fft.rfft(data[:, iW][idx] * win)
            Y = np.fft.rfft(data[:, iY][idx] * win)
            Z = np.fft.rfft(data[:, iZ][idx] * win)
            X = np.fft.rfft(data[:, iX][idx] * win)
            Pw = W.real ** 2 + W.imag ** 2
            env.append(Pw[:, psel].sum(2).astype(np.float32))
            ix.append((W.conj() * X).real[:, psel].sum(2).astype(np.float32))
            iy.append((W.conj() * Y).real[:, psel].sum(2).astype(np.float32))
            iz.append((W.conj() * Z).real[:, psel].sum(2).astype(np.float32))
            L = np.log1p(1e6 * Pw[:, fmask])
            Lp = np.concatenate([prevL[None] if prevL is not None else L[:1], L])
            odf.append(np.maximum(np.diff(Lp, axis=0), 0).sum(1)
                       .astype(np.float32))
            prevL = L[-1]
            carry = data[nwin * hop:].copy()
    odf = np.concatenate(odf)
    return dict(t=(np.arange(len(odf)) * hop + nfft // 2) / fs,
                env=np.concatenate(env), ix=np.concatenate(ix),
                iy=np.concatenate(iy), iz=np.concatenate(iz),
                odf=odf, pfreq=np.atleast_1d(pfreq).astype(np.float64))

cluster_partials(env, mask=None, th=0.75, min_size=3)

Group partial columns into sources: correlate rectified log-envelope derivatives (with ±1-frame jitter tolerance), average-linkage cluster. Returns a list of column-index lists, largest first; singletons and groups below min_size are dropped (assign them afterwards with strike-triggered statistics if needed).

Source code in src/ambiscape/rhythm.py

def cluster_partials(env: np.ndarray, mask=None, th=0.75, min_size=3):
    """Group partial columns into sources: correlate rectified log-envelope
    derivatives (with ±1-frame jitter tolerance), average-linkage cluster.
    Returns a list of column-index lists, largest first; singletons and
    groups below ``min_size`` are dropped (assign them afterwards with
    strike-triggered statistics if needed)."""
    from scipy.cluster.hierarchy import fcluster, linkage
    e = env if mask is None else env[mask]
    dL = np.maximum(np.diff(np.log10(e + EPS), axis=0), 0)
    dL = maximum_filter1d(dL, 3, axis=0)
    C = np.corrcoef(dL.T)
    d = 1 - C
    np.fill_diagonal(d, 0)
    lab = fcluster(linkage(d[np.triu_indices_from(d, 1)], method="average"),
                   th, criterion="distance")
    groups = {}
    for i, l in enumerate(lab):
        groups.setdefault(l, []).append(i)
    out = [v for v in groups.values() if len(v) >= min_size]
    return sorted(out, key=len, reverse=True)

source_odf(env, cols)

Mean rectified log-envelope derivative over one source's partials.

Source code in src/ambiscape/rhythm.py

def source_odf(env: np.ndarray, cols) -> np.ndarray:
    """Mean rectified log-envelope derivative over one source's partials."""
    L = np.log10(env[:, cols] + EPS)
    return np.maximum(np.diff(L, axis=0, prepend=L[:1]), 0).mean(1)

pick_strikes(odf, t, min_sep=0.5, k=1.5, t_max=None)

Adaptive strongest-first onset picking.

Candidates are local maxima exceeding a running median by k MADs; they are accepted strongest-first subject to a min_sep guard (set it just below the shortest true inter-onset interval). Returns strike times.

Source code in src/ambiscape/rhythm.py

def pick_strikes(odf, t, min_sep=0.5, k=1.5, t_max=None):
    """Adaptive strongest-first onset picking.

    Candidates are local maxima exceeding a running median by ``k`` MADs;
    they are accepted strongest-first subject to a ``min_sep`` guard (set it
    just below the shortest true inter-onset interval). Returns strike times.
    """
    dt = float(np.median(np.diff(t)))
    med = median_filter(odf, int(8.0 / dt) | 1)
    mad = median_filter(np.abs(odf - med), int(8.0 / dt) | 1) + 1e-9
    z = (odf - med) / mad
    ismax = odf == maximum_filter1d(odf, 2 * int(0.3 / dt) + 1)
    ok = ismax & (z > k)
    if t_max is not None:
        ok &= t < t_max
    cand = np.flatnonzero(ok)
    taken = np.zeros(len(t), bool)
    guard = int(min_sep / dt)
    keep = []
    for i in cand[np.argsort(z[cand])[::-1]]:
        if not taken[max(0, i - guard):i + guard].any():
            keep.append(i)
            taken[i] = True
    return np.sort(t[np.array(keep, int)]) if keep else np.array([])

acf_structure(odf, dt, t_mask=None, max_lag_s=8.0, rel=0.25)

Cycle period and shortest intra-cycle gap from the ODF autocorrelation.

Returns (cycle_period, min_gap): the lag of the strongest ACF peak, and the shortest peak lag whose value exceeds rel x the strongest — use 0.8 * min_gap as the pick_strikes separation guard. Estimating the cycle from the ACF first keeps the later Rayleigh refinement off subharmonics.

Source code in src/ambiscape/rhythm.py

def acf_structure(odf, dt, t_mask=None, max_lag_s=8.0, rel=0.25):
    """Cycle period and shortest intra-cycle gap from the ODF autocorrelation.

    Returns (cycle_period, min_gap): the lag of the strongest ACF peak, and
    the shortest peak lag whose value exceeds ``rel`` x the strongest —
    use ``0.8 * min_gap`` as the ``pick_strikes`` separation guard. Estimating
    the cycle from the ACF first keeps the later Rayleigh refinement off
    subharmonics."""
    x = odf if t_mask is None else odf[t_mask]
    y = x - x.mean()
    a = np.correlate(y, y, "full")[len(y) - 1:]
    a /= a[0] + EPS
    m = int(max_lag_s / dt)
    pk, props = find_peaks(a[:m], prominence=0.02)
    if not len(pk):
        return None, None
    vals = a[pk]
    near_max = pk[vals >= 0.9 * vals.max()]     # prefer the fundamental over
    cycle = float(near_max.min() * dt)          # its multiples
    strong = pk[vals >= rel * vals.max()]
    return cycle, float(strong.min() * dt)

rayleigh_period(times, grid, harm=2)

Resultant length of the strike point process folded at each candidate period, summed over harm harmonics (multi-strike cycles).

Source code in src/ambiscape/rhythm.py

def rayleigh_period(times, grid, harm=2):
    """Resultant length of the strike point process folded at each candidate
    period, summed over ``harm`` harmonics (multi-strike cycles)."""
    R = np.zeros(len(grid))
    for i, P in enumerate(grid):
        ph = 2 * np.pi * times / P
        R[i] = sum(np.abs(np.exp(1j * h * ph).mean())
                   for h in range(1, harm + 1))
    return R

best_period(times, lo=0.5, hi=8.0, step=0.0005, harm=2)

Grid-search rayleigh_period with parabolic refinement.

Source code in src/ambiscape/rhythm.py

def best_period(times, lo=0.5, hi=8.0, step=5e-4, harm=2) -> float:
    """Grid-search ``rayleigh_period`` with parabolic refinement."""
    grid = np.arange(lo, hi, step)
    R = rayleigh_period(times, grid, harm)
    i = int(np.argmax(R))
    if 0 < i < len(R) - 1:
        d = (R[i - 1] - R[i + 1]) / (2 * (R[i - 1] - 2 * R[i] + R[i + 1]))
        return float(grid[i] + d * step)
    return float(grid[i])

period_track(times, P0, win=150.0, step=30.0, half_range=0.05, harm=2)

Sliding-window period estimates around P0; returns (t, P).

Source code in src/ambiscape/rhythm.py

def period_track(times, P0, win=150.0, step=30.0, half_range=0.05, harm=2):
    """Sliding-window period estimates around ``P0``; returns (t, P)."""
    ts, Ps = [], []
    g = np.arange(P0 - half_range, P0 + half_range, 5e-4)
    for w0 in np.arange(times.min(), times.max() - win, step):
        sel = times[(times >= w0) & (times < w0 + win)]
        if len(sel) < 20:
            continue
        ts.append(w0 + win / 2)
        Ps.append(g[np.argmax(rayleigh_period(sel, g, harm))])
    return np.array(ts), np.array(Ps)

phase_clusters(times, P, width=0.12)

Split one stream into phase clusters at period P (histogram modes). Returns (phase0, list of (center, mask)): phases are relative to the dominant cluster.

Source code in src/ambiscape/rhythm.py

def phase_clusters(times, P, width=0.12):
    """Split one stream into phase clusters at period ``P`` (histogram modes).
    Returns (phase0, list of (center, mask)): phases are relative to the
    dominant cluster."""
    ph = (times / P) % 1.0
    h, e = np.histogram(ph, bins=100)
    main = e[np.argmax(h)] + 0.005
    ph0 = (ph - main) % 1.0
    out = []
    left = np.ones(len(times), bool)
    while left.any():
        h, e = np.histogram(ph0[left], bins=50, range=(0, 1))
        if h.max() < max(3, 0.02 * len(times)):
            break
        c = e[np.argmax(h)] + 0.01
        d = np.minimum(np.abs(ph0 - c), 1 - np.abs(ph0 - c))
        m = left & (d < width)
        out.append((float(c) % 1.0, m))
        left &= ~m
    return ph0, out

cycle_grid(streams, P, t_max)

Per-cycle timing residuals and hit rates for named event streams sharing one cycle period.

streams maps name -> strike times. The cycle phase reference is the first stream. Returns per-stream position (s into cycle), residual array (NaN = missed cycle), and summary stats: timing sd, lag-1 autocorrelation, slow wander vs cycle-to-cycle sd, hit rate.

Source code in src/ambiscape/rhythm.py

def cycle_grid(streams: dict, P: float, t_max: float) -> dict:
    """Per-cycle timing residuals and hit rates for named event streams
    sharing one cycle period.

    ``streams`` maps name -> strike times. The cycle phase reference is the
    first stream. Returns per-stream position (s into cycle), residual array
    (NaN = missed cycle), and summary stats: timing sd, lag-1 autocorrelation,
    slow wander vs cycle-to-cycle sd, hit rate.
    """
    names = list(streams)
    ref = streams[names[0]]
    t0 = float(np.median((ref / P) % 1.0)) * P
    ncyc = int((t_max - t0) / P)
    out = {"P": P, "t0": t0, "ncyc": ncyc, "streams": {}}
    for name in names:
        tk = streams[name]
        pos = float(np.median(((tk - t0) / P) % 1.0)) * P
        res = np.full(ncyc, np.nan)
        for s in tk:
            c = int(round((s - t0 - pos) / P))
            if 0 <= c < ncyc:
                r = s - (t0 + c * P + pos)
                if abs(r) < 0.45 * P and (np.isnan(res[c]) or
                                          abs(r) < abs(res[c])):
                    res[c] = r
        v = res[~np.isnan(res)]
        m = ~np.isnan(res)
        if not m.any():          # cluster never landed on the cycle grid
            continue
        g = ~np.isnan(res[:-1]) & ~np.isnan(res[1:])
        r1 = float(np.corrcoef(res[:-1][g], res[1:][g])[0, 1]) if g.sum() > 3 \
            else np.nan
        xi = np.interp(np.arange(ncyc), np.flatnonzero(m), res[m])
        slow = uniform_filter1d(xi, max(3, int(60.0 / P)))
        out["streams"][name] = dict(
            pos=pos, res=res, hit_rate=float(m.mean()),
            sd_ms=float(np.std(v) * 1e3), lag1=r1,
            slow_sd_ms=float(np.std(slow) * 1e3),
            fast_sd_ms=float(np.std(xi - slow) * 1e3))
    return out

rise_spectrum(take, times, nfft=8192, n_max=150, seed=1)

Strike-triggered mean post/pre log-spectral rise (dB) on W.

The tonal fingerprint of one rhythmic position. If a stream's rise spectrum reproduces another source's partials, that stream is cross-talk rather than a distinct strike. Returns (freqs, rise_db).

Source code in src/ambiscape/rhythm.py

def rise_spectrum(take: Take, times, nfft=8192, n_max=150, seed=1):
    """Strike-triggered mean post/pre log-spectral rise (dB) on W.

    The tonal fingerprint of one rhythmic position. If a stream's rise
    spectrum reproduces another source's partials, that stream is cross-talk
    rather than a distinct strike. Returns (freqs, rise_db).
    """
    fs = take.samplerate
    win = np.hanning(nfft)
    freqs = np.fft.rfftfreq(nfft, 1 / fs)
    rng = np.random.default_rng(seed)
    sel = rng.choice(times, min(n_max, len(times)), replace=False)
    acc, n = np.zeros(len(freqs)), 0
    iW = take.wyzx[0]
    with sf.SoundFile(str(take.path)) as f:
        for s in sel:
            i = int(s * fs)
            if i < nfft or i + nfft > f.frames:
                continue
            f.seek(i - nfft + int(0.02 * fs))
            x = f.read(2 * nfft, dtype="float64", always_2d=True)[:, iW]
            if len(x) < 2 * nfft:
                continue
            pre = np.abs(np.fft.rfft(x[:nfft] * win)) ** 2
            post = np.abs(np.fft.rfft(x[nfft:] * win)) ** 2
            acc += 10 * np.log10((post + EPS) / (pre + EPS))
            n += 1
    return freqs, acc / max(n, 1)

run_session(sess, out_dir, n_sources=2, t_stop=None, verbose=True)

Full rhythm pipeline for a single-take session; writes rhythm_overview.png and rhythm.json, returns the summary dict.

Source code in src/ambiscape/rhythm.py

def run_session(sess, out_dir, n_sources=2, t_stop=None, verbose=True):
    """Full rhythm pipeline for a single-take session; writes
    ``rhythm_overview.png`` and ``rhythm.json``, returns the summary dict."""
    import json
    from .features import load_features
    out_dir = Path(out_dir)
    fdir = out_dir / "features"
    F = load_features(sorted(fdir.glob("*.npz")))
    take = sess.takes[0]
    active, quiet, med = _activity_masks(F)
    if t_stop is None:
        sm = median_filter(10 * np.log10(F["oct_pow"][:, 5] + EPS), 31)
        thr = (np.median(sm[np.repeat(active, 60)[:len(sm)]])
               + np.median(sm[np.repeat(quiet, 60)[:len(sm)]])) / 2
        t_stop = float(np.flatnonzero(sm > thr).max())
    if verbose:
        print(f"  active section ends at {t_stop:.0f} s")
    pfreq, rise = detect_partials(F, active, quiet)
    if verbose:
        print(f"  {len(pfreq)} partials "
              f"({pfreq.min():.0f}-{pfreq.max():.0f} Hz)")
    P = partial_pass(take, pfreq)
    t = P["t"]
    groups = cluster_partials(P["env"], mask=t < t_stop)[:n_sources]
    summary = {"t_stop_s": round(t_stop, 1), "sources": []}
    streams = {}
    for gi, cols in enumerate(groups):
        name = chr(ord("A") + gi)
        odf = source_odf(P["env"], cols)
        dt = float(np.median(np.diff(t)))
        cycle, min_gap = acf_structure(odf, dt, t_mask=t < t_stop)
        if cycle is None:
            continue
        s0 = pick_strikes(odf, t, min_sep=0.8 * min_gap, t_max=t_stop)
        if len(s0) < 20:
            continue
        Pbest = best_period(s0, lo=0.9 * cycle, hi=1.1 * cycle)
        ph0, clusters = phase_clusters(s0, Pbest)
        src = {"name": name,
               "partials_hz": [round(float(pfreq[c]), 1) for c in cols],
               "period_s": round(Pbest, 4),
               "n_strikes": int(len(s0)),
               "phase_clusters": [round(c, 3) for c, _ in clusters]}
        az, el, R = strike_doa(P, s0, cols)
        src["azimuth_deg"], src["elevation_deg"], src["az_R"] = \
            round(az, 1), round(el, 1), round(R, 2)
        for ci, (c, m) in enumerate(clusters):
            streams[f"{name}{ci}"] = s0[m]
        summary["sources"].append(src)
    # circular statistics per stream and inter-source phase locking
    if streams:
        from .circstats import phase_stats, relative_phase
        P0 = summary["sources"][0]["period_s"]
        summary["phase_stats"] = {n: phase_stats(tk, P0)
                                  for n, tk in streams.items()}
        prim = {n[0]: tk for n, tk in sorted(streams.items())
                if n.endswith("0")}
        names = sorted(prim)
        summary["phase_lock"] = {
            f"{b}_vs_{a}": relative_phase(prim[b], prim[a], P0)
            for a, b in zip(names, names[1:])}
    # flag phase clusters that coincide with another source's strikes
    # (cross-talk between partial groups, not independent strikes)
    for n, tk in streams.items():
        others = np.sort(np.concatenate(
            [v for m, v in streams.items() if m[0] != n[0]] or [np.array([])]))
        if not len(others) or not len(tk):
            continue
        i = np.clip(np.searchsorted(others, tk), 1, len(others) - 1)
        near = np.minimum(np.abs(tk - others[i - 1]), np.abs(tk - others[i]))
        frac = float((near < 0.06).mean())
        if frac > 0.5:
            summary.setdefault("crosstalk_suspects", {})[n] = round(frac, 2)
    if streams:
        P0 = summary["sources"][0]["period_s"]
        grid = cycle_grid(streams, P0, t_stop)
        summary["cycle"] = {
            n: {k: (round(v, 3) if isinstance(v, float) else v)
                for k, v in d.items() if k != "res"}
            for n, d in grid["streams"].items()}
        summary["cycle_period_s"] = P0
        _overview_figure(P, streams, grid, t_stop,
                         out_dir / "rhythm_overview.png", title=sess.name)
    (out_dir / "rhythm.json").write_text(json.dumps(summary, indent=2))
    return summary

partial_fm(take, freq, period, t_max, nfft=8192, hop=960)

Frequency modulation of one partial at the cycle rate.

A genuinely swinging bell Doppler-shifts its partials by a few cents at the swing period; a chimed (hammer-struck) bell does not. Tracks the instantaneous frequency of the partial (5-bin spectral centroid around its bin, energy-gated), then complex-demodulates at 1/period and at an off-rate control (1/(1.37*period)). Returns FM depth in cents at both; a swing verdict needs depth well above the control.

Source code in src/ambiscape/rhythm.py

def partial_fm(take: Take, freq: float, period: float, t_max: float,
               nfft=8192, hop=960) -> dict:
    """Frequency modulation of one partial at the cycle rate.

    A genuinely *swinging* bell Doppler-shifts its partials by a few cents
    at the swing period; a chimed (hammer-struck) bell does not. Tracks the
    instantaneous frequency of the partial (5-bin spectral centroid around
    its bin, energy-gated), then complex-demodulates at 1/period and at an
    off-rate control (1/(1.37*period)). Returns FM depth in cents at both;
    a swing verdict needs depth well above the control.
    """
    fs = take.samplerate
    hop = int(hop * fs / 48000)
    win = np.hanning(nfft).astype(np.float32)
    freqs = np.fft.rfftfreq(nfft, 1 / fs)
    b0 = int(round(freq * nfft / fs))
    sel = np.arange(b0 - 2, b0 + 3)
    iW = take.wyzx[0]
    fi, ei, tt = [], [], []
    carry = np.zeros((0,), np.float32)
    t_off = 0.0
    with sf.SoundFile(str(take.path)) as f:
        while t_off < t_max:
            block = f.read(60 * fs, dtype="float32", always_2d=True)[:, iW]
            if block.shape[0] == 0:
                break
            data = np.concatenate([carry, block]) if carry.shape[0] else block
            nwin = (len(data) - nfft) // hop + 1
            if nwin <= 0:
                carry = data
                continue
            idx = np.arange(nfft)[None, :] + hop * np.arange(nwin)[:, None]
            S = np.abs(np.fft.rfft(data[idx] * win))[:, sel] ** 2
            e = S.sum(1)
            fi.append((S * freqs[sel]).sum(1) / (e + EPS))
            ei.append(e)
            tt.append(t_off + (idx[:, 0] + nfft // 2) / fs)
            t_off += nwin * hop / fs
            carry = data[nwin * hop:].copy()
    fi = np.concatenate(fi)
    e = np.concatenate(ei)
    tt = np.concatenate(tt)
    m = (tt < t_max) & (e > np.percentile(e, 60))    # ringing frames only
    cents = 1200 * np.log2(fi[m] / freq)
    w = e[m] / e[m].sum()

    def demod(P):
        # weighted least-squares sinusoid fit: unbiased even though the
        # energy gate samples the cycle unevenly
        ph = 2 * np.pi * tt[m] / P
        A = np.stack([np.sin(ph), np.cos(ph), np.ones_like(ph)], 1)
        Aw = A * w[:, None]
        coef = np.linalg.lstsq(Aw.T @ A, Aw.T @ cents, rcond=None)[0]
        return float(np.hypot(coef[0], coef[1]))

    return {
        "freq_hz": freq, "period_s": period,
        "fm_cents_at_cycle": round(demod(period), 2),
        "fm_cents_control": round(demod(1.37 * period), 2),
        "n_frames": int(m.sum()),
    }

strike_doa(P, times, cols, dur=0.25)

Median per-strike azimuth/elevation (deg) from the pass arrays, energy-integrated over dur seconds after each strike.

Source code in src/ambiscape/rhythm.py

def strike_doa(P: dict, times, cols, dur=0.25):
    """Median per-strike azimuth/elevation (deg) from the pass arrays,
    energy-integrated over ``dur`` seconds after each strike."""
    t = P["t"]
    dt = float(np.median(np.diff(t)))
    az, el = [], []
    for s in times:
        sl = slice(int(s / dt), min(int(s / dt) + int(dur / dt), len(t)))
        Ix = P["ix"][sl][:, cols].sum()
        Iy = P["iy"][sl][:, cols].sum()
        Iz = P["iz"][sl][:, cols].sum()
        az.append(np.degrees(np.arctan2(Iy, Ix)))
        el.append(np.degrees(np.arctan2(Iz, np.hypot(Ix, Iy))))
    az = np.array(az)
    R = float(np.abs(np.exp(1j * np.radians(az)).mean()))
    return float(np.median(az)), float(np.median(el)), R

Spatial dynamics

Spatial dynamics at three time scales.

From the cached per-second spatial features (pseudo-intensity per octave, DOA, diffuseness) — no audio pass:

• direct_diffuse_split — per-octave directness (1 − diffuseness proxy) per second: the spatial analogue of foreground/background;
• passby_events — level events whose azimuth sweeps monotonically through the event: moving sources, with sweep rate and direction;
• azimuth_organization — windowed, energy-weighted circular concentration R(t): how directionally organized the scene is over time.

run_session writes spatial.json + spatial.png.

direct_diffuse_split(F)

Per-octave directness in [0, 1]: |pseudo-intensity| / band power.

Uses the cached I_band (re W*X etc. per octave) and oct_pow. A plane wave scores near 1, a diffuse field near 0. Returns (directness[nsec, nband], per-band medians).

Source code in src/ambiscape/spatial.py

def direct_diffuse_split(F: dict):
    """Per-octave directness in [0, 1]: |pseudo-intensity| / band power.

    Uses the cached ``I_band`` (re W*X etc. per octave) and ``oct_pow``.
    A plane wave scores near 1, a diffuse field near 0. Returns
    (directness[nsec, nband], per-band medians).
    """
    I = np.linalg.norm(F["I_band"], axis=2)
    d = np.clip(I / (F["oct_pow"] + EPS), 0, 1)
    return d, np.median(d, axis=0)

passby_events(F, min_dur_s=4, min_sweep_deg=25.0, min_r2=0.7)

Level events whose azimuth sweeps steadily: moving sources.

Detects events with :func:ambiscape.analysis.detect_events, then fits a line to the unwrapped per-second azimuth across each event lasting

= min_dur_s. A sweep of >= min_sweep_deg with fit R^2 >= min_r2 is a pass-by; the sweep sign gives the direction of travel (mic frame). Returns a list of dicts.

Source code in src/ambiscape/spatial.py

def passby_events(F: dict, min_dur_s=4, min_sweep_deg=25.0, min_r2=0.7):
    """Level events whose azimuth sweeps steadily: moving sources.

    Detects events with :func:`ambiscape.analysis.detect_events`, then fits
    a line to the unwrapped per-second azimuth across each event lasting
    >= ``min_dur_s``. A sweep of >= ``min_sweep_deg`` with fit R^2 >=
    ``min_r2`` is a pass-by; the sweep sign gives the direction of travel
    (mic frame). Returns a list of dicts.
    """
    from .analysis import detect_events
    dt = float(np.median(np.diff(F["t_fast"])))
    events, _bg = detect_events(F["fast_db"], dt)
    t0_abs = float(F["t"][0])
    out = []
    for e in events:
        a = F["t_fast"][e["i0"]] - t0_abs
        b = F["t_fast"][e["i1"]] - t0_abs
        i0, i1 = int(a), int(np.ceil(b))
        if i1 - i0 < min_dur_s or i1 >= len(F["az"]):
            continue
        az = np.unwrap(np.radians(F["az"][i0:i1]))
        x = np.arange(len(az), dtype=float)
        A = np.vstack([x, np.ones_like(x)]).T
        coef, res, *_ = np.linalg.lstsq(A, az, rcond=None)
        tot = ((az - az.mean()) ** 2).sum()
        r2 = 1 - float(res[0]) / (tot + EPS) if len(res) else 0.0
        sweep = float(np.degrees(coef[0]) * len(az))
        if abs(sweep) >= min_sweep_deg and r2 >= min_r2:
            out.append({
                "t0_s": i0, "dur_s": i1 - i0,
                "sweep_deg": round(sweep, 1),
                "rate_deg_s": round(float(np.degrees(coef[0])), 1),
                "direction": "left-to-right" if sweep < 0 else "right-to-left",
                "r2": round(r2, 2),
            })
    return out

azimuth_organization(F, win_s=60.0, step_s=20.0)

Windowed energy-weighted circular concentration of the azimuth.

Returns (t_centers, R): R near 1 = one dominant direction, near 0 = directionally disorganized. Window in seconds (per-second features).

Source code in src/ambiscape/spatial.py

def azimuth_organization(F: dict, win_s=60.0, step_s=20.0):
    """Windowed energy-weighted circular concentration of the azimuth.

    Returns (t_centers, R): R near 1 = one dominant direction, near 0 =
    directionally disorganized. Window in seconds (per-second features).
    """
    p = F["rms_w"].astype(np.float64) ** 2
    az = np.radians(F["az"])
    n, w, s = len(az), int(win_s), int(step_s)
    ts, Rs = [], []
    for i0 in range(0, n - w + 1, s):
        _mu, R = mean_resultant(az[i0:i0 + w], weights=p[i0:i0 + w])
        ts.append(float(F["t"][i0] + w / 2 - F["t"][0]))
        Rs.append(R)
    return np.array(ts), np.array(Rs)

directional_entropy(F, nbins=36)

Normalized Shannon entropy of the energy-weighted azimuth histogram.

"How many directions does this place sound from": 0 = all energy from one bearing, 1 = energy spread evenly around the horizon — the spatial analogue of an acoustic diversity index, and something only an ambisonic corpus can report.

Source code in src/ambiscape/spatial.py

def directional_entropy(F: dict, nbins: int = 36) -> float:
    """Normalized Shannon entropy of the energy-weighted azimuth histogram.

    "How many directions does this place sound from": 0 = all energy from
    one bearing, 1 = energy spread evenly around the horizon — the spatial
    analogue of an acoustic diversity index, and something only an
    ambisonic corpus can report.
    """
    p = np.asarray(F["rms_w"], np.float64) ** 2
    h, _ = np.histogram(F["az"], bins=nbins, range=(-180, 180), weights=p)
    q = h / (h.sum() + EPS)
    return float(-(q * np.log(q + EPS)).sum() / np.log(nbins))

horizon_fractions(F, limit_deg=10.0)

Energy fractions arriving from above / around / below the horizon.

Uses the per-second broadband DOA elevation, energy-weighted. A room heard from a couch splits mechanics on walls (above) from footsteps and structure-borne paths (level/below); outdoors it separates birds and building services from ground traffic.

Source code in src/ambiscape/spatial.py

def horizon_fractions(F: dict, limit_deg: float = 10.0) -> dict:
    """Energy fractions arriving from above / around / below the horizon.

    Uses the per-second broadband DOA elevation, energy-weighted. A room
    heard from a couch splits mechanics on walls (above) from footsteps
    and structure-borne paths (level/below); outdoors it separates birds
    and building services from ground traffic.
    """
    p = np.asarray(F["rms_w"], np.float64) ** 2
    el = np.asarray(F["el"], float)
    tot = p.sum() + EPS
    return {"above": round(float(p[el > limit_deg].sum() / tot), 2),
            "level": round(float(p[np.abs(el) <= limit_deg].sum() / tot), 2),
            "below": round(float(p[el < -limit_deg].sum() / tot), 2)}

fg_bg_az_overlap(F, nbins=36)

Bhattacharyya overlap of foreground vs background azimuth energy.

Foreground = loudest 25 % of seconds, background = quietest 25 % (the corpus convention). 1 = the foreground comes from where the background hums (one-source rooms), 0 = figure and ground occupy different directions (a street heard past a courtyard fountain).

Source code in src/ambiscape/spatial.py

def fg_bg_az_overlap(F: dict, nbins: int = 36) -> float:
    """Bhattacharyya overlap of foreground vs background azimuth energy.

    Foreground = loudest 25 % of seconds, background = quietest 25 % (the
    corpus convention). 1 = the foreground comes from where the background
    hums (one-source rooms), 0 = figure and ground occupy different
    directions (a street heard past a courtyard fountain).
    """
    p = np.asarray(F["rms_w"], np.float64) ** 2
    fg = p >= np.percentile(p, 75)
    bg = p <= np.percentile(p, 25)
    hists = []
    for m in (fg, bg):
        h, _ = np.histogram(F["az"][m], bins=nbins, range=(-180, 180),
                            weights=p[m])
        hists.append(h / (h.sum() + EPS))
    return float(np.sqrt(hists[0] * hists[1]).sum())

summarize_spatial(F)

Spatial descriptors for the analyze summary.

Source code in src/ambiscape/spatial.py

def summarize_spatial(F: dict) -> dict:
    """Spatial descriptors for the analyze summary."""
    hf = horizon_fractions(F)
    return {
        "directional_entropy": round(directional_entropy(F), 3),
        "above_horizon_fraction": hf["above"],
        "below_horizon_fraction": hf["below"],
        "fgbg_az_overlap": round(fg_bg_az_overlap(F), 2),
    }

run_session(sess, out_dir)

CLI driver: split + pass-bys + R(t), figure + spatial.json.

Source code in src/ambiscape/spatial.py

def run_session(sess, out_dir) -> dict:
    """CLI driver: split + pass-bys + R(t), figure + spatial.json."""
    import json
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from .features import load_features, OCT_CENTERS

    out_dir = Path(out_dir)
    F = load_features(sorted((out_dir / "features").glob("*.npz")))
    d, dmed = direct_diffuse_split(F)
    pb = passby_events(F)
    ts, Rs = azimuth_organization(F)
    doc = {
        "directness_median_per_octave": {
            str(int(c)): round(float(v), 2)
            for c, v in zip(OCT_CENTERS, dmed)},
        "azimuth_R_median": round(float(np.median(Rs)), 2),
        "azimuth_R_iqr": round(float(np.percentile(Rs, 75)
                                     - np.percentile(Rs, 25)), 2),
        "passbys": pb,
    }
    (out_dir / "spatial.json").write_text(json.dumps(doc, indent=2))

    fig, ax = plt.subplots(2, 1, figsize=(12.8, 6.4), dpi=130, sharex=True)
    tt = F["t"] - F["t"][0]
    ax[0].pcolormesh(tt, np.arange(len(OCT_CENTERS)), d.T, cmap="magma",
                     vmin=0, vmax=1, shading="auto")
    ax[0].set_yticks(range(len(OCT_CENTERS)),
                     [str(int(c)) for c in OCT_CENTERS], fontsize=7)
    ax[0].set(ylabel="octave (Hz)",
              title=f"{sess.name} — directness per octave (1=plane wave, "
                    "0=diffuse)")
    ax[1].plot(ts, Rs, color="#2a78d6", lw=1.2)
    for e in pb:
        ax[1].axvspan(e["t0_s"], e["t0_s"] + e["dur_s"], color="#d66a2a",
                      alpha=0.3)
    ax[1].set(xlabel="time (s)", ylabel="azimuth R (60 s)", ylim=(0, 1),
              title="directional organization; shaded = pass-by events")
    ax[1].grid(alpha=0.2)
    fig.tight_layout()
    fig.savefig(out_dir / "spatial.png")
    plt.close(fig)
    return doc

Schedule matching

Schedule matching: test event/strike streams against civic time grids.

Bells, chimes, and sirens follow wall-clock schedules — hourly strikes, quarter-hour chimes, fixed evening ringing. Given event times on the session's absolute clock, match_periods folds them at candidate civic periods and scores each with circular statistics; clock_offset turns a known schedule into a recorder-clock correction (the workflow behind clock_offset_s in calibration.json).

Times must be absolute seconds (session clock, i.e. take.start + offset into the take), otherwise grid phases are meaningless.

match_periods(times_abs, periods=CIVIC_PERIODS)

Fold events at each candidate period; score alignment.

Returns one dict per period — phase of the grid the events cluster on (seconds past the grid tick), R, circular SD, Rayleigh p — sorted by R. A meaningful match needs both high R and enough events spread over several grid cycles (n_cycles); R is trivially 1 when all events fall inside one cycle.

Source code in src/ambiscape/schedule.py

def match_periods(times_abs: np.ndarray, periods=CIVIC_PERIODS) -> list[dict]:
    """Fold events at each candidate period; score alignment.

    Returns one dict per period — phase of the grid the events cluster on
    (seconds past the grid tick), R, circular SD, Rayleigh p — sorted by R.
    A meaningful match needs both high R and enough events spread over
    several grid cycles (``n_cycles``); R is trivially 1 when all events
    fall inside one cycle.
    """
    t = np.asarray(times_abs, float)
    out = []
    for P in periods:
        ph = 2 * np.pi * (t / P % 1.0)
        mu, R = mean_resultant(ph)
        out.append({
            "period_s": P,
            "phase_s": round(float((mu / (2 * np.pi)) % 1.0 * P), 1),
            "R": round(R, 3),
            "circ_sd_s": round(circular_sd(R) / (2 * np.pi) * P, 1),
            "rayleigh_p": rayleigh_p(R, len(t)),
            "n": int(len(t)),
            "n_cycles": int(np.ptp(t) // P) + 1,
        })
    return sorted(out, key=lambda d: -d["R"])

grid_scan(F, period_s, phase_s=0.0, band=(300.0, 1500.0), win_s=120.0, min_rise_db=6.0, bg_win_s=300.0)

Targeted scan of every tick of a civic grid for band-limited strikes.

The complement of :func:match_periods: instead of asking which grid an event stream fits, look at each tick of a known grid (every quarter hour, every hour + phase_s) for energy in a band — a church clock in the bell band, whether or not the broadband detector heard it. Uses the cached features (band level above a running bg_win_s low-percentile background), so the scan is instant.

Returns one dict per tick inside the feature timeline: t_tick (absolute seconds), detected, rise_db (peak exceedance within win_s centered on the tick), and offset_s of that peak from the tick — a consistent nonzero offset across ticks is recorder-clock error (see :func:clock_offset).

Source code in src/ambiscape/schedule.py

def grid_scan(F: dict, period_s: float, phase_s: float = 0.0,
              band=(300.0, 1500.0), win_s: float = 120.0,
              min_rise_db: float = 6.0, bg_win_s: float = 300.0) -> list[dict]:
    """Targeted scan of every tick of a civic grid for band-limited strikes.

    The complement of :func:`match_periods`: instead of asking which grid an
    event stream fits, look *at each tick* of a known grid (every quarter
    hour, every hour + ``phase_s``) for energy in a band — a church clock in
    the bell band, whether or not the broadband detector heard it. Uses the
    cached features (band level above a running ``bg_win_s`` low-percentile
    background), so the scan is instant.

    Returns one dict per tick inside the feature timeline: ``t_tick``
    (absolute seconds), ``detected``, ``rise_db`` (peak exceedance within
    ``win_s`` centered on the tick), and ``offset_s`` of that peak from the
    tick — a consistent nonzero offset across ticks is recorder-clock error
    (see :func:`clock_offset`).
    """
    from scipy.ndimage import percentile_filter
    from .states import band_level

    t = np.asarray(F["t"], float)
    lvl = band_level(F, band)
    n = max(3, int(round(bg_win_s)) | 1)
    rise = lvl - percentile_filter(lvl, 10, size=n, mode="nearest")
    first = np.ceil((t[0] - phase_s) / period_s) * period_s + phase_s
    out = []
    for tick in np.arange(first, t[-1], period_s):
        m = np.abs(t - tick) <= win_s / 2
        if not m.any():
            continue
        i = int(np.argmax(rise[m]))
        r = float(rise[m][i])
        out.append({
            "t_tick": float(tick),
            "detected": bool(r >= min_rise_db),
            "rise_db": round(r, 1),
            "offset_s": round(float(t[m][i] - tick), 1),
        })
    return out

clock_offset(observed_abs, true_clock_s)

Recorder-clock correction from one event of known wall-clock time.

observed_abs is the event's time on the recorder clock (absolute session seconds), true_clock_s the known true time (seconds since midnight). Returns the clock_offset_s value for calibration.json (positive = recorder clock was slow).

Source code in src/ambiscape/schedule.py

def clock_offset(observed_abs: float, true_clock_s: float) -> float:
    """Recorder-clock correction from one event of known wall-clock time.

    ``observed_abs`` is the event's time on the recorder clock (absolute
    session seconds), ``true_clock_s`` the known true time (seconds since
    midnight). Returns the ``clock_offset_s`` value for
    ``calibration.json`` (positive = recorder clock was slow).
    """
    return float(true_clock_s - observed_abs % 86400.0)

run_session(sess, out_dir)

CLI driver: match cached event streams against civic periods.

Uses broadband events (always) and rhythm strikes when a prior ambiscape rhythm run left rhythm.json phase data; writes schedule.json.

Source code in src/ambiscape/schedule.py

def run_session(sess, out_dir) -> dict:
    """CLI driver: match cached event streams against civic periods.

    Uses broadband events (always) and rhythm strikes when a prior
    ``ambiscape rhythm`` run left ``rhythm.json`` phase data; writes
    ``schedule.json``.
    """
    import json
    from pathlib import Path
    from .analysis import detect_events
    from .features import load_features

    out_dir = Path(out_dir)
    F = load_features(sorted((out_dir / "features").glob("*.npz")))
    dt = float(np.median(np.diff(F["t_fast"])))
    events, _bg = detect_events(F["fast_db"], dt)
    t_ev = np.array([float(F["t_fast"][e["ipk"]]) for e in events])
    doc = {"events": match_periods(t_ev)[:4] if len(t_ev) >= 3 else []}
    (out_dir / "schedule.json").write_text(json.dumps(doc, indent=2))
    return doc

Event timbre templates

Event timbre templates: recurring event classes without machine learning.

Every transient event gets a spectral fingerprint — the strike-triggered post/pre rise spectrum (what appeared) plus a per-band decay slope (how it faded). Fingerprints are clustered by correlation distance into template classes: "the same sound again" across a whole session, fully transparent and corpus-comparable. Complements PANNs tagging ([ml]).

run_session fingerprints the session's spectral events (see :mod:background), clusters them, and writes timbre.json + timbre.png (class templates + counts + exemplar times).

event_fingerprint(take, t_onset, nfft=8192, decay_s=1.0)

Rise spectrum (dB, mel-ish bands) + per-band decay slope (dB/s).

Source code in src/ambiscape/timbre.py

def event_fingerprint(take, t_onset: float, nfft=8192, decay_s=1.0):
    """Rise spectrum (dB, mel-ish bands) + per-band decay slope (dB/s)."""
    fs = take.samplerate
    iW = take.wyzx[0]
    win = np.hanning(nfft)
    freqs = np.fft.rfftfreq(nfft, 1 / fs)
    W, centers = _melish_matrix(freqs)
    i = int(t_onset * fs)
    n_dec = int(decay_s * fs / nfft)
    with sf.SoundFile(str(take.path)) as f:
        if i < nfft or i + (n_dec + 1) * nfft > f.frames:
            return None, None, centers
        f.seek(i - nfft + int(0.02 * fs))
        x = f.read((n_dec + 2) * nfft, dtype="float64", always_2d=True)[:, iW]
    def spec(seg):
        return W @ (np.abs(np.fft.rfft(seg * win)) ** 2)
    pre = spec(x[:nfft])
    post = spec(x[nfft:2 * nfft])
    rise = 10 * np.log10((post + EPS) / (pre + EPS))
    tail = np.array([10 * np.log10(spec(x[(k + 1) * nfft:(k + 2) * nfft])
                                   + EPS) for k in range(n_dec + 1)])
    slope = np.polyfit(np.arange(n_dec + 1) * nfft / fs, tail, 1)[0]
    return rise, slope, centers

cluster_events(fps, th=0.35, min_size=2)

Average-linkage clustering of fingerprints by correlation distance. Returns labels (−1 = unclustered singleton).

Source code in src/ambiscape/timbre.py

def cluster_events(fps: np.ndarray, th=0.35, min_size=2):
    """Average-linkage clustering of fingerprints by correlation distance.
    Returns labels (−1 = unclustered singleton)."""
    from scipy.cluster.hierarchy import fcluster, linkage
    if len(fps) < 2:
        return np.zeros(len(fps), int)
    C = np.corrcoef(fps)
    d = np.clip(1 - C, 0, 2)
    np.fill_diagonal(d, 0)
    lab = fcluster(linkage(d[np.triu_indices_from(d, 1)], "average"),
                   th, criterion="distance")
    out = np.full(len(fps), -1)
    for l in np.unique(lab):
        m = lab == l
        if m.sum() >= min_size:
            out[m] = l
    return out

run_session(sess, out_dir, max_events=150)

Fingerprint + cluster the session's spectral events.

Source code in src/ambiscape/timbre.py

def run_session(sess, out_dir, max_events=150) -> dict:
    """Fingerprint + cluster the session's spectral events."""
    import json
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from .background import band_background, foreground, spectral_events
    from .features import load_features

    out_dir = Path(out_dir)
    F = load_features(sorted((out_dir / "features").glob("*.npz")))
    take = sess.takes[0]
    bg = band_background(F["logspec"])
    rise_db, _frac = foreground(F["logspec"], bg)
    ev = spectral_events(rise_db, F["logf"])
    # one fingerprint per onset: merge blobs that start within 1 s
    seen, dedup = set(), []
    for e in ev:
        if e["t0_s"] not in seen:
            dedup.append(e)
            seen.update((e["t0_s"] - 1, e["t0_s"], e["t0_s"] + 1))
    ev = dedup
    if len(ev) > max_events:
        keep = np.argsort([-e["peak_rise_db"] for e in ev])[:max_events]
        ev = [ev[i] for i in sorted(keep)]
    fps, slopes, kept = [], [], []
    for e in ev:
        r, s, centers = event_fingerprint(take, float(e["t0_s"]))
        if r is not None:
            fps.append(r)
            slopes.append(s)
            kept.append(e)
    fps = np.array(fps)
    lab = cluster_events(fps) if len(fps) else np.array([], int)
    classes = []
    for l in sorted(set(lab.tolist()) - {-1}):
        m = lab == l
        classes.append({
            "n": int(m.sum()),
            "exemplar_t0_s": [kept[i]["t0_s"]
                              for i in np.flatnonzero(m)[:5]],
            "centroid_hz": int(np.exp((fps[m].mean(0) * np.log(centers)).sum()
                                      / (fps[m].mean(0).sum() + EPS))),
            # decay only over the bands the event actually excited
            "decay_median_db_s": round(float(np.median(
                [np.median(s[r > 6]) for r, s in
                 zip(fps[m], np.array(slopes)[m]) if (r > 6).any()] or
                [np.nan])), 1),
        })
    classes.sort(key=lambda c: -c["n"])
    doc = {"n_events_fingerprinted": len(fps),
           "n_classes": len(classes),
           "n_unclustered": int((lab == -1).sum()),
           "classes": classes}
    (out_dir / "timbre.json").write_text(json.dumps(doc, indent=2))

    if len(fps):
        order = np.argsort(lab)
        fig, ax = plt.subplots(figsize=(11.2, 5.2), dpi=130)
        pc = ax.pcolormesh(np.arange(len(fps)), centers, fps[order].T,
                           cmap="magma", shading="auto")
        for b in np.flatnonzero(np.diff(lab[order]) != 0):
            ax.axvline(b + 0.5, color="w", lw=0.8)
        ax.set(yscale="log", xlabel="event (grouped by class)",
               ylabel="Hz", title=f"{sess.name} — event rise-spectrum "
               f"fingerprints, {len(classes)} classes")
        fig.colorbar(pc, ax=ax, pad=0.01, label="rise (dB)")
        fig.tight_layout()
        fig.savefig(out_dir / "timbre.png")
        plt.close(fig)
    return doc

MIR views (librosa)

Librosa-based tempogram and chromagram (optional [music] extra).

Complements the built-in analyses with the MIR-standard views: the tempogram (onset autocorrelation over time, in BPM — against which the windowed-ACF tempogram in :mod:rhythm can be cross-checked) and the chromagram (12-bin pitch-class energy over time, the time-resolved counterpart of :func:ambiscape.tonality.pitch_class_profile).

Audio is read from the W channel and resampled to 22.05 kHz; long sessions are fine (a 25 min file takes on the order of a minute). Requires pip install "ambiscape[music]".

load_w(take, t0=0.0, dur=None, sr=22050)

W channel of a take, resampled to sr.

Source code in src/ambiscape/music.py

def load_w(take, t0=0.0, dur=None, sr=22050):
    """W channel of a take, resampled to ``sr``."""
    librosa = _require_librosa()
    fs = take.samplerate
    iW = take.wyzx[0]
    with sf.SoundFile(str(take.path)) as f:
        f.seek(int(t0 * fs))
        n = f.frames - int(t0 * fs) if dur is None else int(dur * fs)
        x = f.read(n, dtype="float32", always_2d=True)[:, iW]
    return librosa.resample(x, orig_sr=fs, target_sr=sr), sr

tempogram(y, sr, hop=512, win_s=8.0)

Autocorrelation tempogram of the onset-strength envelope.

Returns (times, bpm_axis, T, tempo_bpm): the tempogram plus librosa's global tempo estimate (which resolves the octave ambiguity a raw tempogram argmax suffers from).

Source code in src/ambiscape/music.py

def tempogram(y, sr, hop=512, win_s=8.0):
    """Autocorrelation tempogram of the onset-strength envelope.

    Returns (times, bpm_axis, T, tempo_bpm): the tempogram plus librosa's
    global tempo estimate (which resolves the octave ambiguity a raw
    tempogram argmax suffers from).
    """
    librosa = _require_librosa()
    onset = librosa.onset.onset_strength(y=y, sr=sr, hop_length=hop)
    win = int(win_s * sr / hop)
    T = librosa.feature.tempogram(onset_envelope=onset, sr=sr,
                                  hop_length=hop, win_length=win)
    times = librosa.frames_to_time(np.arange(T.shape[1]), sr=sr,
                                   hop_length=hop)
    bpm = librosa.tempo_frequencies(T.shape[0], sr=sr, hop_length=hop)
    try:
        from librosa.feature.rhythm import tempo as _tempo
    except ImportError:                      # librosa < 0.10
        _tempo = librosa.beat.tempo
    t_est = float(_tempo(onset_envelope=onset, sr=sr, hop_length=hop)[0])
    return times, bpm, T, t_est

chromagram(y, sr, hop=512)

STFT chromagram (12 pitch classes over time).

Source code in src/ambiscape/music.py

def chromagram(y, sr, hop=512):
    """STFT chromagram (12 pitch classes over time)."""
    librosa = _require_librosa()
    C = librosa.feature.chroma_stft(y=y, sr=sr, hop_length=hop)
    times = librosa.frames_to_time(np.arange(C.shape[1]), sr=sr,
                                   hop_length=hop)
    return times, C

run_session(sess, out_dir, t0=0.0, dur=None)

Tempogram + chromagram figure and summary for the first take.

Source code in src/ambiscape/music.py

def run_session(sess, out_dir, t0=0.0, dur=None) -> dict:
    """Tempogram + chromagram figure and summary for the first take."""
    import json
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    librosa = _require_librosa()
    from .tonality import NOTE

    out_dir = Path(out_dir)
    y, sr = load_w(sess.takes[0], t0=t0, dur=dur)
    tt, bpm, T, bpm_peak = tempogram(y, sr)
    tc, C = chromagram(y, sr)
    chroma_mean = C.mean(1)
    doc = {
        "tempo_bpm_global": round(bpm_peak, 1),
        "tempo_period_s": round(60.0 / bpm_peak, 3),
        "chroma_mean": {NOTE[i]: round(float(v / chroma_mean.sum()), 3)
                        for i, v in enumerate(chroma_mean)},
        "top_pitch_classes": [NOTE[i]
                              for i in np.argsort(chroma_mean)[::-1][:3]],
    }
    (out_dir / "music.json").write_text(json.dumps(doc, indent=2))

    fig, ax = plt.subplots(2, 1, figsize=(12.8, 7.2), dpi=130, sharex=True)
    bmask = (bpm > 5) & (bpm < 300)
    ax[0].pcolormesh(tt, bpm[bmask], T[bmask], cmap="magma", shading="auto")
    ax[0].set(yscale="log", ylabel="tempo (BPM)",
              title=f"{sess.name} — tempogram (librosa); global peak "
                    f"{bpm_peak:.1f} BPM = {60/bpm_peak:.2f} s")
    ax[1].pcolormesh(tc, np.arange(12), C, cmap="magma", shading="auto")
    ax[1].set_yticks(range(12), NOTE, fontsize=7)
    ax[1].set(xlabel="time (s)", ylabel="pitch class",
              title="chromagram (librosa)")
    fig.tight_layout()
    fig.savefig(out_dir / "music.png")
    plt.close(fig)
    return doc

Taxonomy

Schaeffer / Schafer taxonomy figures from per-session annotations.

The annotation file (annotations.json or .yml in the session folder) is hand-authored: instruments detect when things sound, but assigning a sound to Schaeffer's typo-morphology or Schafer's functional categories is an interpretive act. This module turns that interpretation into two figures:

• schaeffer_map — objects on the facture x mass plane, colored by Schafer function (keynote / signal / soundmark / figure);
• schafer_timeline — one lane per object on the session clock, keynote spans as bars, events as markers, lo-fi states shaded.

Annotation schema (JSON; YAML accepted if PyYAML is installed)::

{
  "objects": [
    {"name": "air-pump drone",
     "label": "air-pump drone (130 Hz comb, 9 h)",   # optional
     "kind": "keynote",             # keynote|signal|soundmark|figure
     "soundmark": "dwelling",       # optional: community|dwelling
     "source": "anthropophony",     # optional: ...|biophony|geophony
     "mass": "noise",               # tonic|tonic-complex|complex|noise
     "facture": "unlimited",        # impulse|iteration|sustained|unlimited
     "spans": [["23:01:36", "1 07:53:55"]],   # and/or
     "events": ["1 04:42:51"]},
    ...
  ],
  "states": [
    {"label": "LO-FI (drone masks the field)",
     "span": ["23:01:36", "1 07:53:55"]}
  ]
}

Times are "[D ]HH:MM:SS" where the optional leading integer D is days after the session's first day (or plain seconds as a number).

schaeffer_map(ann, out_path, title='')

Objects on the facture x mass grid, colored by Schafer function.

Source code in src/ambiscape/taxonomy.py

def schaeffer_map(ann: dict, out_path, title=""):
    """Objects on the facture x mass grid, colored by Schafer function."""
    with plt.rc_context(RC):
        fig, ax = plt.subplots(figsize=(9.6, 6.4), dpi=130)
        ax.grid(False)
        for i in range(5):
            ax.axhline(i - 0.5, color=GRID, lw=0.8, zorder=0)
            ax.axvline(i - 0.5, color=GRID, lw=0.8, zorder=0)
        cells: dict[tuple, list] = {}
        for o in ann["objects"]:
            key = (FACTURES.index(o["facture"]), MASSES.index(o["mass"]))
            cells.setdefault(key, []).append(o)
        offsets = [(0, .1), (-.18, -.12), (.18, -.12), (-.18, .3), (.18, .3)]
        kinds_seen, bio_seen, ring_seen = set(), False, False
        for (x, y), objs in cells.items():
            for o, (dx, dy) in zip(objs, offsets):
                ring = "soundmark" in o and o["kind"] != "soundmark"
                ring_seen |= ring
                bio_seen |= o.get("source") == "biophony"
                kinds_seen.add(o["kind"])
                ax.scatter(x + dx, y + dy, s=170, marker=_marker(o),
                           color=KIND_COLOR[o["kind"]], zorder=3,
                           edgecolors=MAGENTA if ring else "none",
                           linewidths=2.2)
                ax.annotate(o.get("label", o["name"]), (x + dx, y + dy),
                            xytext=(0, -15), ha="center",
                            textcoords="offset points", fontsize=8.3,
                            color=INK, zorder=4)
        ax.set_xticks(range(4), FACTURE_LABELS)
        ax.set_yticks(range(4), MASS_LABELS)
        ax.set_xlim(-0.5, 3.5)
        ax.set_ylim(3.5, -0.5)
        ax.set_xlabel("facture / temporal sustainment  (Schaeffer typology) →")
        ax.set_ylabel("← mass  (Schaeffer morphology)")
        ax.set_title(f"{title} — sound objects in Schaeffer's typo-morphology,"
                     " colored by Schafer function", loc="left", fontsize=10.5)
        names = {"keynote": "keynote (ground)", "signal": "signal (figure)",
                 "soundmark": "community soundmark",
                 "figure": "incidental figure"}
        handles = [Line2D([], [], marker="o", ls="none", color=KIND_COLOR[k],
                          label=names[k]) for k in names if k in kinds_seen]
        if ring_seen:
            handles.append(Line2D([], [], marker="o", ls="none", color=SURF,
                                  markeredgecolor=MAGENTA, markeredgewidth=2,
                                  label="dwelling soundmark (ring)"))
        if bio_seen:
            handles.append(Line2D([], [], marker="^", ls="none", color=GREEN,
                                  label="biophony (triangle)"))
        ax.legend(handles=handles, loc="lower left", frameon=False,
                  fontsize=8, ncol=2)
        fig.tight_layout()
        fig.savefig(out_path, bbox_inches="tight")
        plt.close(fig)

schafer_timeline(ann, out_path, title='', session=None)

Lane timeline of annotated objects; lo-fi states shaded.

Source code in src/ambiscape/taxonomy.py

def schafer_timeline(ann: dict, out_path, title="", session=None):
    """Lane timeline of annotated objects; lo-fi states shaded."""
    objects = ann["objects"]
    states = ann.get("states", [])
    lanes = ["state"] + [o["name"] for o in objects] if states else \
            [o["name"] for o in objects]
    panels = _panels(ann, session)
    with plt.rc_context(RC):
        fig, axes = plt.subplots(
            1, len(panels), figsize=(12.8, 0.52 * len(lanes) + 1.6), dpi=130,
            sharey=True, squeeze=False,
            gridspec_kw={"width_ratios": [b - a for a, b in panels],
                         "wspace": 0.03})
        axes = axes[0]
        ny = len(lanes)

        def Y(name):
            return ny - 1 - lanes.index(name)

        for ax, (t0, t1) in zip(axes, panels):
            ax.grid(False)
            for st in states:
                a, b = (parse_time(x) for x in st["span"])
                a, b = max(a, t0), min(b, t1)
                if a >= b:
                    continue
                ax.axvspan(a, b, color=LOFI, zorder=0)
                ax.add_patch(Rectangle((a, Y("state") - 0.3), b - a, 0.6,
                                       color=YELLOW, alpha=0.55, lw=0))
                ax.annotate(st.get("label", "lo-fi"), ((a + b) / 2, Y("state")),
                            ha="center", va="center", fontsize=8,
                            color="#6b4a00")
            for o in objects:
                y = Y(o["name"])
                c = KIND_COLOR[o["kind"]]
                for a, b in o.get("spans", []):
                    a, b = parse_time(a), parse_time(b)
                    a, b = max(a, t0), min(b, t1)
                    if a >= b:
                        continue
                    ax.add_patch(Rectangle((a, y - 0.2), max(b - a, (t1-t0)*0.004),
                                           0.4, color=c,
                                           alpha=0.85 if o["kind"] == "keynote"
                                           else 1.0, lw=0))
                ev = [parse_time(e) for e in o.get("events", [])]
                ev = [e for e in ev if t0 <= e <= t1]
                if ev:
                    mk = _marker(o)
                    if o["kind"] == "figure":
                        mk = "|"
                    ax.plot(ev, [y] * len(ev), ls="none", marker=mk,
                            ms=11 if mk == "|" else 7, mew=1.8, color=c)
            ax.set_ylim(-0.7, ny - 0.3)
            ax.set_xlim(t0, t1)
            span = t1 - t0
            step = 3600 if span > 5400 else (600 if span > 900 else 120)
            ticks = np.arange(np.ceil(t0 / step) * step, t1, step)
            ax.set_xticks(ticks)
            ax.set_xticklabels([f"{int(x % 86400)//3600:02d}:"
                                f"{int(x % 3600)//60:02d}" for x in ticks])
            for sp in ("top", "right"):
                ax.spines[sp].set_visible(False)
            ax.grid(axis="x", color=GRID, lw=0.5)
            if session is not None and len(panels) > 1:
                ax.set_title(session.clock(t0)[:6], loc="left",
                             fontsize=8.5, color=SEC)
        axes[0].set_yticks(range(ny), lanes[::-1])
        for lab in axes[0].get_yticklabels():
            o = next((o for o in objects if o["name"] == lab.get_text()), None)
            if o and o["kind"] == "keynote":
                lab.set_color("#1c5cab")
        fig.suptitle(f"{title} — Schafer soundscape timeline: keynotes (blue "
                     "lanes), signals (green), soundmarks (magenta), "
                     "incidental figures (grey)", x=0.01, ha="left",
                     fontsize=10.5, color=INK)
        fig.tight_layout(rect=(0, 0, 1, 0.96))
        fig.savefig(out_path, bbox_inches="tight")
        plt.close(fig)

render(folder, out_dir=None, session=None)

Load annotations from a session folder and write both figures.

Source code in src/ambiscape/taxonomy.py

def render(folder: str | Path, out_dir=None, session=None):
    """Load annotations from a session folder and write both figures."""
    folder = Path(folder)
    ann = load_annotations(folder)
    out = Path(out_dir) if out_dir else folder / "analysis"
    out.mkdir(parents=True, exist_ok=True)
    if session is None:
        from .io import open_session
        try:
            session = open_session(folder)
        except (FileNotFoundError, ValueError):
            session = None
    name = folder.name
    schaeffer_map(ann, out / "schaeffer_map.png", title=name)
    schafer_timeline(ann, out / "schafer_timeline.png", title=name,
                     session=session)
    return out / "schaeffer_map.png", out / "schafer_timeline.png"

Draft annotations

Draft-annotation generator.

Pre-fills annotations.draft.json for a session from its cached features: steady level regimes become keynote candidates (with spans), detected transient events become one unclassified "figure" object whose entries carry listening hints (clock time, exceedance, azimuth/elevation, diffuseness). The taxonomy fields that require a human ear (mass, facture, kind, soundmark status) are left as "TODO" — edit, rename, delete, then save the result as annotations.json and run ambiscape taxonomy.

Calibration and ISO indicators

ISO 12913-3-style psychoacoustic indicators + level calibration.

Calibration

A session is calibrated by <folder>/calibration.json::

{"dbfs_to_dbspl": 94.0,
 "method": "SPL app next to mic, air pump running, LAeq 42 dB",
 "date": "2026-07-16"}

dbfs_to_dbspl is the offset O such that a signal at −X dBFS corresponds to (O − X) dB SPL. With it, dBFS descriptors become dB SPL and waveforms convert to pascals for psychoacoustic metrics.

The same file may carry clock_offset_s — seconds added to every take's start time when the recorder clock was found to be off (positive = clock was slow; e.g. calibrated against a known external time reference). Applied in :func:ambiscape.io.open_session, so all clock-labeled outputs agree. Both keys are optional.

Indicators (via MoSQITo, optional dependency)

ISO 532-1 time-varying loudness (N5, N50), DIN 45692 sharpness, and Daniel & Weber roughness, computed per ear on a binaural render of the B-format signal. If ambiviz (with its HRIR-based binauralizer) is installed it is used; otherwise a documented fallback renders a back-to-back cardioid pair at ±90° — a pseudo-binaural approximation without pinna/ILD spectral cues. Uncalibrated sessions are computed with an assumed offset and flagged: absolute sone/acum values are then indicative only (their ratios between segments remain meaningful).

apply_calibration(summary, cal)

Add dB SPL versions of the level descriptors to a summary dict.

Source code in src/ambiscape/iso.py

def apply_calibration(summary: dict, cal: dict) -> dict:
    """Add dB SPL versions of the level descriptors to a summary dict."""
    off = float(cal["dbfs_to_dbspl"])
    out = dict(summary)
    for key in ("leq_dbfs", "laeq_dbfs", "L10", "L50", "L90"):
        if key in summary and summary[key] is not None:
            out[key.replace("_dbfs", "") + "_db_spl"] = round(summary[key] + off, 1)
    out["calibration"] = {"dbfs_to_dbspl": off,
                          "method": cal.get("method", "")}
    return out

binaural(x4, fs)

FOA (W,Y,Z,X columns) -> stereo ear signals.

Tries ambiviz's HRIR binauralizer; falls back to a ±90° cardioid pair (no pinna cues). Returns (n, 2) array and the method name.

Source code in src/ambiscape/iso.py

def binaural(x4: np.ndarray, fs: int) -> tuple[np.ndarray, str]:
    """FOA (W,Y,Z,X columns) -> stereo ear signals.

    Tries ambiviz's HRIR binauralizer; falls back to a ±90° cardioid pair
    (no pinna cues). Returns (n, 2) array and the method name.
    """
    try:
        from ambiviz.ambisonics.binauralizer import binauralize  # type: ignore
        y = binauralize(x4.T, fs)  # ambiviz convention: channels first
        return np.asarray(y).T[:, :2], "ambiviz-hrir"
    except Exception:
        w, ych = x4[:, 0], x4[:, 1]
        left = 0.5 * (w + ych)
        right = 0.5 * (w - ych)
        return np.stack([left, right], axis=1), "cardioid-pair-fallback"

indicators(x_pa, fs, rough_dur=10.0)

ISO 532-1 loudness (N5/N50), DIN 45692 sharpness, D&W roughness for one calibrated (pascal) channel.

MoSQITo runs ~5x slower than realtime, and roughness is the costliest metric; it is therefore computed on a central rough_dur-second slice (roughness is a texture measure and stabilizes within seconds).

Source code in src/ambiscape/iso.py

def indicators(x_pa: np.ndarray, fs: int, rough_dur: float = 10.0) -> dict:
    """ISO 532-1 loudness (N5/N50), DIN 45692 sharpness, D&W roughness
    for one calibrated (pascal) channel.

    MoSQITo runs ~5x slower than realtime, and roughness is the costliest
    metric; it is therefore computed on a central `rough_dur`-second slice
    (roughness is a texture measure and stabilizes within seconds).
    """
    from mosqito.sq_metrics import (loudness_zwtv,
                                    sharpness_din_from_loudness,
                                    roughness_dw)
    N, N_spec, _bark, _t = loudness_zwtv(x_pa, fs, field_type="diffuse")
    S = sharpness_din_from_loudness(N, N_spec)
    n_r = int(rough_dur * fs)
    mid = max(0, (len(x_pa) - n_r) // 2)
    R = np.atleast_1d(roughness_dw(x_pa[mid:mid + n_r], fs)[0])
    return {
        "N5_sone": round(float(np.percentile(N, 95)), 2),
        "N50_sone": round(float(np.percentile(N, 50)), 2),
        "sharpness_median_acum": round(float(np.median(S)), 2),
        "roughness_median_asper": round(float(np.median(R)), 3),
    }

segment_indicators(sess, F, folder, dur=30.0, offset=None)

Compute per-ear indicators on representative segments.

Segments come from analysis.pick_segments (typical / quietest / most_active / transition); dur seconds from the start of each.

Source code in src/ambiscape/iso.py

def segment_indicators(sess, F: dict, folder: str | Path,
                       dur: float = 30.0, offset: float | None = None) -> dict:
    """Compute per-ear indicators on representative segments.

    Segments come from analysis.pick_segments (typical / quietest /
    most_active / transition); `dur` seconds from the start of each.
    """
    from .analysis import pick_segments
    from .io import read_span

    cal = load_calibration(folder)
    has_spl = bool(cal and "dbfs_to_dbspl" in cal)
    if offset is None:
        offset = float(cal["dbfs_to_dbspl"]) if has_spl else ASSUMED_OFFSET
    calibrated = has_spl or offset != ASSUMED_OFFSET

    out = {"calibrated": calibrated, "dbfs_to_dbspl": offset,
           "field_type": "diffuse", "segments": {}}
    if not calibrated:
        out["warning"] = (f"no calibration.json — assumed offset "
                          f"{ASSUMED_OFFSET} dB; absolute values indicative only")
    for pick in pick_segments(F, seg_s=dur):
        try:
            x, fs = read_span(sess, pick["t0"], dur)
        except ValueError:
            continue
        ears, method = binaural(x, fs)
        seg = {"t0": sess.clock(pick["t0"]), "dur_s": dur,
               "binaural_method": method}
        for ch, name in ((0, "left"), (1, "right")):
            seg[name] = indicators(to_pascal(ears[:, ch], offset), fs)
        seg["N5_sone_max_ear"] = max(seg["left"]["N5_sone"],
                                     seg["right"]["N5_sone"])
        out["segments"][pick["kind"]] = seg
    return out

room_criteria(oct_spl_db)

NR, NC, and RC ratings of an octave-band SPL spectrum.

oct_spl_db maps octave center frequency (Hz) to band SPL (dB). Ratings are only physically meaningful for calibrated levels (dbfs_to_dbspl in calibration.json); on uncalibrated dBFS they are relative numbers, comparable within one recorder+gain setup only.

• NR (ISO/R 1996 Noise Rating): analytic curves L = a + b*NR; the rating is the highest per-band NR value and NR_governing_hz names the band that sets it.
• NC (ANSI S12.2 Noise Criterion): tangency against the tabulated curves, linearly interpolated per band (63 Hz–8 kHz).
• RC (Blazier Room Criterion, simplified): arithmetic mean of the 500/1000/2000 Hz levels; the reference line has a −5 dB/octave slope through (1 kHz, RC). RC_class is "R" (rumble) when any 31.5–250 Hz band exceeds the line by > 5 dB, "H" (hiss) when any 2–4 kHz band exceeds it by > 3 dB, "RH" for both, "N" (neutral) otherwise.

Source code in src/ambiscape/iso.py

def room_criteria(oct_spl_db: dict) -> dict:
    """NR, NC, and RC ratings of an octave-band SPL spectrum.

    ``oct_spl_db`` maps octave center frequency (Hz) to band SPL (dB).
    Ratings are only physically meaningful for *calibrated* levels
    (``dbfs_to_dbspl`` in ``calibration.json``); on uncalibrated dBFS they
    are relative numbers, comparable within one recorder+gain setup only.

    - **NR** (ISO/R 1996 Noise Rating): analytic curves ``L = a + b*NR``;
      the rating is the highest per-band NR value and
      ``NR_governing_hz`` names the band that sets it.
    - **NC** (ANSI S12.2 Noise Criterion): tangency against the tabulated
      curves, linearly interpolated per band (63 Hz–8 kHz).
    - **RC** (Blazier Room Criterion, simplified): arithmetic mean of the
      500/1000/2000 Hz levels; the reference line has a −5 dB/octave slope
      through (1 kHz, RC). ``RC_class`` is "R" (rumble) when any
      31.5–250 Hz band exceeds the line by > 5 dB, "H" (hiss) when any
      2–4 kHz band exceeds it by > 3 dB, "RH" for both, "N" (neutral)
      otherwise.
    """
    spec = {float(k): float(v) for k, v in oct_spl_db.items()}

    nr_per = {f: (spec[f] - NR_A[f]) / NR_B[f] for f in NR_A if f in spec}
    f_gov = max(nr_per, key=nr_per.get)
    nr = nr_per[f_gov]

    nc = None
    ncs = sorted(NC_TABLE)
    per_band = []
    for i, f in enumerate(NC_FREQS):
        if f not in spec:
            continue
        levels = np.array([NC_TABLE[n][i] for n in ncs], float)
        per_band.append(float(np.interp(spec[f], levels, ncs)))
    if per_band:
        nc = max(per_band)

    rc = None
    rc_class = None
    if all(f in spec for f in (500.0, 1000.0, 2000.0)):
        rc = (spec[500.0] + spec[1000.0] + spec[2000.0]) / 3
        ref = {f: rc + 5 * np.log2(1000.0 / f) for f in spec}
        rumble = any(spec[f] > ref[f] + 5 for f in (31.5, 63.0, 125.0, 250.0)
                     if f in spec)
        hiss = any(spec[f] > ref[f] + 3 for f in (2000.0, 4000.0)
                   if f in spec)
        rc_class = ("RH" if rumble and hiss else
                    "R" if rumble else "H" if hiss else "N")

    return {"NR": round(nr, 1), "NR_governing_hz": int(f_gov),
            "NC": round(nc, 1) if nc is not None else None,
            "RC": round(rc, 1) if rc is not None else None,
            "RC_class": rc_class}

background_octaves_db(F, pct=50.0, offset_db=0.0)

Per-octave percentile level (dB) from cached features, for :func:room_criteria. offset_db is the dBFS→dB SPL calibration offset (0 keeps uncalibrated dBFS).

Source code in src/ambiscape/iso.py

def background_octaves_db(F: dict, pct: float = 50.0,
                          offset_db: float = 0.0) -> dict:
    """Per-octave percentile level (dB) from cached features, for
    :func:`room_criteria`. ``offset_db`` is the dBFS→dB SPL calibration
    offset (0 keeps uncalibrated dBFS)."""
    from .features import OCT_CENTERS
    lv = 10 * np.log10(np.asarray(F["oct_pow"], float) + 1e-20) + offset_db
    return {c: float(np.percentile(lv[:, i], pct))
            for i, c in enumerate(OCT_CENTERS) if c <= 8000}

Machine listening

Machine-listening helpers (optional [ml] extra).

• PANNs (CNN14, AudioSet, 527 classes) tags 10-s windows around detected events and steady states; used by ambiscape draft to suggest object names in annotations.draft.json.
• silero-vad estimates the fraction of speech in a file or span — the privacy gate to run on every excerpt before publishing (Freesound etc.).
• BirdNET (birdnetlib) identifies bird species in 3-s windows — the species layer for biophony, best run on the hi-fi windows the drone-free soundscape exposes (see :mod:ambiscape.biophony for the no-ML structural measures it confirms).

All models are trained on 16/32 kHz mono internet audio: the W (omni) channel is downmixed and resampled, spatial information is not used, and low-SNR domestic material is out of distribution — treat tags as suggestions to confirm by ear, not ground truth.

tag_window(x, fs, top_k=3, min_prob=0.1)

AudioSet tags for one mono window via PANNs CNN14 (32 kHz input).

Source code in src/ambiscape/ml.py

def tag_window(x: np.ndarray, fs: int, top_k: int = 3,
               min_prob: float = 0.10) -> list[dict]:
    """AudioSet tags for one mono window via PANNs CNN14 (32 kHz input)."""
    global _panns_model
    from panns_inference import AudioTagging, labels
    if _panns_model is None:
        _panns_model = AudioTagging(checkpoint_path=None, device="cpu")
    y = _resample(x.astype(np.float32), fs, 32000)
    clip = np.clip(y, -1, 1)[None, :]
    clipwise, _emb = _panns_model.inference(clip)
    probs = np.asarray(clipwise)[0]
    order = np.argsort(probs)[::-1][:top_k]
    return [{"label": labels[i], "p": round(float(probs[i]), 2)}
            for i in order if probs[i] >= min_prob]

birdnet_window(x, fs, lat=None, lon=None, week=-1, min_conf=0.25)

BirdNET species detections for one mono window (48 kHz input).

Analyzes the W channel resampled to 48 kHz. lat/lon and week (1–48, ISO-ish) enable BirdNET's location/season species filter — pass the session's coordinates to cut false positives. Returns [{"species", "common_name", "confidence"}] above min_conf. Requires the [ml] extra plus birdnetlib.

Source code in src/ambiscape/ml.py

def birdnet_window(x: np.ndarray, fs: int, lat: float | None = None,
                   lon: float | None = None, week: int = -1,
                   min_conf: float = 0.25) -> list[dict]:
    """BirdNET species detections for one mono window (48 kHz input).

    Analyzes the W channel resampled to 48 kHz. ``lat``/``lon`` and
    ``week`` (1–48, ISO-ish) enable BirdNET's location/season species
    filter — pass the session's coordinates to cut false positives.
    Returns ``[{"species", "common_name", "confidence"}]`` above
    ``min_conf``. Requires the ``[ml]`` extra plus ``birdnetlib``.
    """
    global _birdnet_analyzer
    import tempfile
    import soundfile as sf
    from birdnetlib import Recording
    from birdnetlib.analyzer import Analyzer
    if _birdnet_analyzer is None:
        _birdnet_analyzer = Analyzer()
    y = _resample(x.astype(np.float32), fs, 48000)
    with tempfile.NamedTemporaryFile(suffix=".wav", delete=True) as tmp:
        sf.write(tmp.name, np.clip(y, -1, 1), 48000)
        kw = dict(min_conf=min_conf, week_48=week)
        if lat is not None and lon is not None:
            kw.update(lat=lat, lon=lon)
        rec = Recording(_birdnet_analyzer, tmp.name, **kw)
        rec.analyze()
    return [{"species": d["scientific_name"],
             "common_name": d["common_name"],
             "confidence": round(float(d["confidence"]), 2)}
            for d in rec.detections]

birdnet_session(sess, F=None, windows=None, win_s=9.0, hifi_max_diffuse=None, lat=None, lon=None, min_conf=0.25)

Run BirdNET across a session, optionally only on hi-fi windows.

windows is an explicit list of absolute start seconds; if omitted, the session is tiled in win_s steps. When F (cached features) and hifi_max_diffuse are given, windows whose median diffuseness exceeds the threshold are skipped — a cheap "is the room masked?" gate so BirdNET runs where birds are actually legible, not under a drone. Returns per-window detections and an aggregated species tally.

Source code in src/ambiscape/ml.py

def birdnet_session(sess, F=None, windows=None, win_s: float = 9.0,
                    hifi_max_diffuse: float | None = None,
                    lat: float | None = None, lon: float | None = None,
                    min_conf: float = 0.25) -> dict:
    """Run BirdNET across a session, optionally only on hi-fi windows.

    ``windows`` is an explicit list of absolute start seconds; if omitted,
    the session is tiled in ``win_s`` steps. When ``F`` (cached features)
    and ``hifi_max_diffuse`` are given, windows whose median diffuseness
    exceeds the threshold are skipped — a cheap "is the room masked?" gate
    so BirdNET runs where birds are actually legible, not under a drone.
    Returns per-window detections and an aggregated species tally.
    """
    from .io import read_span
    if windows is None:
        windows = []
        for tk in sess.takes:
            t = tk.start + 1.0
            while t + win_s <= tk.end:
                windows.append(t)
                t += win_s
    tally: dict[str, dict] = {}
    per_window = []
    for t0 in windows:
        if F is not None and hifi_max_diffuse is not None:
            i0 = int(np.searchsorted(F["t"], t0))
            i1 = int(np.searchsorted(F["t"], t0 + win_s))
            if i1 > i0 and float(np.median(F["diffuse"][i0:i1])) > \
                    hifi_max_diffuse:
                continue
        x, fs = read_span(sess, t0, win_s)
        dets = birdnet_window(x[:, 0], fs, lat=lat, lon=lon,
                              min_conf=min_conf)
        if dets:
            per_window.append({"t0_s": float(t0), "detections": dets})
            for d in dets:
                e = tally.setdefault(d["species"], {
                    "common_name": d["common_name"], "n": 0, "max_conf": 0.0})
                e["n"] += 1
                e["max_conf"] = max(e["max_conf"], d["confidence"])
    species = sorted(({"species": k, **v} for k, v in tally.items()),
                     key=lambda s: (-s["n"], -s["max_conf"]))
    return {"n_windows_analyzed": len(windows),
            "n_windows_with_birds": len(per_window),
            "n_species": len(species),
            "species": species, "windows": per_window}

speech_fraction(x, fs)

silero-vad speech statistics for one mono signal.

Source code in src/ambiscape/ml.py

def speech_fraction(x: np.ndarray, fs: int) -> dict:
    """silero-vad speech statistics for one mono signal."""
    import torch
    from silero_vad import load_silero_vad, get_speech_timestamps
    model = load_silero_vad()
    y = _resample(x.astype(np.float32), fs, 16000)
    ts = get_speech_timestamps(torch.from_numpy(np.ascontiguousarray(y)),
                               model, sampling_rate=16000)
    dur = len(y) / 16000
    speech = sum((t["end"] - t["start"]) for t in ts) / 16000
    return {"duration_s": round(dur, 1),
            "speech_s": round(speech, 1),
            "speech_fraction": round(speech / dur, 4) if dur else 0.0,
            "n_speech_segments": len(ts),
            "first_speech_at_s": round(ts[0]["start"] / 16000, 1) if ts else None}

speech_gate(path, threshold=0.01)

Privacy gate for a WAV file (any channel count; W/ch0 is analyzed).

Returns the speech statistics plus a pass/fail verdict against threshold (default: fail if more than 1 % of the file is speech).

Source code in src/ambiscape/ml.py

def speech_gate(path: str | Path, threshold: float = 0.01) -> dict:
    """Privacy gate for a WAV file (any channel count; W/ch0 is analyzed).

    Returns the speech statistics plus a pass/fail verdict against
    `threshold` (default: fail if more than 1 % of the file is speech).
    """
    import soundfile as sf
    x, fs = sf.read(str(path), dtype="float32", always_2d=True)
    res = speech_fraction(x[:, 0], fs)
    res["file"] = str(path)
    res["passes"] = res["speech_fraction"] <= threshold
    return res

Deposit export

Non-identifying feature export in the StillStanding365 deposit schema.

Writes one TSV per take with the columns used by the StillStanding365 Zenodo deposit (audio/{day}.tsv): per-second Time, level_dbfs, centroid_hz, low_frac (< 250 Hz), high_frac (> 2 kHz). A 1 Hz loudness/spectral envelope is far below speech timescales and carries no intelligible content, so these files are safe to publish where raw audio is not.

Method notes vs. the original extract_audio.py: levels here come from the W (omni) channel at native rate (the original used an ffmpeg 4-channel downmix at 8 kHz — offsets of a few tenths of a dB are expected), and band fractions are power fractions from the cached log-spectrogram (the original used magnitude fractions of an 8 kHz FFT). Trends and dynamics are directly comparable; absolute fraction values differ slightly by construction.