aural.py

#!/usr/bin/env python3
"""
Aural micro-Doppler playback.

Converts a micro-Doppler spectrogram (or raw IQ from blah2) into a WAV file
that captures the rotor signature audibly. The human auditory system is
"particularly robust in the presence of noise" (Chen §9.3.4), so playing
the rotor harmonics over speakers gives judges a memorable second channel
on the demo.

Two modes:

  npz     Take a saved (frequency, time, magnitude) numpy array (from blah2's
          IQ-save mode + scipy.signal.spectrogram) and resynthesize audio.

  iq      Take a raw IQ recording (Pluto .pluto / .iq complex64) and run
          STFT + downshift the Doppler band into the audible range, write WAV.

Usage:
    python3 aural.py iq capture.iq -o rotor.wav --duration 6
    python3 aural.py npz spectrogram.npz -o rotor.wav --shift 200

Requires: numpy, scipy.
"""
import argparse
import sys
from pathlib import Path

import numpy as np
from scipy import signal
from scipy.io import wavfile


def iq_to_audio(iq_path: Path, *, sr_in: float, audio_sr: int, duration: float,
                doppler_band_hz: tuple[float, float] = (-3000, 3000),
                shift_hz: float = 800.0) -> tuple[int, np.ndarray]:
    """Read complex IQ, baseband-filter to the rotor Doppler band, frequency-shift
    into the audible range, and emit mono int16 audio.

    Args:
        iq_path: file with int16 interleaved I/Q (Pluto raw) OR complex64.
        sr_in:   IQ sample rate (Hz).
        audio_sr: output WAV sample rate (44_100 typical).
        duration: max output seconds.
        doppler_band_hz: (low, high) micro-Doppler band of interest.
        shift_hz: how far to lift the band into audio so 0 Hz ≠ silence.
    """
    raw = np.fromfile(iq_path, dtype=np.int16)
    if raw.size == 0:
        raise SystemExit(f"{iq_path}: empty file")
    if raw.size % 2 != 0:
        raw = raw[:-1]
    iq = raw[0::2].astype(np.float32) + 1j * raw[1::2].astype(np.float32)
    iq /= 32768.0

    n_keep = int(min(duration, len(iq) / sr_in) * sr_in)
    iq = iq[:n_keep]

    # 1) low-pass to the Doppler band of interest (centered at zero in baseband)
    bw_lo, bw_hi = doppler_band_hz
    nyq = sr_in / 2
    sos = signal.butter(8, [max(bw_lo, -nyq + 1) / nyq, min(bw_hi, nyq - 1) / nyq],
                        btype="bandpass", output="sos")
    # bandpass on complex IQ: filter I and Q separately
    iq_f = signal.sosfiltfilt(sos, iq.real) + 1j * signal.sosfiltfilt(sos, iq.imag)

    # 2) decimate IQ to audio_sr * 2 (oversample for shift step)
    decim = max(1, int(sr_in // (audio_sr * 4)))
    iq_d = signal.decimate(iq_f, decim, ftype="fir", zero_phase=True)
    sr_d = sr_in / decim

    # 3) shift the rotor band to audible range — multiply by exp(j*2pi*f0*t)
    t = np.arange(len(iq_d)) / sr_d
    iq_shift = iq_d * np.exp(1j * 2 * np.pi * shift_hz * t)

    # 4) take real part as the audio signal, decimate to audio_sr
    audio = iq_shift.real
    decim2 = max(1, int(sr_d // audio_sr))
    audio = signal.decimate(audio, decim2, ftype="fir", zero_phase=True)
    sr_out = int(round(sr_d / decim2))

    # 5) resample to exact audio_sr
    if sr_out != audio_sr:
        n_out = int(round(len(audio) * audio_sr / sr_out))
        audio = signal.resample(audio, n_out)
        sr_out = audio_sr

    # 6) normalize to int16
    peak = np.max(np.abs(audio)) + 1e-9
    audio = (audio / peak * 0.95 * 32767).astype(np.int16)
    return sr_out, audio


def npz_to_audio(npz_path: Path, audio_sr: int, shift_hz: float, duration: float):
    """Resynthesize from a (f, t, |S|) spectrogram. Each time slice's magnitudes
    drive sinusoid amplitudes at frequency = f + shift_hz."""
    data = np.load(npz_path)
    f = data["f"]; t = data["t"]; S = np.abs(data["Sxx"])
    if duration > 0:
        keep = t <= duration
        t = t[keep]; S = S[:, keep]
    n_out = int(t[-1] * audio_sr)
    samples = np.linspace(0, t[-1], n_out, endpoint=False)
    # Linear interpolation of each frequency bin's envelope, then sum.
    audio = np.zeros(n_out, dtype=np.float32)
    for i, fi in enumerate(f):
        env = np.interp(samples, t, S[i])
        audio += env * np.sin(2 * np.pi * (fi + shift_hz) * samples)
    peak = np.max(np.abs(audio)) + 1e-9
    audio = (audio / peak * 0.95 * 32767).astype(np.int16)
    return audio_sr, audio


def main():
    p = argparse.ArgumentParser()
    sub = p.add_subparsers(dest="mode", required=True)

    iq = sub.add_parser("iq", help="convert raw IQ recording to rotor-signature WAV")
    iq.add_argument("input", type=Path)
    iq.add_argument("-o", "--output", type=Path, default=Path("rotor.wav"))
    iq.add_argument("--sr-in", type=float, default=20_000_000)
    iq.add_argument("--audio-sr", type=int, default=44_100)
    iq.add_argument("--duration", type=float, default=6.0)
    iq.add_argument("--shift", type=float, default=800.0,
                    help="Hz to shift the rotor band up into audio")

    npz = sub.add_parser("npz", help="resynthesize from saved spectrogram (f,t,Sxx)")
    npz.add_argument("input", type=Path)
    npz.add_argument("-o", "--output", type=Path, default=Path("rotor.wav"))
    npz.add_argument("--audio-sr", type=int, default=44_100)
    npz.add_argument("--duration", type=float, default=6.0)
    npz.add_argument("--shift", type=float, default=400.0)

    args = p.parse_args()

    if args.mode == "iq":
        sr, audio = iq_to_audio(args.input, sr_in=args.sr_in, audio_sr=args.audio_sr,
                                duration=args.duration, shift_hz=args.shift)
    else:
        sr, audio = npz_to_audio(args.input, args.audio_sr, args.shift, args.duration)

    wavfile.write(str(args.output), sr, audio)
    print(f"# wrote {args.output} ({len(audio)/sr:.2f}s @ {sr}Hz)", file=sys.stderr)


if __name__ == "__main__":
    main()