Audio model with your own MFCC

cd ~/sutraai/data/raw/audio
for f in */*.m4a; do ffmpeg -loglevel error -y -i "$f" -ac 1 -ar 16000 "${f%.m4a}.wav"; done
ls */*.wav | wc -l

import librosa, librosa.display, numpy as np, matplotlib.pyplot as plt
from pathlib import Path
ROOT = Path.home() / "sutraai"
files = sorted((ROOT / "data/raw/audio").glob("*/*.wav"))
by_label = {}
for f in files:
    y, sr = librosa.load(f, sr=16000, mono=True)
    by_label.setdefault(f.parent.name, []).append(len(y) / sr)
print({k: f"{sum(v)/60:.1f} min in {len(v)} files" for k, v in by_label.items()})

fig, axes = plt.subplots(1, len(by_label), figsize=(5 * len(by_label), 3), squeeze=False)
for ax, label in zip(axes[0], by_label):
    f = next(p for p in files if p.parent.name == label)
    y, _ = librosa.load(f, sr=16000, duration=3)
    S = librosa.amplitude_to_db(np.abs(librosa.stft(y, n_fft=512, hop_length=160)), ref=np.max)
    librosa.display.specshow(S, sr=16000, hop_length=160, x_axis="time", y_axis="hz", ax=ax); ax.set_title(label)
plt.show()

import numpy as np

def hz_to_mel(f): return 2595 * np.log10(1 + f / 700)
def mel_to_hz(m): return 700 * (10 ** (m / 2595) - 1)

def mel_filterbank(sr, n_fft, n_mels):
    mels = np.linspace(hz_to_mel(0), hz_to_mel(sr / 2), n_mels + 2)
    bins = np.floor((n_fft + 1) * mel_to_hz(mels) / sr).astype(int)
    fb = np.zeros((n_mels, n_fft // 2 + 1))
    for m in range(1, n_mels + 1):
        l, c, r = bins[m - 1], bins[m], bins[m + 1]
        if c > l: fb[m - 1, l:c] = (np.arange(l, c) - l) / (c - l)
        if r > c: fb[m - 1, c:r] = (r - np.arange(c, r)) / (r - c)
    return fb

def mfcc(x, sr=16000, frame_ms=20, stride_ms=20, n_fft=512, n_mels=32, n_coef=13, preemph=0.97):
    x = np.asarray(x, dtype=np.float64)
    if preemph:
        x = np.append(x[0], x[1:] - preemph * x[:-1])                                  # 1 pre-emphasis
    flen, step = sr * frame_ms // 1000, sr * stride_ms // 1000
    n = 1 + (len(x) - flen) // step
    idx = np.arange(flen)[None, :] + step * np.arange(n)[:, None]
    frames = x[idx] * np.hamming(flen)                                                   # 2 frames + window
    power = np.abs(np.fft.rfft(frames, n_fft)) ** 2 / n_fft                              # 3 power spectrum
    logmel = np.log(power @ mel_filterbank(sr, n_fft, n_mels).T + 1e-10)                 # 4 log mel energies
    k = np.arange(n_mels)
    dct = np.cos(np.pi / n_mels * (k[None, :] + 0.5) * np.arange(n_coef)[:, None])      # 5 DCT-II
    return (logmel @ dct.T).astype(np.float32)                                           # [frames, n_coef]

%load_ext autoreload
%autoreload 2
import sys, numpy as np, librosa, matplotlib.pyplot as plt
from pathlib import Path
ROOT = Path.home() / "sutraai"; sys.path.insert(0, str(ROOT))
from sutra_core.dsp import mfcc

wav = sorted((ROOT / "data/raw/audio/rose_ringed_parakeet").glob("*.wav"))[0]
x, _ = librosa.load(wav, sr=16000, mono=True)
x = x[:16000]                                              # exactly 1 second
M = mfcc(x)
print(M.shape)                                             # expect (50, 13)
plt.imshow(M.T, aspect="auto", origin="lower"); plt.xlabel("frame (20 ms)"); plt.ylabel("coefficient"); plt.colorbar(); plt.show()

ours = mfcc(x, n_fft=320, preemph=0.0)                     # same 320-sample frames as librosa below
ref = librosa.feature.mfcc(y=x, sr=16000, n_mfcc=13, n_fft=320, win_length=320, hop_length=320,
                           n_mels=32, center=False, htk=True, window="hamming").T
n = min(len(ours), len(ref))
corr = [np.corrcoef(ours[:n, k], ref[:n, k])[0, 1] for k in range(13)]
print("per-coefficient correlation:", np.round(corr, 3), "| mean", round(float(np.mean(corr)), 3))

gold = ROOT / "exports/golden"; gold.mkdir(parents=True, exist_ok=True)
rng = np.random.default_rng(0)
for i, start in enumerate(rng.integers(0, max(1, len(x) - 16000), size=3)):
    seg = x[start:start + 16000]
    np.savez(gold / f"mfcc_{i}.npz", x=seg.astype(np.float32), y=mfcc(seg),
             params=np.array([16000, 20, 20, 512, 32, 13], dtype=np.int32))
print(sorted(p.name for p in gold.iterdir()))

%load_ext autoreload
%autoreload 2
import sys, json, random, torch, numpy as np, librosa, torch.nn as nn
from pathlib import Path
ROOT = Path.home() / "sutraai"; sys.path.insert(0, str(ROOT))
from sutra_core.dsp import mfcc
from sutra_core.models import Normalize, Conv1DNet
from sutra_core.train import fit_arrays, accuracy

SR, WIN, HOP = 16000, 16000, 8000                       # 1 s windows, 0.5 s stride
files = sorted((ROOT / "data/raw/audio").glob("*/*.wav"))
labels = sorted({f.parent.name for f in files})
noise_clips = [librosa.load(f, sr=SR)[0] for f in files if f.parent.name == "noise"]

def mix_noise(clip, noise, snr_db):
    noise = np.resize(noise, clip.shape)
    p_c, p_n = np.mean(clip ** 2) + 1e-12, np.mean(noise ** 2) + 1e-12
    return clip + noise * np.sqrt(p_c / (p_n * 10 ** (snr_db / 10)))

random.seed(42)
split_of = {f: random.choices(["train", "val", "test"], weights=[70, 15, 15])[0] for f in files}   # split by FILE

def build(split, augment=False):
    X, y = [], []
    for f in files:
        if split_of[f] != split: continue
        sig, _ = librosa.load(f, sr=SR, mono=True)
        for s in range(0, len(sig) - WIN + 1, HOP):
            w = sig[s:s + WIN]
            X.append(mfcc(w)); y.append(labels.index(f.parent.name))
            if augment and f.parent.name != "noise" and noise_clips:
                X.append(mfcc(mix_noise(w, random.choice(noise_clips), random.uniform(5, 20))))
                y.append(labels.index(f.parent.name))
    return np.stack(X).astype(np.float32), np.array(y)

Xtr, ytr = build("train", augment=True); Xva, yva = build("val"); Xte, yte = build("test")
mu, sd = Xtr.mean(axis=(0, 1)), Xtr.std(axis=(0, 1)) + 1e-6
print(Xtr.shape, Xva.shape, Xte.shape, labels)

torch.manual_seed(42)
model = nn.Sequential(Normalize(mu, sd), Conv1DNet(13, len(labels)))
model, val_acc = fit_arrays(model, Xtr, ytr, Xva, yva, epochs=100, lr=5e-3)
test_acc = accuracy(model, Xte, yte)
print("audio test accuracy:", round(test_acc, 3))

from sklearn.metrics import confusion_matrix
EXP = ROOT / "exports"
model.eval().cpu()
torch.onnx.export(model, torch.randn(1, 50, 13), str(EXP / "audio_v1.onnx"),
                  input_names=["input"], output_names=["logits"], opset_version=17, dynamo=False)
np.save(EXP / "calib_audio.npy", Xtr[np.random.default_rng(0).choice(len(Xtr), 200, replace=False)])
from sutra_core.export import onnx2tf_int8                   # helper from Step 05 (PATH fix + per-tensor int8)
tfl = onnx2tf_int8(EXP / "audio_v1.onnx", EXP / "audio_v1_tf", EXP / "calib_audio.npy", "[[[0.0]]]", "[[[1.0]]]")
print(tfl.name, round(tfl.stat().st_size / 1024, 1), "KB")

with torch.no_grad(): pred = model(torch.tensor(Xte)).argmax(1).numpy()
(ROOT / "reports/audio_v1").mkdir(parents=True, exist_ok=True)
json.dump({"run": "audio_v1", "classes": labels, "val_acc": val_acc, "test_acc": test_acc,
           "confusion_matrix": confusion_matrix(yte, pred).tolist()},
          open(ROOT / "reports/audio_v1/metrics.json", "w"), indent=2)