Análisis espectral calibrado

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:

$${\epsilon }_r[{\hat{G}}_{xx}]=\frac{1}{\sqrt{n_d}},\frac{n{\hat{G}}_{xx}}{{\chi }_{n;\alpha /2}^2}\le G_{xx}\le \frac{n{\hat{G}}_{xx}}{{\chi }_{n;1-\alpha /2}^2},n=2n_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

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:

$${\epsilon }_b[{\hat{G}}_{xx}(f_r)]\approx -\frac{1}{3}{\left(\frac{B_e}{B_r}\right)}^2.$$

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.

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 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):

$${\epsilon }_r[|{\hat{G}}_{xy}|]=\frac{1}{|{\gamma }_{xy}|\sqrt{n_d}},s.d.[{\hat{\theta }}_{xy}]=\frac{{\left(1-{\gamma }_{xy}^2\right)}^{1/2}}{|{\gamma }_{xy}|\sqrt{2n_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.

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.

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_{vv}={\gamma }_{xy}^2{G}_{yy},G_{nn}=\left(1-{\gamma }_{xy}^2\right)G_{yy},\mathrm{SNR}(f)=\frac{{\gamma }_{xy}^2}{1-{\gamma }_{xy}^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:

$${\epsilon }_r[{\hat{G}}_{vv}]=\frac{{\left(2-{\gamma }_{xy}^2\right)}^{1/2}}{|{\gamma }_{xy}|\sqrt{n_d}},{\epsilon }_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:

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

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.

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[γ̂²] ≈ (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.

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

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.

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.

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. 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),{\hat{S}}_k(f)=\Delta t|\sum _{t=1}^N{h}_{t,k}{x}_t{e}^{-i2\pi ft\Delta t}{|}^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.

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},{\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.

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.

Relación con los estimadores H1/H2

Los estimadores de respuesta en frecuencia transfer_function y coherence, la sonda de intensidad acústica 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, 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 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.

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