Análisis espectral calibrado

Referencias: Bendat y Piersol 2010 Welch 1967 Harris 1978 Thomson 1982 Percival y Walden 1993

Un espectro sin su incertidumbre es media medición. Esta página cubre los estimadores espectrales de Welch de phonometry.signals que informan, junto al propio espectro, de la calidad estadística de la estimación siguiendo a Bendat y Piersol, Random Data: Analysis and Measurement Procedures (4.ª ed., 2010): la densidad espectral de potencia y la densidad espectral cruzada con el número efectivo de promedios, el error aleatorio normalizado y los intervalos de confianza chi-cuadrado; el espectro de salida coherente que separa una salida medida en la parte explicada linealmente por la entrada y el resto de ruido, con la relación señal-ruido espectral; un suavizador en fracciones de octava con núcleo de potencia constante; y generadores de ruido de colores con pendiente exacta en ley de potencias para ejercitar todo lo anterior. Un estimador multitaper de Thomson (Percival y Walden, 1993) completa la familia para registros demasiado cortos para segmentar. Cada fórmula de error es una forma cerrada de las fuentes, verificada por Monte Carlo con semilla en la batería de tests.

1. Densidad espectral de potencia con su error estadístico

power_spectral_density estima la densidad autoespectral unilateral $G_{xx} (f)$ por el método de Welch: el registro se divide en segmentos con ventana (Hann por defecto) y solape del 50 % cuyos periodogramas se promedian. No se aplica eliminación de tendencia, así que la calibración absoluta se conserva: una señal en pascales da Pa²/Hz. Hay dos escalados: 'density' (unidades²/Hz, integra a la potencia de la señal) y 'spectrum' (unidades², lee directamente la potencia de tonos discretos).

El precio de conservar la calibración es que cualquier desplazamiento de continua o deriva lenta se queda en el registro y se filtra. El desplazamiento cae en el bin de continua, y los lóbulos laterales de la ventana de Hann reparten una fracción de él por los primeros bins, así que un desplazamiento sin corregir aparece como una subida espuria en baja frecuencia que no elimina ninguna cantidad de promediado. Cuando importen los bins más bajos, resta la media o filtra el registro en paso alto explícitamente antes de la llamada (eso es un cambio en la señal, no un ajuste del estimador), y ten en cuenta que ese mismo desplazamiento es invisible en una representación en fracciones de octava, porque ninguna banda llega a continua. La misma condición reaparece en Cualificación de datos, cuyos estadísticos de cruce de nivel y de picos están escritos para un proceso de media nula.

El estimador es una cadena fija: la ventana fija el ancho de banda de resolución, y los promedios efectivos se derivan de la longitud del registro (el número de segmentos brutos) junto con la ventana y el solape. Con ellos fijados, cada cifra de calidad se deriva de la decisión principal de diseño, la longitud del segmento. El diagrama la recorre con los números del ejemplo de esta página.

Diagrama de bloques de la cadena de la PSD de Welch: un registro de ruido rosa de 20 segundos a 48 kilohercios se divide en segmentos de 4096 muestras con 50 por ciento de solape dando 467 segmentos y una separación de bins de 11,7 hercios, cada segmento se enventana con Hann para un ancho de banda de resolución de 17,6 hercios, los periodogramas unilaterales del cuadrado de la FFT se promedian en 442 promedios efectivos, y el resultado es Gxx de f con intervalo de confianza chi-cuadrado, un error aleatorio de 1 entre la raíz de n d igual al 4,8 por ciento y unos 885 grados de libertad; una nota final enuncia el compromiso de que los segmentos largos compran resolución pero gastan promedios

Promediar $n_{d}$ segmentos independientes da al estimador $2 n_{d}$ grados de libertad chi-cuadrado (ec. 8.162), de donde sale todo lo demás:

ε_{r} [\hat{G}_{xx}] = \frac{1}{n _{d}}, \frac{n G ^ _{xx}}{χ _{n; α /2}^{2}} \leq G_{xx} \leq \frac{n G ^ _{xx}}{χ _{n; 1 - α /2}^{2}}, n = 2 n_{d} .

Con segmentos solapados y con ventana los promedios están correlacionados, así que el resultado informa tanto del número bruto de segmentos (n_segments) como del número efectivo de promedios independientes (n_averages), calculado con la fórmula de correlación de ventana de Welch (1967) que Bendat y Piersol citan en la sección 11.5.2.2: para Hann con solape del 50 %, aproximadamente 0,95 del número bruto. El error aleatorio y el intervalo de confianza usan el valor efectivo. En DC, y en Nyquist con longitud de segmento par, el espectro unilateral tiene una sola componente de Fourier real, así que esos bins llevan la mitad de grados de libertad y un intervalo proporcionalmente más ancho.

from phonometry import power_spectral_density

# record: un canal de una captura calibrada, en pascales (ver la guía de
#   Calibración más abajo); trabajando en dBFS es el vector bruto en [-1, 1] y
#   el resultado sale en FS^2/Hz. fs: su frecuencia de muestreo en Hz.
res = power_spectral_density(record, fs)          # Hann, solape 50 %, IC 95 %
print(res.n_averages, res.random_error)           # nd y 1/sqrt(nd)
print(res.ci_lower[10], res.psd[10], res.ci_upper[10])
res.plot()                                        # PSD en dB con la banda de IC

El sesgo de resolución es la otra mitad del presupuesto de error: un ancho de banda de análisis finito $B_{e}$ (expuesto como resolution_bandwidth, el ancho de banda de ruido efectivo de la ventana) suaviza los rasgos espectrales abruptos, siempre en la dirección de reducir el rango dinámico (ec. 8.139). Para un pico resonante de ancho de banda de media potencia $B_{r}$ , el sesgo normalizado de primer orden es la forma cerrada de la ec. 8.141, expuesta como resolution_bias_error:

ε_{b} [\hat{G}_{xx} (f_{r})] \approx - \frac{1}{3} (\frac{B _{e}}{B _{r}})^{2} .

from phonometry import resolution_bias_error

eps_b = resolution_bias_error(res.resolution_bandwidth, 25.0)  # pico de Br = 25 Hz

Un $B_{e}$ estrecho (segmentos largos) suprime el sesgo pero deja menos promedios y un error aleatorio mayor; los dos requisitos sobre la longitud de segmento tiran en direcciones opuestas, y ese es exactamente el compromiso que los números expuestos hacen visible.

Elegir la longitud del segmento, y después la del registro

La fórmula del sesgo es más útil leída al revés. Mantener el sesgo del pico de una resonancia por debajo de 1 dB significa $ε_{b} \geq - 0, 2$ , y por tanto $B_{e} \leq 0, 775 B_{r}$ ; mantenerlo por debajo de una décima de decibelio significa $B_{e} \leq 0, 263 B_{r}$ . La regla empírica de siempre, tres o cuatro anchos de banda de análisis a lo ancho de la anchura a media potencia del rasgo más estrecho que importe, es exactamente esta fórmula, y ahora lleva un número al lado.

Trabaja hacia delante desde la medición y no desde un nperseg de costumbre. Un modo a 100 Hz con $Q = 50$ tiene $B_{r} = f / Q = 2$ Hz. Mantener su sesgo por debajo del 2 % (0,09 dB) pide $B_{e} \leq B_{r} 3 \times 0, 02 = 0, 49$ Hz y, para la ventana de Hann, $B_{e} = 1, 5 f_{s} / nperseg$ , así que a 48 kHz eso es nperseg $\geq 147000$ , es decir, segmentos de al menos 3,1 s. Pedir después un error aleatorio del 10 % necesita $n_{d} = 100$ promedios efectivos, lo que con un solape del 50 % significa un registro de unos 150 s:

target_bias = 0.02                     # 2 % sobre el pico, unos 0,09 dB
b_r, q = 100.0 / 50.0, 100             # un modo de 100 Hz con Q = 50
b_e = b_r * (3 * target_bias) ** 0.5   # el ancho de banda de análisis que permite
nperseg = 1.5 * 48000 / b_e            # Hann: Be = 1,5 fs / nperseg
record_s = q * 0.5 * nperseg / 48000   # nd = 100 promedios con solape del 50 %
print(round(b_e, 2), int(nperseg), round(record_s))     # 0.49 146969 153

Izquierda: la estimación de Welch de una resonancia de 25 hercios de ancho en 1 kilohercio calculada con longitudes de segmento de 512, 1024, 4096 y 16384 muestras; los segmentos cortos recortan y ensanchan el pico mientras que los largos lo resuelven. Derecha: el sesgo de resolución y el error aleatorio en decibelios frente a la longitud de segmento en eje logarítmico, con el sesgo cayendo abruptamente y el error aleatorio subiendo despacio, cruzándose cerca de 7400 muestras, y el déficit medido del pico siguiendo al sesgo

La misma resonancia de 25 Hz de ancho vista con cuatro longitudes de segmento. Con nperseg = 512 el ancho de banda de análisis es de 140,6 Hz, cinco veces la anchura de la resonancia, y el pico se lee unos 8 dB de menos y mucho más ancho de la cuenta; con 4096 el ancho de banda es de 17,6 Hz y el pico está esencialmente resuelto. El panel derecho pone los dos errores en un mismo eje: el sesgo cae como $B_{e}^{2}$ y el error aleatorio sube como $n_{d}^{- 1/2}$ , y se igualan cerca de nperseg = 7400 para esta resonancia, que es la longitud de segmento a elegir si no hay razón para preferir un error sobre el otro.

Las dos decisiones son independientes: la longitud del segmento la fija el rasgo más agudo que importe, y la del registro la precisión que quieras. Alargar los segmentos sin alargar el registro cambia un error por el otro y no compra nada. Y 150 s de registro son 150 s durante los cuales el proceso tiene que mantenerse estacionario, lo cual es una afirmación y no una suposición, y es la que decide Cualificación de datos. Cuando el registro sencillamente no puede ser tan largo, la salida es el estimador multitaper del §7.

De unidades digitales a dB SPL

power_spectral_density escala lo que se le dé y no aplica ninguna calibración propia, así que el factor de sensibilidad va sobre el registro, antes de la estimación:

calibration_factor = 1.002             # Pa por unidad digital, del calibrador
pressure = calibration_factor * record
psd = power_spectral_density(pressure, fs)      # ahora en Pa^2/Hz

La ordenada es entonces $10 lo g_{10} (G_{xx} / p_{0}^{2})$ dB re (20 µPa)²/Hz, y con scaling="spectrum" esa misma llamada lee directamente el nivel de un tono discreto como $10 lo g_{10} (P / p_{0}^{2})$ dB. Para pasar de ahí al nivel de banda que contiene un informe acústico, integra la densidad en la banda (suma los bins y multiplica por la separación entre bins) y divide por $p_{0}^{2}$ ; fractional_octave_smoothing del §4 es un suavizador de representación y no un sustituto de esa integración, como detalla ese mismo apartado. El factor sale de Calibración y dBFS.

Densidad espectral de potencia de Welch de ruido rosa en dB por Hz entre 20 Hz y 20 kHz, con la banda de confianza chi-cuadrado del 95 por ciento sombreada alrededor de la estimación, la curva suavizada en 1/3 de octava encima y la ley de potencias exacta de -3,01 dB por octava como línea discontinua de referencia

Lo que compran 20 s de registro, leído sobre la figura. 467 segmentos de 4096 muestras dan 442 promedios efectivos, un error aleatorio normalizado de 0,047 y una banda de confianza chi-cuadrado del 95 % de 0,81 dB de ancho, así que la estimación dentada no se aparta de la curva suavizada más de unos 0,8 dB en toda la década, y la afirmación honesta de lo que se sabe es la banda sombreada, no el rizado. Las dos curvas ajustan −3,007 dB/octava frente a los −3,01 dB/octava exactos del ruido rosa. Suavizar no estrecha nada: redibuja la misma estimación, y el intervalo que hay debajo sigue valiendo.

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from phonometry import (
    fractional_octave_smoothing,
    noise_signal,
    power_spectral_density,
)

fs = 48000.0
x = noise_signal(fs, 20.0, color="pink", seed=11)
res = power_spectral_density(x, fs, nperseg=4096)
band = (res.frequencies >= 20.0) & (res.frequencies <= 20000.0)
freqs = res.frequencies[band]
smooth = fractional_octave_smoothing(res.frequencies, res.psd, 3.0)[band]

fig, ax = plt.subplots(figsize=(10, 6))
ax.fill_between(freqs, 10 * np.log10(res.ci_lower[band]),
                10 * np.log10(res.ci_upper[band]), alpha=0.3,
                label="Intervalo de confianza chi-cuadrado del 95 %")
ax.semilogx(freqs, 10 * np.log10(res.psd[band]), lw=1.0,
            label="Estimación de la PSD de Welch")
ax.semilogx(freqs, 10 * np.log10(smooth), lw=2.2,
            label="Suavizado en 1/3 de octava")
ax.set_xlabel("Frecuencia [Hz]")
ax.set_ylabel("PSD [dB re 1/Hz]")
ax.legend()
plt.show()

2. Densidad espectral cruzada

cross_spectral_density estima la $G_{x y} (f)$ compleja entre dos canales con el mismo núcleo de Welch, e informa de la coherencia ordinaria $γ_{x y}^{2} = ∣ G_{x y} ∣^{2} / (G_{xx} G_{y y})$ junto a los errores aleatorios de Bendat y Piersol de la magnitud y la fase (ecs. 9.33 y 9.52, con la coherencia medida en lugar del valor verdadero desconocido, como recomienda el libro para datos medidos):

ε_{r} [∣ \hat{G}_{x y} ∣] = \frac{1}{∣ γ _{x y} ∣ n _{d}}, s.d. [\hat{θ}_{x y}] = \frac{( 1 - γ _{x y}^{2} ) ^{1/2}}{∣ γ _{x y} ∣ 2 n _{d}} .

Ambos se reducen cuando la coherencia se acerca a uno: un par fuertemente coherente necesita muchos menos promedios para la misma confianza. La fase se devuelve desenrollada, así que su pendiente frente a la frecuencia es el retardo de grupo $τ_{g} = - d φ / (2 π d f)$ ; para un camino de retardo puro la fase es lineal y esa pendiente lee directamente el retardo de propagación.

from phonometry import cross_spectral_density

res = cross_spectral_density(x, y, fs)
print(res.magnitude_random_error[100], res.phase_std[100])  # errores del bin 100
res.plot()   # magnitud, fase con banda ±sigma, coherencia

Dos paneles para la densidad espectral cruzada de un camino de dos sensores con un retardo de 2 milisegundos: la estimación de magnitud de Welch fluctuando alrededor de un nivel plano, y debajo la fase desenrollada del espectro cruzado cayendo como una línea recta sobre un eje de frecuencia logarítmico, exactamente sobre la referencia discontinua de menos dos pi f tau con una banda estrecha de una sigma a su alrededor

La densidad espectral cruzada de un camino con retardo de 2 ms: la fase desenrollada es exactamente la recta $- 2 π f τ$ , así que su pendiente lee el retardo de propagación directamente, y la banda ±1 d.e. de la Ec. 9.52 cuantifica cuánto fiarse de ella por frecuencia.

Mostrar el código de esta figura

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

fs = 8000.0
tau = 0.002                                   # 2 ms = 16 muestras
delay = int(tau * fs)
x = noise_signal(fs, 8.0, seed=8)
noise = noise_signal(fs, 8.0, rms=0.3, seed=9)
y = 0.9 * np.concatenate([np.zeros(delay), x[:-delay]]) + noise

res = cross_spectral_density(x, y, fs)

# Una línea: magnitud, fase con su banda ±sigma y coherencia:
res.plot(language="es")
plt.show()

# A mano, desde los campos que lleva el resultado:
band = (res.frequencies >= 20) & (res.frequencies <= 3500)
freqs = res.frequencies[band]
fig, (ax_m, ax_p) = plt.subplots(2, 1, sharex=True)
ax_m.semilogx(freqs, 10 * np.log10(res.magnitude[band]),
              label="|Gxy| (estimación de Welch)")
ax_m.set_ylabel("Magnitud [dB]")
ax_p.semilogx(freqs, res.phase[band], label="Fase desenrollada")
ax_p.fill_between(freqs, res.phase[band] - res.phase_std[band],
                  res.phase[band] + res.phase_std[band], alpha=0.25,
                  label="±1 d.e. (Ec. 9.52)")
ax_p.semilogx(freqs, -2 * np.pi * freqs * tau, "r--",
              label="pendiente -2·pi·f·tau")
ax_p.set(xlabel="Frecuencia [Hz]", ylabel="Fase [rad]")
for ax in (ax_m, ax_p):
    ax.legend()
plt.show()

3. Espectro de salida coherente y SNR espectral

En el modelo de una entrada y una salida, el autoespectro medido de la salida se separa exactamente en la parte explicada linealmente por la entrada y el resto no correlacionado (ecs. 9.55–9.57):

G_{v v} = γ_{x y}^{2} G_{y y}, G_{nn} = (1 - γ_{x y}^{2}) G_{y y}, SNR (f) = \frac{γ _{x y}^{2}}{1 - γ _{x y}^{2}} .

coherent_output_spectrum devuelve los tres espectros, la relación señal-ruido espectral (lineal y en dB) y el error aleatorio del estimador de la salida coherente (ec. 9.73), más la propagación de primer orden del error de la coherencia a través de la SNR:

ε_{r} [\hat{G}_{v v}] = \frac{( 2 - γ _{x y}^{2} ) ^{1/2}}{∣ γ _{x y} ∣ n _{d}}, ε_{r} [SNR] = \frac{2}{∣ γ _{x y} ∣ n _{d}} .

Para ruido aditivo no correlacionado en la salida con nivel conocido, la coherencia tiene la forma cerrada $γ^{2} = SNR / (1 + SNR)$ , lo que hace toda la cadena verificable con una señal sintética:

import numpy as np
from phonometry import coherent_output_spectrum, noise_signal

fs = 48000.0
x = noise_signal(fs, 8.0, color="white", seed=1)
noise = noise_signal(fs, 8.0, color="white", rms=0.5, seed=2)
y = 0.8 * x + noise                      # SNR = 0.64/0.25 en toda frecuencia

res = coherent_output_spectrum(x, y, fs)
print(np.median(res.coherence))          # -> SNR/(1+SNR) = 0.719
print(np.median(res.snr_db))             # -> 10·lg(2.56) = 4.1 dB
res.plot()                               # Gyy, Gvv, Gnn y el panel de SNR

Dos paneles para el espectro de salida coherente de un sistema de ruido blanco con ruido aditivo: el espectro de salida medido con la parte coherente aproximadamente 1,5 decibelios por debajo y el suelo de ruido no correlacionado unos 6 decibelios más abajo, y debajo la relación señal-ruido espectral fluctuando alrededor de la línea discontinua de forma cerrada en 4,1 decibelios

El reparto exacto del modelo del fragmento: $G_{v v} = γ^{2} G_{y y}$ explica la salida salvo el resto de ruido plano $G_{nn}$ , y la SNR espectral se dispersa alrededor de su forma cerrada $10 lo g_{10} (0, 64/0, 25) = 4, 1$ dB en toda frecuencia.

Mostrar el código de esta figura

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

fs = 48000.0
x = noise_signal(fs, 8.0, color="white", seed=1)
noise = noise_signal(fs, 8.0, color="white", rms=0.5, seed=2)
y = 0.8 * x + noise                      # SNR = 0.64/0.25 por banda

res = coherent_output_spectrum(x, y, fs, nperseg=2048)

# Una línea: los tres espectros y el panel de SNR:
res.plot(language="es")
plt.show()

# A mano, desde los campos que lleva el resultado:
band = (res.frequencies >= 20) & (res.frequencies <= 20000)
freqs = res.frequencies[band]
fig, (ax_g, ax_s) = plt.subplots(2, 1, sharex=True)
for values, style, label in ((res.output_psd, "-", "Gyy (medida)"),
                             (res.coherent_psd, "--", "Gvv (coherente)"),
                             (res.noise_psd, ":", "Gnn (ruido)")):
    ax_g.semilogx(freqs, 10 * np.log10(values[band]), style, label=label)
ax_g.set_ylabel("Densidad espectral [dB re 1/Hz]")
ax_s.semilogx(freqs, res.snr_db[band], label="SNR espectral [dB]")
ax_s.axhline(10 * np.log10(0.64 / 0.25), color="r", linestyle="--",
             label="forma cerrada 4,1 dB")
ax_s.set(xlabel="Frecuencia [Hz]", ylabel="SNR [dB]")
for ax in (ax_g, ax_s):
    ax.legend()
plt.show()

El campo coherence_bias informa del pequeño sesgo positivo del estimador de coherencia, $b [\overset{γ}{^}^{2}] \approx (1 - γ^{2})^{2} / n_{d}$ (ec. 9.75): despreciable en cuanto $n_{d}$ llega a unos cientos, y otra razón para promediar con generosidad antes de fiarse de una coherencia baja.

Dos canales en la práctica

Los §2 y §3 parten los dos de una pareja de registros, y todo lo que informan vale lo que valga esa pareja. Dos mitades: cómo adquirirla y cómo leer una coherencia que vuelve por debajo de uno.

Adquisición. Usa una sola interfaz con entradas de muestreo simultáneo, o mide el desfase entre canales y elimínalo. Un convertidor que secuencia sus canales mete un desfase fijo directamente en la fase del espectro cruzado y, por tanto, en todos los retardos de grupo que se lean de ella: una muestra de desfase a 48 kHz son 20,8 µs, que son 15° de fase a 2 kHz y un error plano de 20,8 µs en cualquier estimación de retardo. Fija las dos ganancias antes de la tirada y no las toques, porque el reparto de salida coherente del §3 supone que los dos canales comparten una escala fija. Valida la pareja una vez con una tirada de base cero: mete la misma señal en los dos canales (el generador partido hacia ambas entradas, o los dos micrófonos uno al lado del otro) y confirma que el retardo medido es una fracción pequeña de muestra y que la coherencia vale uno en toda la banda. Esa única tirada separa un problema de cadena de un problema de física para todas las mediciones que vengan después.

Leer una coherencia por debajo de uno. No es automáticamente un problema de relación señal-ruido, y Bendat y Piersol nombran cinco causas con arreglos distintos:

Lo que ves	La causa	El arreglo
Coherencia plana en toda la banda, coherente con $γ^{2} = SNR / (1 + SNR)$	ruido no correlacionado en alguno de los sensores, el caso modelado del §3	promedia más, sube el nivel, acércate
Coherencia que cae progresivamente con la frecuencia, peor cuanto más corto se segmenta el registro	un retardo global de propagación comparable con la longitud del segmento	alinea antes los registros, o sube `nperseg` muy por encima del retardo; ver Correlación y retardo
Hundimientos abruptos justo en las resonancias	sesgo de resolución: el ancho de banda de análisis emborrona el pico de forma distinta en cada canal	alarga el segmento; la fórmula de $ε_{b}$ del §1 dice cuánto
Baja en los armónicos de un tono fuerte mientras la coherencia de banda ancha sigue alta	no linealidad, un traqueteo o un canal recortado	arregla el camino, o trabaja por debajo del nivel al que aparece
Baja en toda una banda sin estructura evidente	una segunda fuente no correlacionada está excitando también la salida	la página de MISO: una entrada no puede explicar dos fuentes

Dos advertencias de cierre. La estimación de coherencia está sesgada al alza en $(1 - γ^{2})^{2} / n_{d}$ , así que una coherencia baja leída con un puñado de promedios es todavía más baja de lo que parece. Y, leído al revés, una coherencia cercana a uno es el criterio de aceptación que permite declarar toda la estimación: decláralo banda a banda junto a la función de transferencia, porque las bandas donde cae son exactamente las bandas donde los números de encima no significan nada.

4. Suavizado en fracciones de octava

fractional_octave_smoothing promedia un espectro sobre una ventana rectangular de anchura relativa constante: 1/n de octava, $[f \cdot 2^{- 1/2 n}, f \cdot 2^{+ 1/2 n}]$ alrededor de cada frecuencia. Es el ancho de banda de resolución de porcentaje constante que Bendat y Piersol recomiendan para espectros de sistemas resonantes (sección 8.5.3), y el estándar de facto para presentar respuestas de altavoces y salas. El promedio se calcula siempre sobre potencia (las amplitudes se elevan al cuadrado primero, los niveles en dB se convierten ida y vuelta), así que se conserva la potencia de banda y no la amplitud, y un espectro plano pasa exactamente sin cambios.

Los tres valores de domain= nombran qué es la ordenada que se le entrega, porque el promedio se hace siempre sobre potencia y la conversión de vuelta tiene que saber de dónde partió: una densidad o un espectro de potencia ("power", el valor por defecto), una magnitud lineal como el $∣ H ∣$ de una respuesta en frecuencia ("amplitude", que se eleva antes al cuadrado), o una curva ya en decibelios ("db", que se convierte y se vuelve a convertir).

from phonometry import fractional_octave_smoothing

freqs = res.frequencies                    # de la estimación del §1
smooth_psd = fractional_octave_smoothing(freqs, res.psd, 3.0)
magnitude = np.sqrt(res.psd)               # aquí valdría cualquier |H| lineal
smooth_mag = fractional_octave_smoothing(freqs, magnitude, 6.0,
                                         domain="amplitude")  # una FRF |H|
levels = 10.0 * np.log10(res.psd)          # aquí valdría cualquier curva en dB
smooth_db = fractional_octave_smoothing(freqs, levels, 3.0, domain="db")

Una sola línea espectral con ordenada de PSD $P$ (unidades²/Hz) en un bin de anchura $Δ f$ se suaviza al nivel en forma cerrada $P Δ f / (f_{0} (2^{1/2 n} - 2^{- 1/2 n}))$ sobre una anchura de núcleo: el oráculo fijado en los tests.

Una densidad suavizada sigue siendo una densidad

El suavizado cambia la resolución de la estimación, nunca sus unidades. La salida de fractional_octave_smoothing es una densidad en unidades²/Hz, exactamente como lo era la entrada (el ruido blanco plano se suaviza a la misma densidad plana), mientras que los niveles de banda de tercio de octava de ese mismo ruido blanco suben 3 dB por banda. Superpón una PSD suavizada sobre un espectro de banco de filtros sin convertir y las dos discreparán en 20 o 30 dB, y ninguna estará equivocada.

La conversión es un solo término. Una banda de fracción de octava de frecuencia central $f_{c}$ tiene un ancho de banda $B_{n} = f_{c} (2^{1/2 n} - 2^{- 1/2 n})$ , es decir, $0, 2316 f_{c}$ en un tercio de octava y $0, 7071 f_{c}$ en una octava completa, así que

L_{banda} = 10 lo g_{10} (\frac{G ( f _{c} ) B _{n}}{p _{0}^{2}}),

que a 1 kHz en tercio de octava añade $10 lo g_{10} (0, 2316 \times 1000) = 23, 65$ dB al nivel de densidad, y 3 dB más por cada octava hacia arriba. Sumar los bins de la PSD a lo ancho de la banda y multiplicar por la separación entre bins es la forma exacta; la expresión de arriba es su aproximación de banda plana.

Qué representación va dónde: la densidad suavizada para comparar respuestas en magnitud de altavoces y salas, donde el ojo quiere una curva a resolución relativa constante; el nivel de banda para todo lo que se enfrente a un criterio, o que tenga que cuadrar con octave_filter y OctaveFilterBank; ver Niveles.

5. Generadores de ruido de colores

noise_signal produce ruido gaussiano cuya PSD sigue $G_{xx} (f) \propto f^{α}$ de forma exacta en esperanza: ruido blanco con semilla se moldea en el dominio de la frecuencia con la respuesta en magnitud exacta $(f / f_{ref})^{α /2}$ bin a bin (un filtro de fase cero aplicado de forma circular), así que una pendiente medida solo se desvía de la ley de potencias por el error aleatorio de la estimación espectral, no como las aproximaciones rosas por tramos o de pocos polos cuya pendiente ondula fracciones de dB. El registro tiene media cero y se reescala exactamente al RMS pedido, y la misma semilla reproduce el mismo registro bit a bit.

color	$α$	pendiente de la PSD
`white`	0	0 dB/octava
`pink`	-1	-3,01 dB/octava
`red` (browniano)	-2	-6,02 dB/octava
`blue`	+1	+3,01 dB/octava
`violet`	+2	+6,02 dB/octava

from phonometry import noise_signal

pink = noise_signal(48000, 10.0, color="pink", seed=7)     # determinista
white = noise_signal(48000, 10.0, color="white", rms=0.5, seed=7)

Medida sobre tres décadas (20 Hz – 20 kHz) con el estimador de la sección 1, la pendiente de regresión de cada color cae a unas milésimas de dB/octava del valor exacto: la batería de conformidad fija la pendiente rosa en -3,0116 frente al exacto -3,0103.

Esas pendientes son pendientes de densidad, y no son lo que muestra una representación en fracciones de octava de la misma señal. La anchura de una banda crece en proporción a su frecuencia central, así que un nivel de banda es la densidad más $10 lo g_{10} B_{n}$ con $B_{n} \propto f_{c}$ : todos los colores se leen 3 dB por octava más altos en una representación por bandas de lo que sugiere su pendiente de densidad. El ruido blanco, cuya densidad es plana, sube por tanto 3 dB por octava a través de un banco de filtros de octava, y el ruido rosa, a −3 dB/octava, sale plano, que es por lo que el rosa es el estímulo de referencia en trabajo de salas y altavoces, y por lo que un espectro en bandas de octava de ruido blanco que sube no es ningún defecto. El ruido rojo cae 3 dB por octava en una representación por bandas y el azul sube 6. Los mismos +3 dB por octava valen para cualquier comparación de densidad frente a banda en esta biblioteca; el §4 de arriba da la conversión exacta, y en Niveles vive el lado de las bandas.

Densidades espectrales de Welch de los cinco colores de ruido entre 20 hercios y 20 kilohercios, cada una normalizada a su propio nivel en 1 kilohercio de modo que las cinco rectas se abren en abanico desde un único punto: el violeta subiendo a 6 dB por octava, el azul a 3, el blanco plano, el rosa cayendo a 3 y el rojo a 6, cada uno con su ley de potencias exacta discontinua debajo y la pendiente de regresión medida en la leyenda

Los cinco generadores a lo largo de tres décadas, cada uno normalizado a su propio nivel en 1 kHz para que las pendientes se abran en abanico desde un solo punto, con la ley de potencias exacta discontinua debajo. Las pendientes de regresión medidas caen a menos de cuatro milésimas de decibelio por octava de los valores exactos, y el residuo es el error aleatorio del estimador, no del generador: la misma semilla reproduce el registro bit a bit.

6. Elegir la ventana

Todos los estimadores de esta página aceptan cualquier ventana que conozca scipy.signal.get_window, pero la elección es un compromiso cuantificado, no una preferencia. window_metrics calcula las figuras de mérito que Harris (1978) tabuló, para cualquier ventana y longitud, muestreadas en modo periódico (DFT-even) exactamente como las aplican los estimadores de Welch:

ENBW (ancho de banda equivalente de ruido, en bins): cuánto más ancho que un bin es el ancho de banda de análisis efectivo. Es el mismo número que el resultado de la PSD informa como resolution_bandwidth (ENBW·fs/nperseg en Hz), y entra directamente en el balance tono/ruido: un suelo de ruido de banda ancha leído de un espectro con ventana queda $10 lo g_{10} (ENBW)$ dB por encima de la densidad verdadera.
Ganancia coherente: la ganancia en continua $\sum w / N$ que escala un tono centrado en un bin.
Pérdida de festoneado: la atenuación en el peor caso de un tono que cae a medio camino entre dos bins (3,92 dB para la rectangular, 1,42 dB para Hann).
Pérdida de proceso en el peor caso: el festoneado más $10 lo g_{10} (ENBW)$ , la reducción de SNR de salida en el peor caso al detectar un tono en ruido blanco.
Lóbulo lateral máximo y anchura a -3 dB del lóbulo principal: suelo de fuga frente a resolución.

from phonometry import window_metrics

m = window_metrics("hann", 2048)
print(m.enbw_bins)              # 1.5, exacto
print(m.scalloping_loss_db)     # 1.42 dB
print(m.highest_sidelobe_db)    # -31.5 dB
m.plot()                        # ventana + espectro con las métricas

Las formas cerradas anclan los tests: el ENBW es exactamente 1 para la rectangular, 3/2 para Hann, 1987/1458 para Hamming y 1523/882 para Blackman (muestreo DFT-even), y la pérdida de festoneado de la rectangular es $20 lo g_{10} (N sin (π /2 N))$ , el núcleo de Dirichlet evaluado medio bin fuera del centro.

Espectros de las ventanas rectangular, Hann, Hamming y Blackman a lo largo de 16 bins de la DFT, mostrando el compromiso entre anchura del lóbulo principal y nivel de los lóbulos laterales, con el ancho de banda equivalente de ruido y el lóbulo lateral máximo de cada ventana en la leyenda

El compromiso en una imagen, y 0,73 de bin cuestan 44,8 dB de suelo de fuga. Pasar de la rectangular a la Blackman baja el lóbulo lateral máximo de −13,3 dB a −58,1 dB mientras el ENBW se ensancha de 1,000 a 1,727 bins y el lóbulo principal a −3 dB de 0,886 a 1,644 bins. La Hann se queda donde está el valor por defecto del módulo: −31,5 dB por 1,500 bins. Lee los dos ejes juntos: un lóbulo principal estrecho separa dos tonos de nivel parecido, un suelo de lóbulos laterales bajo encuentra un tono débil al lado de uno fuerte, y ninguna ventana hace las dos cosas.

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from phonometry import window_metrics

n, oversample = 1024, 256
fig, ax = plt.subplots(figsize=(10, 6.2))
for name in ("boxcar", "hann", "hamming", "blackman"):
    res = window_metrics(name, n)
    spectrum = np.abs(np.fft.rfft(res.taps, n=n * oversample))
    level = 20.0 * np.log10(spectrum / spectrum[0])
    bins = np.arange(level.size) / oversample
    shown = bins <= 16.0
    ax.plot(bins[shown], level[shown],
            label=(f"{name}: ENBW {res.enbw_bins:.2f} bins, "
                   f"lóbulo lateral {res.highest_sidelobe_db:.1f} dB"))
ax.set_xlim(0.0, 16.0)
ax.set_ylim(-100.0, 5.0)
ax.set_xlabel("Desplazamiento en frecuencia [bins de la DFT]")
ax.set_ylabel("Nivel re lóbulo principal [dB]")
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()

La ventana Hann por defecto de este módulo es la elección equilibrada: la caída rápida de los lóbulos laterales (-18 dB/octava) protege los espectros de ruido frente a la fuga, su ENBW de 1,5 solo cuesta 1,76 dB frente a la rectangular, y su correlación de solape al 50 % conserva casi toda la información de los segmentos en el número efectivo de promedios. Recurre a una ventana de lóbulos más bajos (Blackman, Kaiser con beta alta) cuando haya que encontrar un tono débil junto a uno fuerte, aceptando el lóbulo principal más ancho; y a la rectangular solo para registros que se autoenventanan (transitorios que decaen dentro del segmento) o síntesis centrada en bins.

7. Estimación multitaper para registros cortos

El método de Welch compra estabilidad con longitud de registro: cada segmento independiente añade dos grados de libertad, así que un registro en el que solo caben un par de segmentos deja una estimación apenas mejor que un periodograma. multitaper_psd implementa la alternativa de Thomson (1982), tal y como la desarrollan Percival y Walden (1993, capítulo 7): el registro entero se multiplica por $K$ ventanas ortogonales esferoidales prolatas discretas (de Slepian) - las secuencias que concentran la mayor energía de su ventana espectral dentro de una banda de diseño $[- W, W]$ elegida - y se promedian los $K$ autoespectros propios resultantes:

\hat{S}^{(m t)} (f) = \frac{1}{K} k = 0 \sum K - 1 \hat{S}_{k} (f), \hat{S}_{k} (f) = Δ t t = 1 \sum N h_{t, k} x_{t} e^{- i 2 π f t Δ t}^{2} .

Como las ventanas son ortogonales, los autoespectros están casi no correlacionados, así que el promedio acarrea unos $2 K$ grados de libertad chi-cuadrado de un único registro - la misma maquinaria estadística que el resultado de Welch (error aleatorio, intervalo de confianza chi-cuadrado), sin segmentar. El parámetro de diseño es el producto adimensional de duración por semiancho de banda $p = N W$ (4 por defecto), y para un registro de $N_{r}$ muestras a $f_{s}$ fija el semiancho $W = p f_{s} / N_{r}$ en hercios. $2 W$ es la resolución de la estimación (informada como resolution_bandwidth), y solo las ventanas por debajo del número de Shannon $2 p$ mantienen la energía de su ventana espectral dentro de la banda de diseño - sus concentraciones $λ_{k}$ se informan como eigenvalues, y el número de ventanas por defecto es $K = 2 p - 1$ , todas las de concentración casi unidad. Un $p$ mayor admite más ventanas (menos varianza) a costa de resolución.

from phonometry import multitaper_psd

res = multitaper_psd(record, fs)                 # p = NW = 4, K = 7, adaptativo
print(res.degrees_of_freedom.mean())             # ~2K de un solo registro
print(res.eigenvalues)                           # concentraciones
res.plot(language="es")                          # densidad con la banda de IC

Por defecto los autoespectros se combinan con los pesos adaptativos de Thomson (P&W Ecs. 368a/370a, iterados hasta converger). El peso de cada ventana en cada frecuencia equilibra el espectro local frente a la fuga de banda ancha que esa ventana podría acarrear:

b_{k} (f) = \frac{S ( f )}{λ _{k} S ( f ) + ( 1 - λ _{k} ) σ ^{2} Δ t}, \hat{S}^{(am t)} (f) = \frac{\sum _{k} b _{k}^{2} ( f ) λ _{k} S ^ _{k} ( f )}{\sum _{k} b _{k}^{2} ( f ) λ _{k}},

de modo que las ventanas de orden alto, con más fuga, pierden peso exactamente donde el espectro es localmente débil, y no se pierde nada donde es localmente blanco (para ruido blanco los pesos son uniformes). El precio queda contabilizado con honestidad: los grados de libertad equivalentes pasan a depender de la frecuencia, $ν (f) = 2 (\sum_{k} d_{k})^{2} / \sum_{k} d_{k}^{2}$ con $d_{k} = b_{k}^{2} λ_{k}$ (P&W Ec. 370b), y el intervalo de confianza se ensancha allí donde la protección contra la fuga los gastó. adaptive=False selecciona en su lugar el promedio ponderado por autovalores.

La calibración coincide exactamente con los estimadores de Welch: sin eliminación de tendencia, 'density' integra a la potencia de la señal y 'spectrum' lee $A^{2} /2$ en el pico de una sinusoide de amplitud $A$ (la potencia de un tono en escalado 'density' se reparte sobre la banda $2 W$ ). Las propias ventanas de Slepian proceden de scipy.signal.windows.dpss; sus concentraciones reproducen la tabla en cuádruple precisión de Percival y Walden (tabla 382) hasta la precisión de máquina, que es el oráculo ancla de la batería de tests.

Izquierda: la densidad espectral multitaper de un registro de ruido rosa de 171 milisegundos con su banda de confianza del 95 por ciento, frente a una estimación de Welch con solo 6,7 promedios efectivos y una banda visiblemente más ancha, una estimación de una sola ventana como línea gris dentada y la ley exacta de menos 3,01 dB por octava discontinua. Derecha: un tono de 60 dB sobre un suelo rosa, donde la estimación de Welch con ventana de Hann deja una falda más ancha alrededor del tono que la multitaper adaptativa, y los grados de libertad equivalentes bajan de unos nueve a cinco allí donde los pesos adaptativos los gastan

A la izquierda, el caso en el que el texto sostiene que gana la multitaper: un registro de 171 ms. Welch con un segmento de 2048 muestras solo mete $n_{d} = 6, 7$ promedios efectivos dentro de él, así que su banda de confianza es visiblemente la más ancha de las dos y su curva la más rugosa, mientras que la estimación adaptativa de 7 ventanas acarrea 13,7 grados de libertad equivalentes del mismo registro. A la derecha, la otra razón para echar mano de ella: un tono de 60 dB sobre un suelo rosa. La ventana de Hann deja una falda más ancha alrededor del tono que los pesos adaptativos, y el precio está impreso en el eje derecho: los grados de libertad equivalentes caen allí donde los pesos compraron esa protección, que es exactamente lo que ensancha ahí el intervalo de confianza.

Mostrar el código de esta figura

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

fs = 48000.0
x = noise_signal(fs, 8192 / fs, color="pink", seed=11)   # registro de 171 ms
single = multitaper_psd(x, fs, n_tapers=1, adaptive=False)
res = multitaper_psd(x, fs)                              # NW = 4, K = 7
band = (res.frequencies >= 20.0) & (res.frequencies <= 20000.0)
freqs = res.frequencies[band]

fig, ax = plt.subplots(figsize=(10, 6))
ax.semilogx(freqs, 10 * np.log10(single.psd[band]), color="gray",
            alpha=0.45, lw=0.7, label="Una sola ventana de Slepian (K = 1)")
ax.fill_between(freqs, 10 * np.log10(res.ci_lower[band]),
                10 * np.log10(res.ci_upper[band]), alpha=0.3,
                label="Intervalo de confianza chi-cuadrado del 95 %")
ax.semilogx(freqs, 10 * np.log10(res.psd[band]), lw=1.2,
            label="Estimación multitaper (K = 7, adaptativa)")
ax.set_xlabel("Frecuencia [Hz]")
ax.set_ylabel("PSD [dB re 1/Hz]")
ax.legend()
plt.show()

Recurre a multitaper_psd cuando el registro sea demasiado corto para segmentar (colas de respuestas al impulso de salas, capturas de transitorios, ciclos sueltos de máquina) o cuando un espectro de gran rango dinámico necesite una protección contra la fuga que un promedio de Welch con ventana Hann no puede dar; quédate con power_spectral_density para registros largos, donde promediar segmentos es más barato que $K$ FFT de longitud completa y ambos estimadores coinciden.

Congruencia con los estimadores de respuesta en frecuencia y de intensidad

Los estimadores de respuesta en frecuencia forman los mismos espectros cruzados en una función de transferencia ( $H_{1} = G_{x y} / G_{xx}$ , insesgado cuando el ruido está en la salida, y $H_{2} = G_{y y} / G_{y x}$ , insesgado cuando está en la entrada), y la sonda de intensidad acústica de dos micrófonos forma la parte imaginaria de $G_{x y}$ . Esos estimadores, transfer_function y coherence, la sonda de intensidad y estos estimadores comparten un único núcleo de Welch (misma ventana, misma política de solape y calibración sin eliminación de tendencia), así que una PSD, una coherencia y una $H_{1}$ calculadas con la misma longitud de segmento son mutuamente congruentes bin a bin. La misma matriz de espectros cruzados sustenta la coherencia múltiple y parcial, que extiende la coherencia ordinaria a varias entradas correlacionadas y una salida.

Qué cubre esta guía

Cubierto
Los estimadores de Welch de Bendat y Piersol (Random Data, 4.ª ed., 2010, secciones 5.2, 8.5, 9.1-9.2 y 11.5): power_spectral_density y cross_spectral_density con el número efectivo de promedios, los intervalos de confianza chi-cuadrado y el error de sesgo por resolución; el espectro de salida coherente y la relación señal-ruido espectral; las figuras de mérito de las ventanas de Harris (1978) de window_metrics; y el estimador multitaper de Thomson (1982) de multitaper_psd, siguiendo los capítulos 7 y 8 de Percival y Walden. fractional_octave_smoothing y noise_signal completan la caja de herramientas.
No cubierto
Los estimadores de respuesta en frecuencia transfer_function y coherence, y la sonda de intensidad acústica, comparten el núcleo de Welch de esta página pero se documentan en las páginas de electroacústica y intensidad acústica, no aquí. Lo mismo ocurre con la coherencia múltiple y parcial, en la página de coherencia MISO. Esta página implementa los estimadores del libro de Bendat y Piersol, no una norma de certificación, así que no lleva números de apartado ni límites de cumplimiento que comprobar.

Véase también

Cualificación de datos: la estacionariedad que da por supuesta todo intervalo chi-cuadrado de esta página.
Calibración y dBFS: de dónde sale el factor que convierte las densidades de esta página en Pa²/Hz.
Coherencia MISO: la misma matriz de espectros cruzados con varias entradas correlacionadas.
Análisis tiempo-frecuencia: la misma segmentación de Welch representada en lugar de promediada.
Niveles: el lado de los niveles de banda de la conversión densidad/banda del §4.
Referencia de la API: signals.spectra, signals.windows y signals.multitaper.
Teoría: Resolución frecuencial vs separación de bins FFT: por qué la separación de bins de una FFT no es la resolución de la estimación, y qué le hace la ventana a ambas.

Referencias

Bendat, J. S. y Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4.ª ed.). Wiley. https://doi.org/10.1002/9781118032428Secciones 5.2 y 8.5 (autoespectros y sus errores aleatorios y de sesgo, intervalos chi-cuadrado), 9.1-9.2 (espectros cruzados, espectro de salida coherente y sus errores) y 11.5 (procesado de Welch, ventanas y solape). 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.10837Las figuras de mérito de las ventanas (tabla 1): ancho de banda equivalente de ruido, ganancia coherente, pérdida de festoneado, pérdida de proceso en el peor caso, nivel del lóbulo lateral máximo y anchura del lóbulo principal, calculadas por window_metrics para cualquier ventana de scipy.
Percival, D. B. y Walden, A. T. (1993). Spectral analysis for physical applications: Multitaper and conventional univariate techniques. Cambridge University Press. https://doi.org/10.1017/CBO9780511622762Capítulo 7 (estimación multitaper: autoespectros propios, ponderación adaptativa, grados de libertad equivalentes) y capítulo 8 (cálculo de las secuencias de Slepian); los autovalores de la tabla 382 anclan el oráculo de las ventanas en la batería de tests. ISBN 978-0-521-43541-3.
Thomson, D. J. (1982). Spectrum estimation and harmonic analysis. Proceedings of the IEEE, 70(9), 1055-1096. https://doi.org/10.1109/PROC.1982.12433El método multitaper: ventanas de Slepian, autoespectros propios y los pesos adaptativos que implementa multitaper_psd.
Welch, P. D. (1967). The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2), 70-73. https://doi.org/10.1109/TAU.1967.1161901La fórmula de varianza con segmentos solapados tras el número efectivo de promedios (Bendat y Piersol, sección 11.5.2.2, ref. 11).

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