Análisis tiempo-frecuencia

Referencias: Bendat y Piersol 2010 Harris 1978

Un espectro estacionario esconde todo lo que ocurre en el tiempo: una sirena que pasa, un impacto, una máquina arrancando. Esta página cubre los dos estimadores tiempo-frecuencia de phonometry.signals, ambos con la disciplina de calibración de la página de análisis espectral: el espectrograma calibrado (la vista de transformada de Fourier de tiempo corto de Bendat y Piersol, sección 12.6.4.2, en unidades absolutas: dB SPL para una señal en pascales) y la FFT con zoom (sección 11.5.4), que calcula el espectro de una banda estrecha en una malla arbitrariamente fina para separar tonos más próximos que un bin práctico de FFT de banda completa. Donde la página de niveles ofrece el espectrograma en bandas de fracción de octava de un sonómetro, esta es su equivalente de banda fina y ancho de banda constante.

Un solo número gobierna todo en esta página: la longitud del segmento, que trocea el plano tiempo-frecuencia en celdas de área fija. El diagrama muestra el mismo registro teselado por una ventana corta y una larga.

Dos teselados del plano tiempo-frecuencia del mismo registro a 16 kilohercios lado a lado: con una ventana corta de 256 muestras las celdas miden 16 milisegundos de ancho y 62,5 hercios de alto, así que un clic queda como una columna estrecha y nítida mientras un tono estable se emborrona en una banda alta; con una ventana larga de 1024 muestras las celdas miden 64 milisegundos de ancho y 15,6 hercios de alto, así que el tono se afila en una banda fina mientras el clic se emborrona en una columna ancha; los pies anotan que cada celda es una única estimación sin promediar con error aleatorio unidad y que el producto de resoluciones temporal y frecuencial se mantiene en torno a uno

1. El espectrograma calibrado

spectrogram divide el registro en segmentos con ventana (Hann por defecto) y solapados - exactamente la segmentación de power_spectral_density - y conserva el periodograma unilateral de cada segmento como una columna de la representación tiempo-frecuencia, en lugar de promediarlos. Como la calibración es el escalado exacto del módulo de Welch sin eliminación de tendencia, se cumplen tres identidades:

una señal en pascales da Pa²/Hz ('density') o Pa² ('spectrum') por celda, así que $10 lo g_{10} (power / p_{0}^{2})$ con escalado 'spectrum' lee directamente el nivel de presión acústica de un tono en cualquier columna que abarque;
la media de las columnas en el tiempo reproduce power_spectral_density bin a bin, con la misma ventana, solape y escalado;
con una ventana cuyo cuadrado suma una constante al solaparse (Hann al 75 % de solape), la potencia 'density' integrada en el tiempo es exactamente la energía del registro en el dominio temporal (Parseval más la identidad COLA, una de las comprobaciones de conformidad).

from phonometry import spectrogram

res = spectrogram(x, fs, nperseg=1024, overlap=0.75, scaling="spectrum")
print(res.power.shape)             # (frecuencias, tiempos)
print(res.time_resolution)         # T_B = nperseg/fs, en s
print(res.resolution_bandwidth)    # Be del segmento con ventana, en Hz
res.plot(language="es")            # imagen en dB en el plano tiempo-frecuencia

La longitud del segmento es toda la decisión de diseño: $T_{B} = nperseg / f_{s}$ de resolución temporal contra $B_{e} = ENBW / T_{B}$ de resolución frecuencial, así que su producto $T_{B} B_{e}$ lo fija solo la ventana: 1 para una ventana rectangular y 1,5 para la Hann por defecto que se usa aquí (sección 12.6.4.2, y las figuras de mérito de ventana de la página de análisis espectral). El diagrama de teselado de arriba dibuja el caso sin ventana, donde ese producto vale exactamente uno. Partir por la mitad la longitud del segmento compra el doble de resolución temporal y cuesta exactamente la mitad de resolución frecuencial; nada compra las dos. Los segmentos largos fijan las frecuencias y emborronan los transitorios; los cortos, lo contrario. Y como cada celda es una única estimación sin promediar, los datos aleatorios llevan un error aleatorio normalizado de 1 por celda (ec. 8.158 con $n_{d} = 1$ ; Bendat y Piersol citan $2 /1, 25 \approx 1, 13$ para la representación en magnitud): el espectrograma es una herramienta para la estructura determinista (tonos, barridos, transitorios), mientras que la estimación de Welch promediada es la herramienta de baja varianza para el fondo estacionario.

Ponle un número. Una única celda de periodograma de datos aleatorios es una ji-cuadrado con dos grados de libertad, así que su nivel tiene una desviación típica de unos 5,6 dB y una dispersión al 90 % de unos 17 dB: la textura moteada del fondo de ruido de la figura de abajo, y la razón de que una sola celda brillante no signifique absolutamente nada. De ahí sale la regla de lectura: un rasgo es real cuando persiste en columnas o bins vecinos, porque la fluctuación de celda a celda es independiente y la estructura no, de modo que el ojo hace el promediado que el estimador se negó a hacer. Los dos remedios cuestan algo: promediar columnas contiguas cambia resolución temporal por varianza y converge hacia la estimación de Welch, y suavizar en frecuencia con el núcleo de fracción de octava de la página de análisis espectral cambia resolución frecuencial por lo mismo. Y por eso mismo un fondo de ruido se cuantifica con una PSD y solo se representa en un espectrograma.

Elegir el segmento a partir del evento

El compromiso está fijado, así que la elección tiene que venir de aquello que se mide.

Para un transitorio, mantén $T_{B}$ por debajo de la caída que quieres ver. Los 12 ms de caída del impacto de la figura de abajo piden segmentos de unos pocos milisegundos, unas 256 muestras a 48 kHz, y la resolución frecuencial sale entonces en unos 190 Hz, la quieras o no.

Para una componente cuya frecuencia se mueve a un ritmo $r$ en Hz/s, el tono se sale de su propio ancho de banda de resolución dentro de un mismo segmento salvo que $B_{e} \geq r T_{B}$ . Con $B_{e} = 1, 5/ T_{B}$ eso da un óptimo $T_{B} = 1, 5/ r$ . La sirena de abajo barre hasta unos 950 Hz/s, así que segmentos de menos de 40 ms mantienen la cresta afilada, y la elección de 1024 puntos a 16 kHz son 64 ms, ya pasado ese límite, que es por lo que la cresta engorda en la parte más empinada del barrido. Una subida de régimen que barre un orden de 500 Hz a 100 Hz/s pide segmentos de unos 0,12 s y nada más largo.

Los dos requisitos entran en conflicto a menudo, como pasa en la propia escena de la figura. La receta honesta: elige $T_{B}$ a partir de lo más rápido que tenga que quedar nítido, contrasta el $B_{e}$ resultante con las dos componentes más próximas que tengan que quedar separadas y, cuando no se puedan cumplir las dos cosas, dilo en el pie en lugar de partir la diferencia. res.time_resolution y res.resolution_bandwidth informan del par realmente entregado, así que la elección se puede declarar en un informe.

Una nota de adquisición, porque es el fallo que parece un resultado: captura los transitorios con un predisparo de al menos un segmento, fija la ganancia de antemano y verifica que ninguna muestra recorta; un impacto recortado dibuja una franja de altura completa que parece exactamente un evento de banda ancha.

Espectrograma calibrado en dB SPL de una escena exterior sintética de cuatro segundos: una sirena que barre sinusoidalmente entre 600 y 1200 hercios traza una cresta oscilante brillante a 70 decibelios, un impacto de banda ancha a los dos segundos y medio dibuja una franja vertical y un fondo de ruido rosa cuyo nivel de banda ancha es de 45 decibelios rellena la imagen con un nivel por celda mucho más bajo, con la barra de color leyendo el nivel de presión acústica absoluto

Cómo leerla: la cresta de la sirena marca de verdad 70 dB SPL, porque un tono mete toda su potencia dentro de un ancho de banda de resolución y la celda 'spectrum' es por tanto su valor cuadrático medio. El fondo rosa no marca 45 dB en ninguna parte. Esos 45 dB son su nivel de banda ancha repartido por toda la banda, y cada celda solo lleva la parte que cae dentro de un $B_{e}$ : unos 25 dB por debajo aquí, para un segmento de 1024 puntos a 16 kHz. Para recuperar un nivel a partir del ruido, convierte las celdas en densidad restando $10 lo g_{10} B_{e}$ e integra en la banda que interese, o usa 'density' e integra directamente. Y la celda de pico del impacto tampoco es el nivel del impacto: un transitorio más corto que el segmento reparte su energía por toda la columna, así que esa lectura depende de nperseg. La barra de color es una escala absoluta para la estructura tonal y una escala dependiente del ancho de banda para todo lo demás; cuando lo que se busca son niveles de banda absolutos a lo largo del tiempo, la herramienta correcta es el espectrograma de bandas de Niveles.

Ver el código de esta figura

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

fs = 16000.0
t = np.arange(int(4.0 * fs)) / fs
p_ref = 2e-5

siren_rms = p_ref * 10.0 ** (70.0 / 20.0)          # sirena a 70 dB SPL
x = siren_rms * np.sqrt(2.0) * np.cos(
    2.0 * np.pi * 900.0 * t - 600.0 * np.cos(np.pi * t)
)
rng = np.random.default_rng(9)
n_imp = int(0.06 * fs)                              # impacto en t = 2,5 s
x[int(2.5 * fs):int(2.5 * fs) + n_imp] += 0.4 * (
    rng.standard_normal(n_imp) * np.exp(-np.arange(n_imp) / (0.012 * fs))
)
x += noise_signal(fs, 4.0, color="pink",            # fondo a 45 dB SPL
                  rms=p_ref * 10.0 ** (45.0 / 20.0), seed=10)

res = spectrogram(x, fs, nperseg=1024, overlap=0.75, scaling="spectrum")
level = 10.0 * np.log10(res.power / p_ref**2)

fig, ax = plt.subplots(figsize=(10, 6))
img = ax.imshow(level, cmap="magma", vmin=level.max() - 55.0,
                vmax=level.max(), aspect="auto", origin="lower",
                extent=(res.times[0], res.times[-1], 0.0,
                        res.frequencies[-1]))
fig.colorbar(img, ax=ax, label="Nivel de presión acústica [dB SPL]")
ax.set_ylim(0.0, 3000.0)
ax.set_xlabel("Tiempo [s]")
ax.set_ylabel("Frecuencia [Hz]")
plt.show()

res.plot() dibuja la misma representación directamente desde el resultado (una única imagen ráster, 80 dB por debajo de la celda más fuerte por defecto; pasa vmin/vmax para cambiar el rango, o ax para dibujar en un panel existente).

2. FFT con zoom

Dos tonos separados 3 Hz - bandas laterales de engranajes, máquinas gemelas, zumbido de red contra un armónico de rotor - son invisibles para una FFT de 1024 puntos a 8192 Hz: sus bins miden 8 Hz. La solución clásica de los analizadores es la transformada con zoom de Bendat y Piersol, sección 11.5.4: filtrar la banda, bajarla a frecuencia cero por demodulación compleja con $exp (- j 2 π f_{1} t)$ (ecs. 11.123-11.126), diezmar por la razón de anchos de banda y transformar el registro diezmado (ecs. 11.128-11.130), obteniendo un espaciado fino de bins sobre la banda sin un bloque de FFT gigante (ec. 11.127).

zoom_fft calcula el equivalente digital exacto en una sola pasada - la evaluación chirp-Z de la DFT del registro con ventana sobre la malla del zoom - que da las mismas muestras de la DFT que la cadena demodular-diezmar; la batería de tests fija ambas entre sí a precisión de máquina. Las amplitudes se calibran con la ganancia coherente de la ventana, así que una sinusoide de amplitud de pico $A$ sobre una frecuencia de análisis lee amplitude $= A$ y power $= A^{2} /2$ exactamente:

from phonometry import zoom_fft

res = zoom_fft(x, fs, 980.0, 1016.0)     # malla a la resolución del registro
print(res.bin_spacing)                   # fs/N por defecto
print(res.resolution_bandwidth)          # Be del registro con ventana
peak = res.amplitude.argmax()
print(res.frequencies[peak], res.amplitude[peak])
res.plot(language="es")                  # espectro de potencia en dB

Una distinción importa y el resultado la declara: la malla puede hacerse arbitrariamente fina (n_points), pero la resolución - la capacidad de separar dos tonos - la fijan la longitud del registro y la ventana, y se informa como resolution_bandwidth ( $B_{e} = f_{s} \sum w^{2} / (\sum w)^{2}$ , es decir, $1/ T_{r}$ sin ventana y $1, 5/ T_{r}$ con Hann; el mismo factor de ventana del apartado 1, aplicado aquí a la duración del registro $T_{r}$ en lugar de a un segmento, porque la FFT con zoom transforma el registro entero). El zoom refina el muestreo del mismo espectro subyacente; solo un registro más largo separa tonos más próximos.

El rango dinámico es otra cuestión, y el zoom no ayuda con ella. La chirp-Z sigue evaluando la DFT del registro entero con su ventana, así que una componente fuerte fuera de la banda ampliada se cuela dentro por los lóbulos laterales de la ventana exactamente igual que en una FFT corriente: la banda ampliada compra resolución de la malla, no rechazo. El primer lóbulo lateral de la Hann está a −31,5 dB, que no da ni de lejos para ver una banda lateral 50 dB por debajo de una portadora a unos pocos hercios. El control es el argumento window (por defecto "hann"): pasa window="blackman", o una Kaiser con beta alta, para cambiar anchura de lóbulo principal por rechazo de lóbulos laterales, recordando que un lóbulo principal más ancho sube $B_{e}$ y con él la separación mínima que puedes resolver, de modo que los dos objetivos tiran en sentidos opuestos. Las figuras de mérito de las ventanas ponen números a ese compromiso. Cuando ni la mejor ventana basta, filtra o elimina con un notch la componente dominante antes de ampliar: ninguna ventana sustituye a quitar la energía.

Y la máquina tiene que estarse quieta. Separar dos rayas a 3 Hz una de otra cerca de 1 kHz pide un segundo de registro aproximadamente, y durante ese segundo el eje tiene que mantener su velocidad con un error menor que 3 partes en 1000, el 0,3 %, o la propia raya del rotor se emborrona más ancha que la separación que intentas resolver, y ningún zoom la recupera. En general la estabilidad de velocidad exigida es la resolución objetivo dividida por la frecuencia de la raya, lo que se vuelve severo enseguida: 0,5 Hz a 5 kHz es el 0,01 %. Compruébalo antes de fiarte del zoom: pasa un espectrograma del mismo registro y confirma que la raya sale horizontal a lo largo de la ventana de análisis, o lee el periodo del eje con un tacómetro durante el registro. Cuando la máquina sí deriva, remuestrea angularmente contra el tacómetro (seguimiento de órdenes) para que las rayas queden estacionarias en orden en vez de en frecuencia.

FFT con zoom de un registro de un segundo con tonos a 997 y 1000 hercios: la FFT gruesa de 1024 puntos con bins de 8 hercios muestra un único bulto ancho, mientras que la FFT con zoom entre 980 y 1016 hercios dibuja dos lóbulos principales separados cuyos picos caen exactamente sobre las frecuencias verdaderas punteadas

Lo que separa los dos tonos es el segundo de registro, no el zoom. La vista gruesa usa 1024 de las 8192 muestras, así que sus bins de 8 Hz no pueden sostener dos rayas separadas 3 Hz y enseña un solo bulto en 1000 Hz. El zoom transforma el registro entero, cuyo ancho de banda de resolución con Hann es de 1,5 Hz, y le pone encima una malla de 0,25 Hz: los picos caen en 997,00 y 1000,00 Hz con −1,94 y −6,02 dB, las amplitudes de los dos tonos. Una malla más fina sobre un registro de 1024 muestras no habría cambiado nada; esta es la distinción que hace el apartado entre densidad de muestreo y resolución.

Ver el código de esta figura

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

fs = 8192.0
t = np.arange(8192) / fs                 # registro de 1 s: resolución de 1 Hz
x = (0.8 * np.cos(2.0 * np.pi * 997.0 * t)
     + 0.5 * np.cos(2.0 * np.pi * 1000.0 * t))

w = signal.get_window("hann", 1024)      # la vista gruesa: bins de 8 Hz
coarse = 2.0 * np.abs(np.fft.rfft(x[:1024] * w)) / np.sum(w)
coarse_f = np.fft.rfftfreq(1024, 1.0 / fs)
band = (coarse_f >= 950.0) & (coarse_f <= 1050.0)

res = zoom_fft(x, fs, 980.0, 1016.0, n_points=145)   # malla de 0,25 Hz

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(coarse_f[band], 20.0 * np.log10(coarse[band]), "o--",
        label="FFT de 1024 puntos (bins de 8 Hz)")
ax.plot(res.frequencies, 20.0 * np.log10(res.amplitude),
        label="FFT con zoom del mismo registro")
for f0 in (997.0, 1000.0):
    ax.axvline(f0, color="k", ls=":", lw=1.0, alpha=0.6)
ax.set_ylim(-70.0, 5.0)
ax.set_xlabel("Frecuencia [Hz]")
ax.set_ylabel("Amplitud [dB]")
ax.legend()
plt.show()

Relación con el resto de la biblioteca

El espectrograma comparte ventana, segmentación y escalado con power_spectral_density, así que ambos son mutuamente congruentes bin a bin; el espectrograma en bandas de octava de OctaveFilterBank.spectrogram es el equivalente de ancho de banda porcentual constante con balística de sonómetro; y para seguir en el tiempo la frecuencia de una sola componente, la frecuencia instantánea de Hilbert de envelope complementa la cresta de la STFT.

Qué cubre esta guía

Cubierto
El espectrograma calibrado de Bendat y Piersol (sección 12.6.4.2): segmentos STFT enventanados y solapados conservados como columnas de tiempo separadas en lugar de promediados, implementado por spectrogram con el escalado exacto del módulo de Welch y el error aleatorio de la Ec. 8.158 de una celda sin promediar. La FFT con zoom de la sección 11.5.4 (ecs. 11.122-11.130): espectros de banda limitada en una malla arbitrariamente fina, calculados como el equivalente chirp-Z de la cadena demodulación-diezmado-DFT del libro, implementada por zoom_fft. Ambas comparten el catálogo de ventanas de Harris (1978) tras el enventanado del segmento.
No cubierto
Esta es solo la familia de ancho de banda constante basada en STFT. El equivalente de porcentaje constante con balística de sonómetro es OctaveFilterBank.spectrogram, en la página de Niveles. Seguir en el tiempo la frecuencia instantánea de una sola componente es la envolvente de Hilbert de envelope, en la página de Correlación y retardo, no en esta. Como la página de análisis espectral, esta implementa los estimadores de Bendat y Piersol, no una norma de certificación.

Véase también

Análisis espectral: la estimación de Welch promediada para el fondo estacionario, y las figuras de mérito de las ventanas tras el enventanado del segmento.
Niveles: el espectrograma en bandas de fracción de octava con balística de sonómetro.
Correlación y retardo: la frecuencia instantánea de Hilbert para seguir una componente.
Referencia de la API: signals.time_frequency.
Teoría: Resolución frecuencial vs separación de bins FFT: la separación de bins frente a la resolución efectiva, el compromiso al que esta página dedica toda su extensión.

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 12.6.4.2 (espectrogramas y sus errores aleatorios), sección 11.5.4 (procedimientos de transformada con zoom, ecs. 11.122-11.130) y secciones 8.5.1/8.5.4 (ancho de banda de resolución y errores estadísticos de estimaciones sin promediar). ISBN 978-0-470-24877-5.
Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings of the IEEE, 66(1), 51-83. https://doi.org/10.1109/PROC.1978.10837El catálogo de ventanas tras el enventanado del segmento: anchura del lóbulo principal frente a fuga de lóbulos laterales, el compromiso que fija qué puede resolver una columna del espectrograma a una longitud de segmento dada.

Creado y mantenido por José Manuel Requena Plens ORCID Todos los proyectos