<!-- canonical: https://jmrplens.github.io/phonometry/es/guides/spectral-analysis/ -->
Source: https://jmrplens.github.io/phonometry/es/guides/spectral-analysis/

Un espectro sin su incertidumbre es media medición. Esta página cubre los
estimadores espectrales de Welch de `phonometry.metrology` 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
`Gxx(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).

Promediar `nd` segmentos independientes da al estimador `2·nd` grados de
libertad chi-cuadrado (ec. 8.162), de donde sale todo lo demás:

$$
\varepsilon_r[\hat{G}_{xx}] = \frac{1}{\sqrt{n_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_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

res = power_spectral_density(signal, 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 `Be` (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 `Br`, el sesgo normalizado de primer orden es la forma
cerrada de la ec. 8.141, expuesta como `resolution_bias_error`:

$$
\varepsilon_b[\hat{G}_{xx}(f_r)] \approx -\frac{1}{3}\left(\frac{B_e}{B_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 `Be` 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.

<details>
<summary>Mostrar el código de esta figura</summary>

```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()
```

</details>

## 2. Densidad espectral cruzada

`cross_spectral_density` estima la `Gxy(f)` compleja entre dos canales con
el mismo núcleo de Welch, e informa de la coherencia ordinaria
`γ²xy = |Gxy|²/(Gxx·Gyy)` 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_r[|\hat{G}_{xy}|] = \frac{1}{|\gamma_{xy}|\sqrt{n_d}}, \qquad
\mathrm{s.d.}[\hat{\theta}_{xy}] =
\frac{\left(1-\gamma^2_{xy}\right)^{1/2}}{|\gamma_{xy}|\sqrt{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π·df)`; 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π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.*

<details>
<summary>Mostrar el código de esta figura</summary>

```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()
```

</details>

## 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_r[\hat{G}_{vv}] =
\frac{\left(2-\gamma^2_{xy}\right)^{1/2}}{|\gamma_{xy}|\sqrt{n_d}}, \qquad
\varepsilon_r[\widehat{\mathrm{SNR}}] = \frac{\sqrt{2}}{|\gamma_{xy}|\sqrt{n_d}} .
$$

Para ruido aditivo no correlacionado en la salida con nivel conocido, la
coherencia tiene la forma cerrada `γ² = SNR/(1+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: `Gvv = γ²·Gyy` explica la
salida salvo el resto de ruido plano `Gnn`, y la SNR espectral se dispersa
alrededor de su forma cerrada 10·lg(0,64/0,25) = 4,1 dB en toda frecuencia.*

<details>
<summary>Mostrar el código de esta figura</summary>

```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()
```

</details>

El campo `coherence_bias` informa del pequeño sesgo positivo del estimador
de coherencia, `b[γ̂²] ≈ (1-γ²)²/nd` (ec. 9.75): despreciable en cuanto `nd`
llega a unos cientos, y otra razón para promediar con generosidad antes de
fiarse de una coherencia baja.

## 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·2^(-1/2n), f·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.

```python
from phonometry import fractional_octave_smoothing

smooth_psd = fractional_octave_smoothing(res.frequencies, res.psd, 3.0)
smooth_mag = fractional_octave_smoothing(freqs, np.abs(response), 6.0,
                                         domain="amplitude")  # una FRF |H|
smooth_db = fractional_octave_smoothing(freqs, levels, 3.0, domain="db")  # curva en 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₀·(2^{1/2n} - 2^{-1/2n}))` sobre una anchura de núcleo: el oráculo
fijado en los tests.

## 5. Generadores de ruido de colores

`noise_signal` produce ruido gaussiano cuya PSD sigue `Gxx(f) ∝ 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 |

```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.

## 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·lg(ENBW)` dB por encima de la densidad verdadera.
- **Ganancia coherente**: la ganancia en continua `Σ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·lg(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·lg(N·sin(π/2N))`, el núcleo de Dirichlet evaluado medio bin fuera del
centro.

<details>
<summary>Mostrar el código de esta figura</summary>

```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()
```

</details>

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. Recurra 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
incorrelados, 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 semiancho `W = NW·fs/N` se fija mediante el producto
duración x semiancho de banda `NW` (4 por defecto). `2W` es la resolución de
la estimación (informada como `resolution_bandwidth`), y solo las ventanas
por debajo del número de Shannon `2·NW` 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·NW - 1`, todas las de concentración casi unidad. Un `NW` mayor
admite más ventanas (menos varianza) a costa de resolución.

```python
from phonometry import multitaper_psd

res = multitaper_psd(signal, fs)                 # 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, `ν(f) = 2·(Σk dk)²/Σk dk²` con
`dk = b²k·λ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` 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.

<details>
<summary>Mostrar el código de esta figura</summary>

```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()
```

</details>

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.

## Relación con los estimadores H1/H2

Los [estimadores de respuesta en frecuencia](/phonometry/es/guides/electroacoustics/)
`transfer_function` y `coherence`, la sonda de
[intensidad acústica](/phonometry/es/guides/intensity/) de dos micrófonos 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 H1 calculadas con la misma longitud de segmento
son mutuamente consistentes bin a bin. La misma matriz de espectros cruzados
sustenta la [coherencia múltiple y parcial](/phonometry/es/guides/miso-coherence/),
que extiende la coherencia ordinaria a varias entradas correladas 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](/phonometry/es/guides/electroacoustics/) y
[intensidad acústica](/phonometry/es/guides/intensity/), no aquí. Lo mismo
ocurre con la coherencia múltiple y parcial, en la página de
[coherencia MISO](/phonometry/es/guides/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.
