Cepstro, ecos y espectro de la envolvente

Los estimadores espectrales describen qué frecuencias contiene una señal; esta página cubre lo que se esconde en la forma de ese espectro. El cepstro - la transformada de Fourier inversa del espectro logarítmico - vive en phonometry.metrology y convierte dos problemas espectrales difíciles en simple localización de picos: el rizado espectral periódico (un eco, una familia de armónicos) colapsa en un único pico en la quefrencia de su período, y la envolvente espectral suave se separa de la estructura fina mediante un simple enventanado - liftering - en el dominio de la quefrencia. La misma maquinaria extiende la envolvente de Hilbert con un espectro de la envolvente en el que las modulaciones de amplitud se vuelven líneas discretas.

1. El cepstro y sus tres variantes

Como el logaritmo convierte la convolución x = h * u en la suma ln X = ln H + ln U, componentes que se solapan en el espectro se suman - y se separan - en el dominio cepstral (Havelock cap. 27, Ecs. (22)-(23)). cepstrum calcula las tres variantes estándar sobre el eje de quefrencia:

• 'power' (por defecto): la DFT inversa de ln|X|² (Fig. 21 de Milner). Real, par y ciega a la fase - el caballo de batalla para detectar ecos y familias de armónicos. Es el cepstro con signo del espectro logarítmico de potencia según Milner; el “cepstro de potencia” original de Bogert (1963) eleva al cuadrado una vez más y es no negativo - la biblioteca sigue a Milner en todo, así que los rahmónicos negativos conservan su signo;
• 'real': la DFT inversa de ln|X| - exactamente la mitad del cepstro de potencia, y la cantidad cuyo plegado causal es la reconstrucción de fase mínima (más abajo);
• 'complex': la DFT inversa de ln|X| + j·arg X con la fase desenrollada y su componente lineal eliminada (Havelock cap. 87, Ec. (14)). Conserva la fase, así que es invertible: la puerta de entrada a la deconvolución homomórfica.

import numpy as np
from phonometry import cepstrum

fs = 48000.0
rng = np.random.default_rng(1)
x = rng.standard_normal(4096)

res = cepstrum(x, fs, kind="power")
print(res.quefrencies[:3], res.cepstrum.shape)   # eje de quefrencia, en s
res.plot(language="es")

Las tres variantes de un mismo registro con eco (una ondícula de banda limitada más una reflexión en 8 ms, a = 0,5). Las tres llevan los rahmónicos en 8 y 16 ms con alturas a y −a²/2 (convención con signo de Milner); el recuadro muestra el cepstro real exactamente a la mitad del de potencia, y la ondícula fuente se concentra por debajo de 2 ms.

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sp_signal
from phonometry import cepstrum

fs = 48000.0
# Una ondícula fuente de banda limitada más una reflexión en 8 ms (a = 0,5)
b, a = sp_signal.butter(2, 0.3)
s = np.zeros(4096)
s[37:37 + 256] = sp_signal.lfilter(b, a, np.r_[1.0, np.zeros(255)])
x = s + 0.5 * np.roll(s, 384)

# Una línea por variante: cada CepstrumResult se dibuja solo:
cepstrum(x, fs, kind="power").plot(language="es")
plt.show()

# Las tres variantes superpuestas a mano en un mismo eje de quefrencia:
fig, ax = plt.subplots()
for kind, style in (("power", "-"), ("real", "--"), ("complex", ":")):
    res = cepstrum(x, fs, kind=kind)
    q_ms = 1e3 * res.quefrencies
    mask = (q_ms > 0.5) & (q_ms <= 20.0)
    ax.plot(q_ms[mask], res.cepstrum[mask], style, label=f"cepstro {kind}")
ax.set(xlabel="Quefrencia [ms]", ylabel="Cepstro")
ax.legend()
plt.show()

El resultado lleva el eje de quefrencia periódico completo (0 .. (nfft-1)/fs); las quefrencias por encima de nfft/(2·fs) son las quefrencias negativas especulares, donde los cepstros de potencia y real (que son pares) se repiten y el cepstro complejo guarda su contenido anticausal (de fase no mínima). El relleno de ceros vía nfft reduce el aliasing temporal cepstral cuando el espectro logarítmico tiene rasgos abruptos, exactamente igual que el relleno oversample de minimum_phase.

2. Detección de ecos: el tren de picos rahmónicos

Una reflexión única x(t) = s(t) + a·s(t-t0) multiplica el espectro por 1 + a·e^{-j2πft0} - un rizado de período 1/t0 en toda la banda. Su logaritmo se expande, para |a| < 1, en la serie exactamente sumable

$$\ln \left(1+a{e}^{-j\theta }\right)=\sum _{n\ge 1}(-1{)}^{n+1}\frac{a^n}{n}{e}^{-jn\theta },$$

así que el cepstro lleva un tren de picos en los rahmónicos n·t0 con amplitudes a, -a²/2, a³/3, ... (su suma es ln(1+a)), independientemente del espectro de la propia s, que se concentra en las quefrencias bajas. En el cepstro con signo del espectro logarítmico de potencia (kind='power', la convención de Milner) la altura del primer pico es el coeficiente de reflexión a - con su signo - más lo que el cepstro de la fuente aporte en esa quefrencia (despreciable para fuentes de banda ancha, cuyo cepstro se concentra en las quefrencias bajas); sobre un impulso ideal con eco la identidad es una forma cerrada que los tests y la batería de conformidad fijan a 1e-10. echo_detection automatiza la lectura: busca el mayor pico de |cepstrum| en la banda de búsqueda (así una reflexión inversora, a < 0, se encuentra en su retardo verdadero en vez de perderse), refina el retardo por interpolación cuadrática a través del pico y sus vecinos, y devuelve como coeficiente de reflexión el valor del pico con su signo. Cuando el retardo verdadero cae entre muestras, el rahmónico se reparte entre bins de quefrencia vecinos: el retardo interpolado sigue acertando, pero el coeficiente devuelto subestima |a| (hasta cerca del 65 % de su valor a mitad de camino entre muestras).

import numpy as np
from phonometry import echo_detection

fs = 48000.0
rng = np.random.default_rng(2)
s = rng.standard_normal(12000)               # fuente de banda ancha
x = s + 0.5 * np.roll(s, 384)                # eco: 8 ms, a = 0.5

res = echo_detection(x, fs, min_quefrency=0.002)
print(res.delay, res.reflection_coefficient)  # 0.008 s, ~0.5
res.plot(language="es")

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sp_signal
from phonometry import echo_detection, noise_signal

fs = 48000.0
n = 12000
impulse = np.zeros(n)
impulse[0] = 1.0
b, a = sp_signal.butter(2, [0.004, 0.9], btype="bandpass")
direct = sp_signal.lfilter(b, a, impulse)              # clic de banda ancha
ir = direct + 0.5 * np.roll(direct, int(0.008 * fs))   # eco a 8 ms
ir += noise_signal(fs, n / fs, color="white", rms=1e-4, seed=13)

res = echo_detection(ir, fs, min_quefrency=0.002)

fig, ax = plt.subplots(figsize=(10, 6))
half = res.nfft // 2 + 1
ax.plot(1e3 * res.quefrencies[:half], res.cepstrum[:half], lw=1.1)
ax.axvline(8.0, ls="--", color="k", label="Retardo verdadero del eco")
ax.plot([1e3 * res.delay], [res.reflection_coefficient], "v", ms=10,
        label="Pico detectado (altura = reflexión a)")
ax.set_xlim(0.0, 30.0)
ax.set_xlabel("Quefrencia [ms]")
ax.set_ylabel("Cepstro")
ax.legend()
plt.show()

La banda de búsqueda empieza por encima de la región de quefrencias bajas que ocupa la envolvente espectral de la propia fuente (min_quefrency, por defecto 16 muestras) y termina en la mitad no ambigua del eje (max_quefrency). Los trenes de picos por reverberación de la sismología (Havelock cap. 87) son la misma firma a escala geofísica. Obsérvese el segundo rahmónico negativo en 2·t0 de la figura: el término -a²/2 de la serie, una confirmación útil de que un pico es de verdad un eco y no una periodicidad espectral sin relación.

3. Liftering: envolvente frente a estructura fina

Filtrar en el dominio de la quefrencia se llama liftering (Havelock cap. 27, sec. 4.3). Un lifter paso bajo conserva las quefrencias por debajo del corte y devuelve la envolvente log-espectral suave con el rizado eliminado; un lifter paso alto conserva el complemento - solo el rizado. Los dos modos son exactamente complementarios en dB, porque la partición es lineal en el dominio logarítmico:

import numpy as np
from phonometry import lifter

fs = 48000.0
rng = np.random.default_rng(3)
s = rng.standard_normal(12000)
x = s + 0.5 * np.roll(s, 384)                # el mismo eco de 8 ms

low = lifter(x, fs, cutoff=0.004, mode="lowpass")    # envolvente de ln|X|
high = lifter(x, fs, cutoff=0.004, mode="highpass")  # el rizado del eco
print(np.allclose(low.liftered_db + high.liftered_db, low.spectrum_db))
low.plot(language="es")

El reparto del lifter de 4 ms sobre un registro de eco puro: el lado paso bajo devuelve la envolvente espectral suave, el lado paso alto aísla el rizado de 125 Hz del eco, oscilando exactamente entre las formas cerradas 20·lg(1±a) = +3,5 y −6,0 dB.

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sp_signal
from phonometry import lifter

fs = 48000.0
# Una ondícula de banda limitada con un eco puro de 8 ms (a = 0,5), de modo
# que el rizado del paso alto tiene las cotas exactas de forma cerrada.
b, a = sp_signal.butter(2, 0.3)
s = np.zeros(4096)
s[37:37 + 256] = sp_signal.lfilter(b, a, np.r_[1.0, np.zeros(255)])
x = s + 0.5 * np.roll(s, 384)

low = lifter(x, fs, cutoff=0.004, mode="lowpass")
high = lifter(x, fs, cutoff=0.004, mode="highpass")

# Una línea: el cepstro con el corte y los dos espectros logarítmicos:
low.plot(language="es")
plt.show()

# El reparto envolvente/rizado a mano, ampliado a 500-2000 Hz:
band = (low.frequencies >= 500) & (low.frequencies <= 2000)
fig, axes = plt.subplots(2, 1, sharex=True)
axes[0].semilogx(low.frequencies[band], low.spectrum_db[band], "0.6",
                 lw=0.7, label="Espectro logarítmico")
axes[0].semilogx(low.frequencies[band], low.liftered_db[band], lw=2,
                 label="Lifter paso bajo: envolvente")
axes[1].semilogx(high.frequencies[band], high.liftered_db[band], "r",
                 label="Lifter paso alto: rizado")
for bound in (20 * np.log10(1.5), 20 * np.log10(0.5)):
    axes[1].axhline(bound, color="g", linestyle="--")
axes[1].set_xlabel("Frecuencia [Hz]")
for ax in axes:
    ax.set_ylabel("Magnitud [dB]")
    ax.legend()
plt.show()

Para la señal de eco puro el rizado del paso alto oscila entre las formas cerradas 20·log10(1+a) y 20·log10(1-a) dB, otro oráculo que fijan los tests. En análisis de voz la misma operación separa la envolvente del tracto vocal (formantes) de los armónicos de la excitación; aquí es la herramienta general para separar “lo suave de lo periódico” en cualquier respuesta en magnitud medida.

4. El cepstro complejo y la conexión con la fase mínima

El cepstro complejo conserva la fase desenrollada, así que la transformada es de ida y vuelta: CepstrumResult.invert() restaura el registro con precisión de máquina, incluida la componente de fase lineal (retardo puro) que la transformada directa elimina y guarda en linear_phase_samples:

import numpy as np
from scipy import signal as sp_signal
from phonometry import cepstrum

fs = 48000.0
x = np.zeros(2048)
b, a = sp_signal.butter(2, 0.3)
x[37:293] = sp_signal.lfilter(b, a, np.r_[1.0, np.zeros(255)])

res = cepstrum(x, fs, kind="complex")
print(res.linear_phase_samples)              # negativo: retardo eliminado
print(np.max(np.abs(res.invert() - x)))      # ~1e-14
res.plot(language="es")   # el cepstro complejo frente a la quefrencia

Entre el logaritmo y la transformada inversa se puede editar lo que se quiera - eso es la deconvolución homomórfica (Havelock cap. 87, sec. 3.3): anular los rahmónicos elimina un eco, conservar solo las quefrencias bajas extrae la ondícula de la fuente. Una señal de fase mínima tiene un cepstro complejo causal, y por eso plegar el cepstro real sobre las quefrencias positivas reconstruye la fase mínima a partir de |H| únicamente: minimum_phase y phase_decomposition corren sobre ese mismo núcleo de plegado (Bendat y Piersol sec. 13.1.4; Tohyama, en Havelock cap. 75, edita la reverberación manipulando exactamente estas partes causal y anticausal).

5. El espectro de la envolvente: modulaciones como líneas

Donde el cepstro encuentra periodicidades del espectro, el espectro de la envolvente encuentra periodicidades de la amplitud. La sección 13.3 de Bendat y Piersol (Fig. 13.11) formaliza la estructura: un detector de envolvente, un eliminador de continua y una vista espectral de lo que queda. envelope_spectrum pasa la envolvente de Hilbert (kind="magnitude", el valor práctico por defecto) o el detector cuadrático del libro (kind="squared") por exactamente esa cadena, escalada por la ganancia coherente de la ventana para que una modulación sinusoidal cuya frecuencia cae en un bin de análisis se lea como una línea con su amplitud exacta (las líneas fuera de bin leen de menos por la pérdida de scalloping de la ventana, hasta cerca de 1,4 dB con la Hann por defecto). El argumento opcional band=(low, high) reproduce el filtro paso banda de entrada de la figura - la cadena clásica de envolvente para rodamientos: aislar con un paso banda de fase cero la banda de resonancia estructural que excitan los impactos del defecto y después obtener su envolvente - de modo que un interferente fuera de banda llega muy atenuado al detector (la caída de un Butterworth de cuarto orden aplicado hacia delante y hacia atrás; el rechazo es finito, y el paso de fase cero deja pequeños transitorios en los bordes del registro).

Para un tono AM A0·(1 + m·cos(2πfm·t))·cos(2πfc·t) con fm en un bin de análisis las formas cerradas son:

kind	nivel medio	línea en fm	línea en 2fm
'magnitude'	A0	A0·m	-
'squared'	A0²·(1 + m²/2)	2·A0²·m	A0²·m²/2

import numpy as np
from phonometry import envelope_spectrum

fs = 8192.0
t = np.arange(int(4 * fs)) / fs
x = (1.0 + 0.4 * np.cos(2 * np.pi * 25.0 * t)) * np.cos(2 * np.pi * 1000.0 * t)

res = envelope_spectrum(x, fs)
k = int(round(25.0 * res.nfft / fs))
print(res.mean_level, res.amplitude[k])      # ~1.0 y ~0.4
res.plot(language="es")

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from phonometry import envelope_spectrum, noise_signal

fs = 8192.0
seconds = 4.0
t = np.arange(int(seconds * fs)) / fs
x = (1.0 + 0.4 * np.cos(2 * np.pi * 25.0 * t)) * np.cos(2 * np.pi * 1000.0 * t)
x += noise_signal(fs, seconds, color="white", rms=0.03, seed=8)

res = envelope_spectrum(x, fs)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(res.frequencies, res.amplitude, lw=1.4,
        label="Espectro de la envolvente")
ax.axvline(25.0, ls="--", color="k", label="Frecuencia de modulación")
ax.axhline(0.4, ls=":", color="r",
           label=r"Amplitud exacta de la línea $A_0 m$")
ax.set_xlim(0.0, 100.0)
ax.set_xlabel("Frecuencia [Hz]")
ax.set_ylabel("Amplitud de modulación")
ax.legend()
plt.show()

Los defectos de rodamientos y engranajes, el zumbido de la red eléctrica y la modulación de amplitud de los aerogeneradores aparecen así: líneas en la frecuencia de modulación y sus armónicos, limpiamente separadas del espectro de la propia portadora. La media de la envolvente que elimina la etapa de continua se conserva en mean_level (la amplitud de la portadora para el detector de magnitud), y remove_dc=False omite el eliminador cuando importa la línea de continua absoluta.

Qué cubre esta guía

Cubierto. Las tres variantes del cepstro y el liftering del Handbook of Signal Processing in Acoustics de Havelock, Kuwano y Vorländer (cepstrum, lifter, capítulos 27 y 87), el retardo de un único eco y el coeficiente de reflexión leídos en el cepstro de potencia (echo_detection), el cepstro complejo invertible y su ida y vuelta homomórfica (CepstrumResult.invert), y el espectro de la envolvente del capítulo 13 de Bendat y Piersol (envelope_spectrum) con sus líneas de forma cerrada para un tono AM.

No cubierto. echo_detection solo elige el mayor pico cepstral de la banda de búsqueda, así que una respuesta con varios ecos solapados necesita localización manual de picos o llamadas repetidas sobre bandas más estrechas: no es un separador multi-eco. El módulo tampoco tiene un cepstro en escala mel ni de tipo MFCC para rasgos de audio perceptual: lifter y cepstrum trabajan solo sobre el espectro logarítmico en frecuencia lineal.

Relación con los demás estimadores

El cepstro parte de las mismas convenciones FFT de registro único que los estimadores espectrales calibrados, y su núcleo de plegado es literalmente el que hay dentro de minimum_phase - la refactorización queda fijada bit a bit en los tests. El espectro de la envolvente es la vista en frecuencia de la misma señal analítica que la envolvente de Hilbert devuelve en el tiempo, y un preanálisis natural antes de las métricas dedicadas de modulación de amplitud de aerogeneradores: el espectro de la envolvente dice si hay modulación y a qué ritmo, las métricas del dominio la cuantifican normativamente.

Referencias

• Bendat, J. S. y Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4.ª ed.). Wiley. https://doi.org/10.1002/9781118032428Sección 13.1.4 (la relación de Hilbert entre log-magnitud y fase tras el plegado de fase mínima) y sección 13.3 con la figura 13.11 (detección de envolvente seguida de eliminación de la continua, la estructura del espectro de la envolvente). ISBN 978-0-470-24877-5.
• Havelock, D., Kuwano, S. y Vorländer, M. (eds.). (2008). Handbook of signal processing in acoustics. Springer. https://doi.org/10.1007/978-0-387-30441-0Capítulo 27 (Milner: la transformada cepstral como DFT inversa del espectro logarítmico de potencia, la quefrencia, liftering paso bajo/paso alto), capítulo 87 (Neelamani: el cepstro complejo, deconvolución homomórfica, trenes de picos periódicos por reverberación) y capítulo 75 (Tohyama: manipulación fase mínima/pasa-todo en el dominio cepstral). ISBN 978-0-387-77698-9.