Source: https://jmrplens.github.io/phonometry/es/guides/cepstrum-echoes/
Los [estimadores espectrales](/phonometry/es/guides/spectral-analysis/) 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](/phonometry/es/guides/correlation-delay/) 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.
```python
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
```python
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`](/phonometry/es/guides/swept-sine-distortion/).
## 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).
```python
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
```python
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:
```python
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
```python
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`:
```python
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`](/phonometry/es/guides/swept-sine-distortion/) 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](/phonometry/es/guides/correlation-delay/)
(`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` |
```python
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
```python
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](/phonometry/es/guides/spectral-analysis/),
y su núcleo de plegado es literalmente el que hay dentro de
[`minimum_phase`](/phonometry/es/guides/swept-sine-distortion/) - 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](/phonometry/es/guides/correlation-delay/) devuelve en
el tiempo, y un preanálisis natural antes de las métricas dedicadas de
[modulación de amplitud de aerogeneradores](/phonometry/es/guides/wind-turbine-noise/):
el espectro de la envolvente dice *si hay modulación y a qué ritmo*, las
métricas del dominio la cuantifican normativamente.