Audio: del sonido a los números

import re
import tarfile
import textwrap
import urllib.request
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import soundfile as sf
import torch
import torchaudio
from torch import nn

torch.manual_seed(0)
device = "cuda" if torch.cuda.is_available() else "cpu"
DATOS = Path("datos")


def bajar(url, destino):
    """Baja un archivo de internet, solo si todavía no está."""
    destino = Path(destino)
    if not destino.exists():
        destino.parent.mkdir(parents=True, exist_ok=True)
        pedido = urllib.request.Request(
            url, headers={"User-Agent": "ruta-obia/1.0 (curso de IA)"})
        with urllib.request.urlopen(pedido) as respuesta:
            destino.write_bytes(respuesta.read())
    return destino


# Don Quijote, cap. I, leído por Nessi para LibriVox (dominio público)
wav = DATOS / "audio-ruta/quijote-cap1-55s-16k.wav"
if not wav.exists():
    url = ("https://archive.org/download/don_quijote_vol1_0706_librivox/"
           "quijote_vol1_03_cervantes_64kb.mp3")
    mp3 = bajar(url, DATOS / "audio-ruta/quijote-vol1-03-librivox.mp3")
    # Del segundo 36,0 al 90,8: antes viene el aviso de LibriVox
    trozo, sr_mp3 = sf.read(mp3, dtype="float32",
                            start=int(36.0 * 22050),
                            stop=int(90.8 * 22050))
    trozo = torchaudio.functional.resample(torch.from_numpy(trozo),
                                           sr_mp3, 16000)
    sf.write(wav, trozo.numpy(), 16000)

onda, sr = sf.read(wav, dtype="float32")
print("Forma:", onda.shape, "  tipo:", onda.dtype)
print("Frecuencia de muestreo:", sr, "muestras por segundo")
print(f"Duración: {len(onda) / sr:.1f} s")
print("Seis muestras seguidas:", np.round(onda[16000:16006], 3))

t = np.arange(len(onda)) / sr              # el tiempo de cada muestra
tres = t < 3
zoom = (t >= 1.30) & (t < 1.32)            # 20 milisegundos

fig, (arriba, abajo) = plt.subplots(2, 1, figsize=(4.8, 4.8))
arriba.plot(t[tres], onda[tres], color="C0", linewidth=0.5)
arriba.axvspan(1.30, 1.32, color="C1")
arriba.set_xlabel("Tiempo (s)")
arriba.set_ylabel("Amplitud")
arriba.set_title("Los primeros 3 segundos")
abajo.plot(t[zoom] * 1000, onda[zoom], color="C2", linewidth=1.2)
abajo.scatter(t[zoom] * 1000, onda[zoom], s=12, color="C0",
              edgecolors="none")
abajo.set_xlabel("Tiempo (ms)")
abajo.set_ylabel("Amplitud")
abajo.set_title("Zoom: 20 ms, una muestra por punto")
plt.show()

f_tono = 900
t_fino = np.linspace(0, 0.01, 2000)        # «continuo»: 10 ms
t_lento = np.arange(0, 0.01, 1 / 1000)     # 1000 muestras por segundo

plt.plot(t_fino * 1000, np.cos(2 * np.pi * f_tono * t_fino),
         color="C2", linewidth=1, label="Tono real de 900 Hz")
plt.plot(t_fino * 1000, np.cos(2 * np.pi * 100 * t_fino),
         color="C0", linestyle="--", label="Lo que parece: 100 Hz")
plt.scatter(t_lento * 1000, np.cos(2 * np.pi * f_tono * t_lento),
            color="C1", s=60, zorder=3, label="Muestras a 1000 Hz")
plt.xlabel("Tiempo (ms)")
plt.ylabel("Amplitud")
plt.title("Muestrear muy lento engaña")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.3), fontsize=13)
plt.show()

# ¿Las muestras de 900 Hz y de 100 Hz son las mismas?
n = np.arange(1000)
print(np.allclose(np.cos(2 * np.pi * 900 * n / 1000),
                  np.cos(2 * np.pi * 100 * n / 1000)))

sr_m = 8000
t1 = np.arange(sr_m) / sr_m                # un segundo
rng = np.random.default_rng(0)
mezcla = (1.0 * np.sin(2 * np.pi * 220 * t1)
          + 0.5 * np.sin(2 * np.pi * 440 * t1)
          + 0.25 * np.sin(2 * np.pi * 660 * t1)
          + 0.05 * rng.normal(size=sr_m))

espectro = np.fft.rfft(mezcla)             # números complejos
frecuencias = np.fft.rfftfreq(len(mezcla), d=1 / sr_m)
amplitud = 2 * np.abs(espectro) / len(mezcla)
print("Muestras:", len(mezcla), "  frecuencias:", len(espectro))
print(f"De {frecuencias[0]:.0f} a {frecuencias[-1]:.0f} Hz, "
      f"cada {frecuencias[1]:.0f} Hz")
for i in np.argsort(amplitud)[::-1][:3]:
    print(f"{frecuencias[i]:5.0f} Hz: amplitud {amplitud[i]:.2f}")

plt.plot(frecuencias, amplitud, color="C0", linewidth=1.5)
plt.scatter([220], [amplitud[220]], color="C1", s=60, zorder=3)
plt.xlim(0, 1000)
plt.xlabel("Frecuencia (Hz)")
plt.ylabel("Amplitud")
plt.title("La FFT encuentra los tres tonos")
plt.show()

def stft(senal, n_fft, salto):
    """Magnitud de la STFT: una columna por ventana."""
    hann = np.hanning(n_fft)
    inicios = range(0, len(senal) - n_fft + 1, salto)
    cuadros = np.array([senal[i:i + n_fft] * hann for i in inicios])
    return np.abs(np.fft.rfft(cuadros, axis=1)).T   # (frec., tiempo)


def en_db(magnitud):
    """Decibeles: 20·log10 de la amplitud, con un piso para no tomar
    log de 0."""
    return 20 * np.log10(np.maximum(magnitud, 1e-4))

from matplotlib.colors import LinearSegmentedColormap

tinta = LinearSegmentedColormap.from_list("tinta", ["white", "C0"])

senal = np.sin(2 * np.pi * 500 * t1) + np.sin(2 * np.pi * 560 * t1)
for inicio in [0.25, 0.5, 0.75]:
    i = int(inicio * sr_m)
    senal[i:i + 8] += 6                    # un chasquido de 1 ms

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
ejes[0, 0].plot(t1, senal, color="C0", linewidth=0.4)
ejes[0, 0].set_title("La señal")
ejes[0, 0].set_xlabel("Tiempo (s)")
for eje, n_fft in zip(ejes.flat[1:], [64, 256, 2048]):
    S = en_db(stft(senal, n_fft, n_fft // 4))
    eje.imshow(S, origin="lower", aspect="auto", cmap=tinta,
               vmin=S.max() - 40, extent=[0, 1, 0, sr_m / 2])
    eje.set_ylim(0, 1500)
    eje.set_yticks([500, 1000, 1500])
    eje.set_title(f"n_fft = {n_fft}")
    eje.set_xlabel("Tiempo (s)")
    eje.grid(False)
ejes[0, 0].set_yticks([])
ejes[1, 1].set_yticklabels([])
ejes[1, 0].set_ylabel("Frecuencia (Hz)")
plt.show()

for n_fft in [64, 256, 2048]:
    print(f"n_fft = {n_fft:4}: ventana de {1000 * n_fft / sr_m:5.1f} ms,"
          f" una fila cada {sr_m / n_fft:6.2f} Hz")

for f in [250, 500, 1000, 2000, 4000, 8000]:
    print(f"{f:5} Hz -> {2595 * np.log10(1 + f / 700):6.0f} mels")

a_mel = torchaudio.transforms.MelSpectrogram(
    sample_rate=sr, n_fft=400, hop_length=160, n_mels=64)
a_db = torchaudio.transforms.AmplitudeToDB()
voz = onda[:3 * sr]                        # los primeros 3 segundos
S_lineal = en_db(stft(voz, 400, 160))
S_mel = a_db(a_mel(torch.from_numpy(voz))).numpy()
print("STFT:", S_lineal.shape, "  Mel:", S_mel.shape)

fig, (arriba, abajo) = plt.subplots(2, 1, figsize=(4.8, 4.8))
arriba.imshow(S_lineal, origin="lower", aspect="auto", cmap=tinta,
              vmin=S_lineal.max() - 70, extent=[0, 3, 0, sr / 2])
arriba.set_title("Espectrograma (Hz)")
arriba.set_ylabel("Frecuencia (Hz)")
abajo.imshow(S_mel, origin="lower", aspect="auto", cmap=tinta,
             vmin=S_mel.max() - 70, extent=[0, 3, 0, 64])
abajo.set_title("Espectrograma de Mel (64 bandas)")
abajo.set_ylabel("Banda de Mel")
abajo.set_xlabel("Tiempo (s)")
for eje in (arriba, abajo):
    eje.grid(False)
plt.show()

# Free Spoken Digit Dataset, de Zohar Jackson y colaboradores,
# con licencia CC BY-SA 4.0 (github.com/Jakobovski/free-spoken-digit-dataset)
carpeta = DATOS / "fsdd/recordings"
if not carpeta.exists():
    url = ("https://github.com/Jakobovski/free-spoken-digit-dataset/"
           "archive/refs/tags/v1.0.10.tar.gz")
    comprimido = bajar(url, DATOS / "fsdd/fsdd-v1.0.10.tar.gz")
    with tarfile.open(comprimido) as tar:
        for miembro in tar.getmembers():
            if "/recordings/" in miembro.name:
                miembro.name = "recordings/" + Path(miembro.name).name
                tar.extract(miembro, DATOS / "fsdd", filter="data")

audios, digitos, hablantes, tomas = [], [], [], []
for archivo in sorted(carpeta.glob("*.wav")):
    # El nombre dice todo: «7_jackson_32.wav»
    digito, hablante, toma = archivo.stem.split("_")
    x, sr_d = sf.read(archivo, dtype="float32")
    audios.append(x)
    digitos.append(int(digito))
    hablantes.append(hablante)
    tomas.append(int(toma))
digitos, tomas = np.array(digitos), np.array(tomas)

duraciones = np.array([len(x) for x in audios]) / sr_d
print("Audios:", len(audios), "  muestreo:", sr_d, "Hz")
print("Hablantes:", sorted(set(hablantes)))
print(f"Duración: mínima {duraciones.min():.2f} s, "
      f"mediana {np.median(duraciones):.2f} s, "
      f"máxima {duraciones.max():.2f} s")
print("Más largos que 1 s:", (duraciones > 1).sum())

LARGO = sr_d                               # 1 segundo = 8000 muestras


def a_un_segundo(x):
    """Recorta o rellena con silencio hasta LARGO muestras."""
    if len(x) >= LARGO:
        return x[:LARGO]
    return np.pad(x, (0, LARGO - len(x)))


ondas = torch.from_numpy(np.stack([a_un_segundo(x) for x in audios]))
mel_d = torchaudio.transforms.MelSpectrogram(
    sample_rate=sr_d, n_fft=256, hop_length=80, n_mels=40)
log_mel = torch.log(mel_d(ondas) + 1e-6)   # (audios, bandas, cuadros)
print("Ondas:", tuple(ondas.shape), "  log-Mel:", tuple(log_mel.shape))

# Los diez dígitos dichos por la misma persona
fig, ejes = plt.subplots(5, 2, figsize=(4.8, 4.8))
for d, eje in enumerate(ejes.T.flat):
    k = np.flatnonzero((digitos == d) & (np.array(hablantes) == "jackson")
                       & (tomas == 0))[0]
    eje.imshow(log_mel[k], origin="lower", aspect="auto", cmap=tinta)
    eje.text(0.97, 0.9, str(d), transform=eje.transAxes, ha="right",
             va="top", fontsize=15, fontweight="bold", color="C0")
    eje.set_xticks([])
    eje.set_yticks([])
fig.suptitle("Log-Mel de los dígitos 0 a 9")
plt.show()

from torch.utils.data import DataLoader, TensorDataset

es_prueba = torch.from_numpy(tomas < 5)
es_val = torch.from_numpy((tomas >= 5) & (tomas < 10))
es_ent = torch.from_numpy(tomas >= 10)
X = log_mel.unsqueeze(1)                   # un canal: (N, 1, 40, 101)
y = torch.from_numpy(digitos)
media, desvio = X[es_ent].mean(), X[es_ent].std()
X = (X - media) / desvio

cargador_ent = DataLoader(
    TensorDataset(X[es_ent], y[es_ent]), batch_size=64,
    shuffle=True, generator=torch.Generator().manual_seed(0))
cargador_val = DataLoader(TensorDataset(X[es_val], y[es_val]),
                          batch_size=300)
cargador_pru = DataLoader(TensorDataset(X[es_prueba], y[es_prueba]),
                          batch_size=300)
print("Entrenamiento:", len(cargador_ent.dataset),
      "  validación:", len(cargador_val.dataset),
      "  prueba:", len(cargador_pru.dataset))

def bloque(entran, salen):
    """Convolución 3×3 + BatchNorm + ReLU."""
    return nn.Sequential(nn.Conv2d(entran, salen, 3, padding=1),
                         nn.BatchNorm2d(salen), nn.ReLU())


torch.manual_seed(0)
cnn = nn.Sequential(
    bloque(1, 16), nn.MaxPool2d(2),        # (16, 20, 50)
    bloque(16, 32), nn.MaxPool2d(2),       # (32, 10, 25)
    bloque(32, 32),                        # (32, 10, 25)
    nn.AdaptiveAvgPool2d((5, 6)),          # (32, 5, 6)
    nn.Flatten(),                          # 960
    nn.Linear(32 * 5 * 6, 10),
)
print("Parámetros:", sum(p.numel() for p in cnn.parameters()))

def evaluar(modelo, cargador, funcion_perdida, dispositivo="cpu"):
    """Pérdida media y exactitud de `modelo` sobre `cargador`."""
    modelo.eval()
    suma, aciertos, total = 0.0, 0, 0
    with torch.no_grad():
        for x, y in cargador:
            x, y = x.to(dispositivo), y.to(dispositivo)
            salida = modelo(x)
            suma += funcion_perdida(salida, y).item() * len(y)
            aciertos += (salida.argmax(dim=1) == y).sum().item()
            total += len(y)
    return suma / total, aciertos / total


def entrenar(modelo, cargador_ent, cargador_val, epocas, optimizador,
             funcion_perdida=None, planificador=None, dispositivo="cpu",
             mostrar=True):
    """Entrena `modelo` y devuelve el historial de cada época."""
    if funcion_perdida is None:
        funcion_perdida = nn.CrossEntropyLoss()
    modelo.to(dispositivo)
    historial = {"perdida_ent": [], "exactitud_ent": [],
                 "perdida_val": [], "exactitud_val": []}
    for epoca in range(1, epocas + 1):
        modelo.train()
        suma, aciertos, total = 0.0, 0, 0
        for x, y in cargador_ent:
            x, y = x.to(dispositivo), y.to(dispositivo)
            optimizador.zero_grad()
            salida = modelo(x)
            perdida = funcion_perdida(salida, y)
            perdida.backward()
            optimizador.step()
            suma += perdida.item() * len(y)
            aciertos += (salida.argmax(dim=1) == y).sum().item()
            total += len(y)
        if planificador is not None:
            planificador.step()
        perdida_val, exactitud_val = evaluar(
            modelo, cargador_val, funcion_perdida, dispositivo)
        historial["perdida_ent"].append(suma / total)
        historial["exactitud_ent"].append(aciertos / total)
        historial["perdida_val"].append(perdida_val)
        historial["exactitud_val"].append(exactitud_val)
        if mostrar:
            print(f"Época {epoca}: pérdida {suma / total:.3f}, "
                  f"validación {perdida_val:.3f}, "
                  f"exactitud de validación {exactitud_val:.3f}")
    return historial

torch.manual_seed(0)
hist = entrenar(cnn, cargador_ent, cargador_val, epocas=8,
                optimizador=torch.optim.Adam(cnn.parameters(), lr=1e-3),
                dispositivo=device)
# La prueba, una sola vez, con la red ya entrenada
_, exactitud_prueba = evaluar(cnn, cargador_pru, nn.CrossEntropyLoss(),
                              device)
print(f"Exactitud de prueba: {exactitud_prueba:.3f}")

from transformers import AutoFeatureExtractor, AutoModel

extractor = AutoFeatureExtractor.from_pretrained(
    "facebook/hubert-base-ls960")
hubert = AutoModel.from_pretrained("facebook/hubert-base-ls960")
hubert = hubert.to(device).eval()

indices_prueba = np.flatnonzero(tomas < 5)
embeddings = []
with torch.no_grad():
    for k in indices_prueba:
        x16 = torchaudio.functional.resample(
            torch.from_numpy(audios[k]), sr_d, 16000)
        entrada = extractor(x16.numpy(), sampling_rate=16000,
                            return_tensors="pt").to(device)
        cuadros = hubert(**entrada).last_hidden_state  # (1, T, 768)
        embeddings.append(cuadros.mean(dim=1)[0].cpu())
E = torch.stack(embeddings).numpy()
y_h = digitos[indices_prueba]

print(f"Último audio: {len(x16) / 16000:.2f} s ->",
      tuple(cuadros.shape))
print("Embeddings:", E.shape)

import umap
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

plano = umap.UMAP(n_neighbors=15, min_dist=0.1, metric="cosine",
                  random_state=0, n_jobs=1).fit_transform(E)

plt.scatter(plano[:, 0], plano[:, 1], s=0)  # solo fija los límites
for (a, b), d in zip(plano, y_h):
    plt.text(a, b, str(d), ha="center", va="center", fontsize=10,
             color="C1" if d == 7 else "C0", alpha=0.8)
plt.title("HuBERT: 300 audios, cada uno con su dígito")
plt.xticks([])
plt.yticks([])
plt.show()

# ¿Alcanza el embedding para clasificar? Validación cruzada, 5 pliegues
clasificador = LogisticRegression(max_iter=2000)
exactitudes = cross_val_score(clasificador, E, y_h, cv=5)
print("Exactitud por pliegue:", np.round(exactitudes, 2))
print(f"Promedio: {exactitudes.mean():.3f}")

from transformers import AutoModelForSpeechSeq2Seq, AutoProcessor

whisper = {}
for nombre in ["tiny", "base"]:
    ruta = f"openai/whisper-{nombre}"
    procesador = AutoProcessor.from_pretrained(ruta)
    modelo = AutoModelForSpeechSeq2Seq.from_pretrained(ruta)
    whisper[nombre] = (procesador, modelo.to(device).eval())
    millones = sum(p.numel() for p in modelo.parameters()) / 1e6
    print(f"whisper-{nombre}: {millones:.0f} millones de parámetros")


@torch.no_grad()
def transcribir(audio, nombre, idioma):
    """Texto de un audio a 16 000 Hz, con el Whisper `nombre`."""
    procesador, modelo = whisper[nombre]
    largo = len(audio) > 30 * 16000
    if largo:
        # Más de 30 s: Whisper lo procesa en tramos de 30 s
        x = procesador(audio, sampling_rate=16000, return_tensors="pt",
                       truncation=False, padding="longest",
                       return_attention_mask=True)
    else:
        x = procesador(audio, sampling_rate=16000, return_tensors="pt")
    ids = modelo.generate(**x.to(device), language=idioma,
                          task="transcribe", return_timestamps=largo)
    return procesador.batch_decode(ids, skip_special_tokens=True)[0]


def mostrar(texto):
    """Imprime el texto en renglones de 72 caracteres."""
    print(textwrap.fill(texto.strip(), 72))


texto_tiny = transcribir(onda, "tiny", "es")
mostrar(texto_tiny)

import jiwer

REFERENCIA = (
    "En un lugar de la Mancha, de cuyo nombre no quiero acordarme, no ha "
    "mucho tiempo que vivía un hidalgo de los de lanza en astillero, "
    "adarga antigua, rocín flaco y galgo corredor. Una olla de algo más "
    "vaca que carnero, salpicón las más noches, duelos y quebrantos los "
    "sábados, lantejas los viernes, algún palomino de añadidura los "
    "domingos, consumían las tres partes de su hacienda. El resto della "
    "concluían sayo de velarte, calzas de velludo para las fiestas, con "
    "sus pantuflos de lo mesmo, y los días de entresemana se honraba con "
    "su vellorí de lo más fino. Tenía en su casa una ama que pasaba de "
    "los cuarenta, y una sobrina que no llegaba a los veinte, y un mozo "
    "de campo y plaza, que así ensillaba el rocín como tomaba la "
    "podadera.")


def normalizar(texto):
    """Minúsculas, sin puntuación y con un espacio entre palabras."""
    texto = re.sub(r"[^\w\s]", " ", texto.lower())
    return " ".join(texto.split())


ref = normalizar(REFERENCIA)
hip = normalizar(texto_tiny)
cuentas = jiwer.process_words(ref, hip)
print("Palabras de la referencia:", len(ref.split()))
print("Sustituidas:", cuentas.substitutions,
      "  borradas:", cuentas.deletions,
      "  insertadas:", cuentas.insertions)
print(f"WER: {cuentas.wer:.3f}")

import io

import pandas as pd
from huggingface_hub import hf_hub_download


def con_ruido(audio, snr_db, semilla=0):
    """Suma ruido gaussiano con la relación señal-ruido pedida."""
    rng = np.random.default_rng(semilla)
    potencia_ruido = np.mean(audio ** 2) / 10 ** (snr_db / 10)
    ruido = rng.normal(0, np.sqrt(potencia_ruido), len(audio))
    return (audio + ruido).astype(np.float32)


# 10 frases de LibriSpeech (Panayotov et al., 2015), licencia CC BY 4.0
parquet = hf_hub_download("hf-internal-testing/librispeech_asr_dummy",
                          "clean/validation-00000-of-00001.parquet",
                          repo_type="dataset")
frases = pd.read_parquet(parquet).head(10)
ingles = [sf.read(io.BytesIO(a["bytes"]), dtype="float32")[0]
          for a in frases["audio"]]
refs_en = [normalizar(texto) for texto in frases["text"]]
print("Frases en inglés:", len(ingles), " duración total:",
      f"{sum(len(x) for x in ingles) / 16000:.0f} s")
mostrar("Ejemplo: " + refs_en[0])

textos, wer = {}, {}
for nombre in ["tiny", "base"]:
    textos[nombre, "es"] = transcribir(onda, nombre, "es")
    textos[nombre, "es + ruido"] = transcribir(con_ruido(onda, 10),
                                               nombre, "es")
    textos[nombre, "en"] = [transcribir(x, nombre, "en")
                            for x in ingles]
    textos[nombre, "en + ruido"] = [
        transcribir(con_ruido(x, 10, semilla=i), nombre, "en")
        for i, x in enumerate(ingles)]
    for caso in ["es", "es + ruido"]:
        wer[nombre, caso] = jiwer.wer(ref, normalizar(textos[nombre, caso]))
    for caso in ["en", "en + ruido"]:
        hips = [normalizar(t) for t in textos[nombre, caso]]
        wer[nombre, caso] = jiwer.wer(refs_en, hips)

casos = ["es", "es + ruido", "en", "en + ruido"]
print("modelo " + "".join(f"{c:>12}" for c in casos))
for nombre in ["tiny", "base"]:
    print(f"{nombre:6} " + "".join(f"{wer[nombre, c]:12.3f}"
                                  for c in casos))

posiciones = np.arange(len(casos))
for desplazar, nombre, color in [(-0.2, "tiny", "C2"),
                                 (0.2, "base", "C0")]:
    plt.bar(posiciones + desplazar, [wer[nombre, c] for c in casos],
            width=0.4, color=color, label=f"whisper-{nombre}")
plt.xticks(posiciones, ["Español", "Español\ncon ruido", "Inglés",
                        "Inglés\ncon ruido"], fontsize=13)
plt.ylabel("WER (menos es mejor)")
plt.title("Whisper con y sin ruido (SNR 10 dB)")
plt.legend()
plt.show()