#!/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')
        e = (1, 2)[ch-1:]
        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 db_spec_mono(x, nfft=16384, hop=2048, sr=SR):
    # averaged magnitude spectrum in dB
    w = np.hanning(nfft)
    frames = [x[t:t+nfft] for t in range(0, len(x)-nfft, hop)]
    if not frames: return None
    X = np.stack([np.fft.rfft(f*w) for f in frames])
    return 20*np.log10(np.abs(X).mean(0) + 1e-12), np.fft.rfftfreq(nfft, 1/sr)

def main():
    inp = sys.argv[1]
    outs = sys.argv[2:]
    x, ch = read_wav(inp)
    bin, fr = db_spec_mono(x)
    NFFT = 16384
    # collect peaks (sines) of test tone from input
    masks = {}
    peakset = set()
    for f in range(1, len(fr)):
        if bin[f] > -120 and bin[f] == np.max(bin[max(0,f-8):f+9]):
            peakset.add(f)
    print("peaks in input (Hz):", sorted(round(fr[f]) for f in peakset if fr[f] < 18000))
    colw = 22
    hdr = f"{'Hz':>6} {'IN dB':>7}" + ''.join(f" {os.path.basename(p)[:colw]:>{colw}}" for p in outs)
    print(hdr)
    for f in sorted(peakset):
        if fr[f] < 18000:
            row = f"{fr[f]:6.0f} {bin[f]:7.1f}"
            refbase = None
            for p in outs:
                y, _ = read_wav(p)
                b, f2 = db_spec_mono(y)
                # nearest bin to fr[f]
                idx = int(round(fr[f]*(len(f2)-1)/f2[-1]))
                idx = int(np.argmin(np.abs(f2 - fr[f])))
                if p == outs[0]:
                    refbase = b[idx]
                row += f" {b[idx]:{colw}.2f}"
            print(row)

if __name__ == '__main__':
    main()