Source: https://jmrplens.github.io/phonometry/es/signals/spectra/spectral-analysis/
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](/phonometry/es/signals/metrology/data-qualification/),
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.
Promediar $n_\mathrm{d}$ segmentos independientes da al estimador $2 n_\mathrm{d}$ grados de
libertad chi-cuadrado (ec. 8.162), de donde sale todo lo demás:
$$
\varepsilon_\mathrm{r}[\hat{G}_{xx}] = \frac{1}{\sqrt{n_\mathrm{d}}}, \qquad
\frac{n\,\hat{G}_{xx}}{\chi^2_{n;\,\alpha/2}} \le G_{xx} \le
\frac{n\,\hat{G}_{xx}}{\chi^2_{n;\,1-\alpha/2}}, \quad n = 2 n_\mathrm{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.
```python
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_\mathrm{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_\mathrm{r}$, el sesgo normalizado de primer orden es la forma
cerrada de la ec. 8.141, expuesta como `resolution_bias_error`:
$$
\varepsilon_\mathrm{b}[\hat{G}_{xx}(f_\mathrm{r})] \approx -\frac{1}{3}\left(\frac{B_\mathrm{e}}{B_\mathrm{r}}\right)^2 .
$$
```python
from phonometry import resolution_bias_error
eps_b = resolution_bias_error(res.resolution_bandwidth, 25.0) # pico de Br = 25 Hz
```
Un $B_\mathrm{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 $\varepsilon_\mathrm{b} \ge -0{,}2$, y por
tanto $B_\mathrm{e} \le 0{,}775\,B_\mathrm{r}$; mantenerlo por debajo de una décima de decibelio
significa $B_\mathrm{e} \le 0{,}263\,B_\mathrm{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_\mathrm{r} = f/Q = 2$ Hz. Mantener su sesgo por debajo
del 2 % (0,09 dB) pide $B_\mathrm{e} \le B_\mathrm{r}\sqrt{3 \times 0{,}02} = 0{,}49$ Hz y, para la
ventana de Hann, $B_\mathrm{e} = 1{,}5\,f_\mathrm{s}/\texttt{nperseg}$, así que a 48 kHz eso es
`nperseg` $\ge 147\,000$, es decir, segmentos de al menos **3,1 s**. Pedir
después un error aleatorio del 10 % necesita $n_\mathrm{d} = 100$ promedios efectivos, lo
que con un solape del 50 % significa un registro de unos **150 s**:
```python
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
```
*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_\mathrm{e}^2$ y el error aleatorio sube como $n_\mathrm{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](/phonometry/es/signals/metrology/data-qualification/).
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:
```python
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\log_{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\log_{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](/phonometry/es/signals/metrology/calibration/).
*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
```python
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_{xy}(f)$ compleja entre dos canales con
el mismo núcleo de Welch, e informa de la coherencia ordinaria
$\gamma^2_{xy} = |G_{xy}|^2/(G_{xx}\,G_{yy})$ 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):
$$
\varepsilon_\mathrm{r}[|\hat{G}_{xy}|] = \frac{1}{|\gamma_{xy}|\sqrt{n_\mathrm{d}}}, \qquad
\mathrm{s.d.}[\hat{\theta}_{xy}] =
\frac{\left(1-\gamma^2_{xy}\right)^{1/2}}{|\gamma_{xy}|\sqrt{2 n_\mathrm{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 $\tau_\mathrm{g} = -\mathrm{d}\varphi/(2\pi\,\mathrm{d}f)$; para un
camino de retardo puro la fase es lineal y esa pendiente lee directamente el
retardo de propagación.
```python
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
```
*La densidad espectral cruzada de un camino con retardo de 2 ms: la fase
desenrollada es exactamente la recta $-2\pi f\tau$, 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
```python
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_{vv} = \gamma^2_{xy}\,G_{yy}, \qquad
G_{nn} = \left(1-\gamma^2_{xy}\right) G_{yy}, \qquad
\mathrm{SNR}(f) = \frac{\gamma^2_{xy}}{1-\gamma^2_{xy}} .
$$
`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:
$$
\varepsilon_\mathrm{r}[\hat{G}_{vv}] =
\frac{\left(2-\gamma^2_{xy}\right)^{1/2}}{|\gamma_{xy}|\sqrt{n_\mathrm{d}}}, \qquad
\varepsilon_\mathrm{r}[\widehat{\mathrm{SNR}}] = \frac{\sqrt{2}}{|\gamma_{xy}|\sqrt{n_\mathrm{d}}} .
$$
Para ruido aditivo no correlacionado en la salida con nivel conocido, la
coherencia tiene la forma cerrada $\gamma^2 = \mathrm{SNR}/(1+\mathrm{SNR})$,
lo que hace toda la
cadena verificable con una señal sintética:
```python
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
```
*El reparto exacto del modelo del fragmento: $G_{vv} = \gamma^2 G_{yy}$ explica la
salida salvo el resto de ruido plano $G_{nn}$, y la SNR espectral se dispersa
alrededor de su forma cerrada $10\log_{10}(0{,}64/0{,}25) = 4{,}1$ dB en toda
frecuencia.*
Mostrar el código de esta figura
```python
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[\hat{\gamma}^2] \approx (1-\gamma^2)^2/n_\mathrm{d}$ (ec. 9.75):
despreciable en cuanto $n_\mathrm{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 $\gamma^2 = \mathrm{SNR}/(1+\mathrm{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](/phonometry/es/signals/spectra/correlation-delay/) |
| 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 $\varepsilon_\mathrm{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](/phonometry/es/signals/spectra/miso-coherence/): una entrada no puede explicar dos fuentes |
Dos advertencias de cierre. La estimación de coherencia está sesgada **al alza**
en $(1-\gamma^2)^2/n_\mathrm{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/2n},\ f\cdot 2^{+1/2n}]$ 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).
```python
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 $\Delta f$ se suaviza al nivel en forma cerrada
$P\,\Delta f / (f_0(2^{1/2n} - 2^{-1/2n}))$ 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_\mathrm{c}$ tiene un ancho de banda $B_n = f_\mathrm{c}\,(2^{1/2n} - 2^{-1/2n})$, es
decir, $0{,}2316\,f_\mathrm{c}$ en un tercio de octava y $0{,}7071\,f_\mathrm{c}$ en una octava
completa, así que
$$
L_\text{banda} = 10\log_{10}\!\left(\frac{G(f_\mathrm{c})\,B_n}{p_0^2}\right),
$$
que a 1 kHz en tercio de octava añade $10\log_{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](/phonometry/es/signals/levels/levels/).
## 5. Generadores de ruido de colores
`noise_signal` produce ruido gaussiano cuya PSD sigue
$G_{xx}(f) \propto f^\alpha$ 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_\text{ref})^{\alpha/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 | $\alpha$ | 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 |
```python
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\log_{10} B_n$ con $B_n \propto f_\mathrm{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](/phonometry/es/signals/levels/levels/) vive el lado de las bandas.
*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\log_{10}(\mathrm{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\log_{10}(\mathrm{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.
```python
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\log_{10}(N\sin(\pi/2N))$, el núcleo de Dirichlet evaluado medio bin fuera del
centro.
*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
```python
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}^{(mt)}(f) = \frac{1}{K}\sum_{k=0}^{K-1} \hat{S}_k(f), \qquad
\hat{S}_k(f) = \Delta t\,\Bigl|\sum_{t=1}^{N} h_{t,k}\,x_t\,
e^{-i 2\pi f t \Delta t}\Bigr|^2 .
$$
Como las ventanas son ortogonales, los autoespectros están casi
no correlacionados, así que el promedio acarrea unos $2K$ 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 = NW$ (4 por defecto), y para un registro de $N_\mathrm{r}$
muestras a $f_\mathrm{s}$ fija el semiancho $W = p\,f_\mathrm{s}/N_\mathrm{r}$ en hercios. $2W$ es la
resolución de la estimación (informada como `resolution_bandwidth`), y solo las
ventanas por debajo del número de Shannon $2p$ mantienen la energía de su ventana
espectral dentro de la banda de diseño - sus concentraciones $\lambda_k$ se
informan como `eigenvalues`, y el número de ventanas por defecto es
$K = 2p - 1$, todas las de concentración casi unidad. Un $p$ mayor
admite más ventanas (menos varianza) a costa de resolución.
```python
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)}{\lambda_k S(f) + (1-\lambda_k)\,\sigma^2 \Delta t},
\qquad
\hat{S}^{(amt)}(f) =
\frac{\sum_k b_k^2(f)\,\lambda_k\,\hat{S}_k(f)}
{\sum_k b_k^2(f)\,\lambda_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,
$\nu(f) = 2\left(\sum_k d_k\right)^2/\sum_k d_k^2$ con
$d_k = b_k^2\lambda_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 $2W$).
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.
*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_\mathrm{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
```python
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](/phonometry/es/devices/electroacoustics/electroacoustics/) forman los
mismos espectros cruzados en una función de transferencia ($H_1 = G_{xy}/G_{xx}$,
insesgado cuando el ruido está en la salida, y $H_2 = G_{yy}/G_{yx}$, insesgado
cuando está en la entrada), y la sonda de
[intensidad acústica](/phonometry/es/devices/emission/intensity/) de dos micrófonos
forma la parte imaginaria de $G_{xy}$. 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](/phonometry/es/signals/spectra/miso-coherence/),
que extiende la coherencia ordinaria a varias entradas correlacionadas y una
salida.
## Qué cubre esta guía
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.
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](/phonometry/es/devices/electroacoustics/electroacoustics/) y
[intensidad acústica](/phonometry/es/devices/emission/intensity/), no aquí. Lo
mismo ocurre con la coherencia múltiple y parcial, en la página de [coherencia
MISO](/phonometry/es/signals/spectra/miso-coherence/). 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](/phonometry/es/signals/metrology/data-qualification/): la estacionariedad que da por supuesta todo intervalo chi-cuadrado de esta página.
- [Calibración y dBFS](/phonometry/es/signals/metrology/calibration/): de dónde sale el factor que convierte las densidades de esta página en Pa²/Hz.
- [Coherencia MISO](/phonometry/es/signals/spectra/miso-coherence/): la misma matriz de espectros cruzados con varias entradas correlacionadas.
- [Análisis tiempo-frecuencia](/phonometry/es/signals/spectra/time-frequency/): la misma segmentación de Welch representada en lugar de promediada.
- [Niveles](/phonometry/es/signals/levels/levels/): el lado de los niveles de banda de la conversión densidad/banda del §4.
- Referencia de la API: [`signals.spectra`](/phonometry/es/reference/api/signals/spectra/), [`signals.windows`](/phonometry/es/reference/api/signals/windows/) y [`signals.multitaper`](/phonometry/es/reference/api/signals/multitaper/).
- Teoría: [Resolución frecuencial vs separación de bins FFT](/phonometry/es/reference/theory/signal-analysis/#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.