Coherencia múltiple y parcial

Cuando varias fuentes parcialmente correladas excitan una misma respuesta, la coherencia ordinaria de cada fuente con la salida engaña: una fuente que solo correla con la causa real hereda una coherencia espuria a través de ella, de modo que leer solo las coherencias ordinarias puede atribuir el mérito a la fuente equivocada. Bendat y Piersol, Random Data (4.a ed., 2010, capítulo 7), lo resuelven para un sistema de entradas múltiples y salida única (MISO) con las funciones de coherencia múltiple y parcial. miso_coherence las calcula desde el mismo núcleo de espectros cruzados de Welch que el resto de phonometry.metrology, para varias entradas correladas y una salida.

Ver el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from phonometry import miso_coherence, noise_signal

fs = 8192.0
# La entrada 1 excita un camino de baja frecuencia; la entrada 2 =
# 0.7*x1 + ruido independiente excita un camino de alta frecuencia, de modo
# que la entrada 2 esta correlada con la entrada 1.
x1 = noise_signal(fs, 32.0, color="white", seed=1)
x2 = 0.7 * x1 + noise_signal(fs, 32.0, color="white", seed=2)
low = signal.butter(4, 400.0, fs=fs, output="sos")
high = signal.butter(4, 1500.0, btype="high", fs=fs, output="sos")
noise = noise_signal(fs, 32.0, color="white", rms=0.05, seed=3)
y = signal.sosfilt(low, x1) + signal.sosfilt(high, x2) + noise

res = miso_coherence([x1, x2], y, fs, nperseg=2048)
f = res.frequencies
band = (f >= 20.0) & (f <= 4000.0)

fig, (ax_top, ax_bot) = plt.subplots(2, 1, figsize=(10, 7.4), sharex=True)
db = lambda v: 10 * np.log10(v)
ax_top.semilogx(f[band], db(res.output_psd[band]), color="gray",
                label="Salida medida")
for i, color in ((0, "#1f77b4"), (1, "#2ca02c")):
    ax_top.semilogx(f[band], db(res.coherent_output_spectra[i][band]),
                    color=color, label=f"Contribución de la entrada {i + 1}")
ax_top.semilogx(f[band], db(res.noise_psd[band]), "--", color="#d62728",
                label="Ruido residual")
ax_top.set_ylabel("Salida coherente [dB re 1/Hz]")
ax_top.legend()

ax_bot.semilogx(f[band], res.ordinary_coherence[1][band], ":", color="#2ca02c",
                label="Entrada 2 ordinaria (inflada por x1)")
ax_bot.semilogx(f[band], res.partial_coherence[1][band], color="#2ca02c",
                label="Entrada 2 parcial (x1 eliminada)")
ax_bot.semilogx(f[band], res.multiple_coherence[band], color="black",
                label="Múltiple")
ax_bot.set_xlabel("Frecuencia [Hz]")
ax_bot.set_ylabel("Coherencia")
ax_bot.set_ylim(0, 1.05)
ax_bot.legend()
plt.show()

El método .plot() dibuja los mismos dos paneles en una sola llamada:

res.plot()          # etiquetas en inglés; res.plot(language="es") en español

1. Coherencia ordinaria, múltiple y parcial

Para un sistema con q entradas x1..xq y salida y, miso_coherence estima todos los autoespectros y espectros cruzados por el método de Welch y reporta tres funciones de coherencia.

La coherencia ordinaria de la entrada i con la salida, por sí sola (Ec. 7.109), es la magnitud clásica de una sola entrada:

$${\gamma }_{iy}^2(f)=\frac{|G_{iy}(f){|}^2}{G_{ii}(f)G_{yy}(f)}.$$

La coherencia múltiple (Ec. 7.35) es la fracción del autoespectro de salida explicada linealmente por todas las entradas conjuntamente. A partir de la matriz hermítica de espectros cruzados de entrada Gxx y el vector entrada-salida Gxy es

$${\gamma }_{y:x}^2(f)=\frac{G_{xy}^H{G}_{xx}^{-1}{G}_{xy}}{G_{yy}}=1-\frac{G_{nn}}{G_{yy}},$$

donde Gnn es el espectro de salida residual no correlado con ninguna entrada. Si la salida lleva ruido aditivo no correlado de relación señal-ruido SNR por banda, esto vale exactamente $\mathrm{SNR}/(1+\mathrm{SNR})$, el oráculo de forma cerrada con que se verifica la implementación.

La coherencia parcial de la entrada i (Ec. 7.87) es su coherencia con la salida una vez eliminado el efecto lineal de las entradas anteriores a ella en el orden de condicionamiento:

$${\gamma }_{iy\cdot (i-1)!}^2(f)=\frac{|G_{iy\cdot (i-1)!}(f){|}^2}{G_{ii\cdot (i-1)!}(f)G_{yy}(f)}.$$

La definición de la 4.a edición mantiene la salida total Gyy en el denominador (no la salida condicionada). Esa elección hace que las coherencias parciales de las entradas ordenadas sumen la coherencia múltiple, ${\gamma }_{y:x}^2={\sum }_i{\gamma }_{iy\cdot (i-1)!}^2$ (Ec. 7.116), y se reduzcan exactamente a las coherencias ordinarias cuando las entradas son mutuamente no correladas (Ec. 7.117).

res = miso_coherence([x1, x2], y, fs)
res.ordinary_coherence   # forma (q, F): cada entrada por sí sola
res.multiple_coherence   # forma (F,):  todas las entradas conjuntamente
res.partial_coherence    # forma (q, F): condicionada a las entradas previas

2. Condicionamiento: separar causa de correlación

Las coherencias parciales provienen de los espectros condicionados Gij·r!, calculados por la recursión de eliminación de Gauss de la sección 7.3. Eliminar el efecto lineal de una entrada pivote r de cada registro restante es un paso de complemento de Schur (Ec. 7.94):

$$G_{ij\cdot r!}=G_{ij\cdot (r-1)!}-\frac{G_{rj\cdot (r-1)!}{G}_{ir\cdot (r-1)!}}{G_{rr\cdot (r-1)!}},$$

aplicado a cada entrada en el orden de condicionamiento hasta que solo queda el espectro de salida residual Gyy·q!. Este es el sentido del método: en la figura anterior la entrada 2 es 0.7·x1 + ruido independiente, así que correla con la entrada 1 pero no excita ningún camino de baja frecuencia. Su coherencia ordinaria con la salida marca por ello unos 0,3 en la banda baja, prestada por completo a través de la entrada 1, mientras que su coherencia parcial se desploma a cero ahí una vez condicionada fuera la entrada 1.

low = (res.frequencies > 100.0) & (res.frequencies < 300.0)
res.ordinary_coherence[1][low].mean()   # ~0.32  (inflada por x1)
res.partial_coherence[1][low].mean()    # ~0.00  (x1 eliminada)

3. Qué fuente domina cada banda

El condicionamiento también reparte la potencia de salida fuente por fuente. El espectro de salida coherente parcial de la entrada i (Ec. 7.86) es la porción del autoespectro de salida que aporta, y las porciones más el ruido residual reconstruyen la salida exactamente (Ecs. 7.88/7.121):

$$G_{v_i}={\gamma }_{iy\cdot (i-1)!}^2{G}_{yy},\sum _i{G}_{v_i}+G_{nn}=G_{yy}.$$

Comparar las porciones banda por banda responde a “¿qué fuente domina aquí?”. dominant_input() devuelve, para cada frecuencia, el índice de la entrada con la mayor porción:

res.coherent_output_spectra        # forma (q, F): Gvi por entrada
dominant = res.dominant_input()    # índice de la fuente más fuerte por bin
f = res.frequencies
dominant[np.argmin(abs(f - 200.0))]    # 0 (la entrada 1 excita la banda baja)
dominant[np.argmin(abs(f - 2500.0))]   # 1 (la entrada 2 excita la banda alta)

4. Calidad estadística y el orden de condicionamiento

Los errores aleatorios siguen la sección 9.3. Condicionar a las i-1 entradas previas cuesta i-1 grados de libertad, de modo que la i-ésima entrada ordenada lleva nd-(i-1) promedios efectivos (Ecs. 9.100/9.101) y la coherencia múltiple de q entradas lleva nd-(q-1) (Ec. 9.98). El resultado expone multiple_coherence_random_error y coherent_output_random_error junto al número de promedios efectivos n_averages.

Las coherencias ordinaria y múltiple no dependen del orden de condicionamiento, pero las parciales y la descomposición de la salida coherente sí: cada entrada se condiciona a lo que la precede. A falta de un orden físico, Bendat y Piersol (sección 7.2.4) recomiendan ordenar las entradas por coherencia ordinaria descendente con la salida. Pasa order para elegirlo:

res = miso_coherence([x1, x2, x3], y, fs, order=(2, 0, 1))
res.order            # (2, 0, 1): el orden aplicado
res.plot(language="es")   # los dos paneles, recalculados en el orden aplicado

Merece la pena nombrar dos escollos prácticos. Promediado: cada paso de condicionamiento gasta un grado de libertad, así que con pocos segmentos las coherencias parciales de las últimas entradas del orden son los números menos fiables de la página; promedia con generosidad antes de leer una coherencia parcial pequeña como cero. Sesgo por retardo: como la coherencia ordinaria, la estimación de Welch queda sesgada a la baja cuando un retardo global entre una entrada y la salida se convierte en una fracción apreciable de la longitud del segmento; elimina primero los retardos de propagación conocidos (ver Correlación y retardo) y mantén nperseg bien por encima del mayor retardo restante.

Los estimadores comparten el núcleo de Welch de la página de análisis espectral calibrado (mismo enventanado de Hann, solape del 50 % y calibración sin quitar tendencia), de modo que una coherencia MISO y un power_spectral_density calculados con la misma longitud de segmento son consistentes bin a bin.

Qué cubre esta guía

Cubierto. Las funciones de coherencia ordinaria, múltiple y parcial del capítulo 7 de Bendat y Piersol (ecuaciones 7.109, 7.35, 7.87), la recursión de condicionamiento por eliminación de Gauss (ecuación 7.94), los espectros de salida coherente parciales y su descomposición exacta de la salida (ecuaciones 7.86/7.116/7.121), y los errores aleatorios de la sección 9.3 para las estimaciones múltiples y condicionadas (ecuaciones 9.98-9.101), todo expuesto por miso_coherence y MISOCoherenceResult.

No cubierto. Las estimaciones de respuesta en frecuencia de entradas múltiples del capítulo (los factores de ganancia Liy), el modelo lineal que predice la salida a partir de las entradas condicionadas, se calculan internamente para construir los espectros de salida coherente pero no se devuelven: MISOCoherenceResult solo expone el lado de la coherencia y la descomposición espectral, no una función de transferencia multientrada ajustada. El método también es de salida única (MISO); un sistema con varias salidas correladas necesita una llamada a miso_coherence por cada salida.

Véase también

• Análisis espectral: el espectro de salida coherente de una entrada que esta página generaliza, y el núcleo de Welch compartido.
• Correlación y retardo: estimar y eliminar los retardos globales que sesgan la coherencia a la baja.
• Multicanal y rendimiento: la ruta por canal que este análisis cruzado complementa.
• Referencia de la API: metrology.miso.

Referencias

• Bendat, J. S. y Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4th ed.). Wiley. https://doi.org/10.1002/9781118032428Capítulo 7 (relaciones de entradas y salidas múltiples: coherencia ordinaria Ec. 7.109, coherencia múltiple Ec. 7.35, coherencia parcial Ec. 7.87, la recursión de eliminación de Gauss de los espectros condicionados Ec. 7.94, la descomposición de la salida coherente Ecs. 7.86/7.116/7.121) y sección 9.3 (errores estadísticos de las estimaciones múltiples y condicionadas, Ecs. 9.98-9.101). ISBN 978-0-470-24877-5.