#!/usr/bin/env python3
import struct, os, sys, math

SR = 44100

def read24(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]
        if cid == b'data':
            data_off = i + 8; data_sz = sz
        i += 8 + sz + (sz & 1)
    raw = d[data_off:data_off+data_sz]
    L = []
    for k in range(data_sz // (ch*3)):
        v = int.from_bytes(raw[k*ch*3:k*ch*3+3], 'little', signed=True)
        L.append(v / 8388608.0)
    return L, ch

def spectral_profile(samples):
    # Welch-ish: split into 128-sample frames, hann window, count per-bin RMS
    n = len(samples)
    win = [0.5 - 0.5*math.cos(2*math.pi*k/(128)) for k in range(128)]
    bins = [0.0]*65
    frames = 0
    for s in range(0, n-128, 64):
        fr = [samples[s+k]*win[k] for k in range(128)]
        re = [0.0]*65; im = [0.0]*65
        for k in range(128):
            f_c = fr[k]
            for b in range(65):
                ang = 2*math.pi*b*k/128
                re[b] += f_c*math.cos(ang)
                im[b] += f_c*math.sin(ang)
        for b in range(65):
            bins[b] += re[b]*re[b] + im[b]*im[b]
        frames += 1
    return [math.sqrt(b/frames) for b in bins]

def profile_db(prof):
    return [20*math.log10(p+1e-12) for p in prof]

def main(inp, outs):
    Li, ch = read24(inp)
    pi = spectral_profile(Li)
    pids = {}
    for _, path in outs:
        L, _ = read24(path)
        pids[path] = spectral_profile(L)
    # clamp frames for profile in case lengths differ
    return pi, pids

if __name__ == '__main__':
    inp = sys.argv[1]
    outs = [(os.path.basename(p).split('_')[1] if False else p, p) for p in sys.argv[2:]]
    pi, pids = main(inp, outs)
    freqs = [SR*k/256 for k in range(65)]
    print(f"{'freq':>7} {'IN':>8}" + ''.join(f" {os.path.basename(p):>18}" for _,p in outs))
    print("LOG2-ish band graph")
    # print compact dB difference
    dbin = profile_db(pi)
    for b in range(2, 64, 2):
        row = f"{freqs[b]:7.0f} {dbin[b]:8.1f}"
        for _, p in outs:
            pp = pids[p]
            d = 20*math.log10(pp[b]+1e-12) - dbin[b]
            row += f" {d:18.2f}"
        print(row)