#!/usr/bin/env python3
import sys, os, struct
import numpy as np

SR = 44100

def read_wav(path):
    d = open(path, 'rb').read()
    i = 12
    ch = bps = None
    data_off = data_sz = None
    while i + 8 <= len(d):
        cid = d[i:i+4]; sz = struct.unpack('<I', d[i+4:i+8])[0]
        if cid == b'fmt ':
            ch = struct.unpack('<H', d[i+8:i+10])[0]
            bps = struct.unpack('<H', d[i+16:i+18])[0]
        if cid == b'data':
            data_off = i + 8; data_sz = sz
        i += 8 + sz + (sz & 1)
    raw = d[data_off:data_off+data_sz]
    if bps == 32:
        a = np.frombuffer(raw, dtype='<f4')
        return a.reshape(-1, ch)[:, 0], ch
    if bps == 24:
        b = np.frombuffer(raw, dtype=np.uint8).reshape(-1, ch*3)
        v = b[:, 0] | (b[:,1].astype(np.int32) << 8) | (b[:,2].astype(np.int32) << 16)
        v = (v ^ (1 << 23)) - (1 << 23)
        return v / 8388608.0, ch
    a = np.frombuffer(raw, dtype='<i2')
    return a.reshape(-1, ch)[:, 0] / 32768.0, ch

def load(path):
    return read_wav(path)[0]

def band_db(x, fl, fh, sr=SR):
    # band-integrated power between fl..fh via Goertzel on windowed segments
    n = len(x)
    seg = 4096
    hop = 2048
    total = 0.0
    frames = 0
    for t in range(0, n-seg, hop):
        fr = x[t:t+seg]
        w = np.hanning(seg)
        X = np.fft.rfft(fr*w)
        frq = np.fft.rfftfreq(seg, 1/sr)
        m = (frq >= fl) & (frq <= fh)
        total += float(np.sum(np.abs(X[m])**2))
        frames += 1
    if frames == 0: return -300.0
    return 10*np.log10(total/frames + 1e-30)

def main():
    # probe band around the carrier at flexible freqs
    center = float(sys.argv[1])
    outs = sys.argv[2:]
    half_w = float(os.environ.get('PROBE_BW', 800.0))
    widths = []
    step = float(os.environ.get('PROBE_STEP', 25.0))
    f = center - half_w
    while f <= center + half_w:
        widths.append(f)
        f += step
    datas = {p: load(p) for p in outs}
    colw = 11
    print(f"{'freq':>6}" + ''.join(f" {os.path.basename(p)[:colw]:>{colw}}" for p in outs))
    prevbase = None
    for fl in widths:
        fh = fl + step
        row = f"{fl + step/2:6.0f}"
        base = None
        for p in outs:
            d = band_db(datas[p], fl, fh)
            if base is None: base = d
            row += f" {d:>{colw}.1f}"
        print(row)
    # also print relative-to-first file per bin is tricky; print a second table relative to first col
    print("--- relative dB to first file ---")
    for fl in widths:
        fh = fl + step
        row = f"{fl + step/2:6.0f}"
        first = None
        for p in outs:
            d = band_db(datas[p], fl, fh)
            if first is None: first = d
            row += f" {d - first:>+{colw}.1f}"
        print(row)

if __name__ == '__main__':
    main()