Coherencia múltiple y parcial

Cuando varias fuentes parcialmente correlacionadas excitan una misma respuesta, la coherencia ordinaria de cada fuente con la salida engaña: una fuente que solo se correlaciona 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.ª 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.signals, para varias entradas correlacionadas y una salida.

Toda la razón de ser de esta guía en una sola medición. La entrada 2 no excita nada por debajo de 400 Hz y, aun así, su coherencia ordinaria con la salida promedia ahí 0,324, heredada por entero a través de su correlación con la entrada 1. Si se condiciona fuera la entrada 1, su coherencia parcial en esa misma banda cae por debajo de 1×10⁻⁴: a 200 Hz, 2,6×10⁻⁶. La coherencia múltiple se mantiene en 0,997 en toda esa banda, así que el modelo no tiene la culpa: quien leyera solo las coherencias ordinarias se gastaría el dinero en silenciar la fuente equivocada.

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 está correlacionada 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 $x_1..x_q$ y salida $y$, miso_coherence estima todos los autoespectros y espectros cruzados por el método de Welch y da tres funciones de coherencia.

La historia del condicionamiento cabe en un diagrama de bloques: dos entradas correlacionadas excitan sus propios caminos hacia una única salida medida, y la matriz de espectros cruzados de Welch es la que convierte la maraña en números por fuente.

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 $G_{xx}$ y el vector entrada-salida $G_{xy}$ es

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

donde $G_{nn}$ es el espectro de salida residual no correlacionado con ninguna entrada. Si la salida lleva ruido aditivo no correlacionado de relación señal-ruido $\mathrm{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. La notación de los subíndices es la del libro, y el signo de admiración no es un factorial: un $\cdot r!$ final significa condicionado a las entradas 1 a $r$, en el orden elegido. Así que $G_{ij\cdot r!}$ es el espectro cruzado de $i$ y $j$ con el efecto lineal de las $r$ primeras entradas eliminado, y $G_{yy\cdot q!}$, con el que termina el §2, es el espectro de salida con todas las entradas eliminadas, es decir, el ruido residual $G_{nn}$ de la coherencia múltiple de arriba. En símbolos:

$${\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.ª edición mantiene la salida total $G_{yy}$ 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 correlacionadas (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 $G_{ij\cdot 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 $G_{yy\cdot q!}$. Este es el sentido del método: en la figura anterior la entrada 2 es $0,7x_1+\text{ruido independiente}$, así que se correlaciona 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)

Cómo se instrumenta una medición MISO

Todo lo anterior corre sobre arrays sintéticos, y todo el valor del método está en máquinas reales, en una sala real, donde las fuentes no se pueden apagar de una en una, que es justo para lo que existe la técnica. Lo que separa es lo que separan las referencias, así que los sensores de referencia son la medición.

• Una referencia por fuente, colocada de modo que oiga la suya y lo menos posible de las demás. Un acelerómetro atornillado a la pata de una máquina, un micrófono a 0,2 o 0,5 m de la superficie de la carcasa y no en medio de la sala. Una referencia que recoge a su vecina es la correlación entre entradas que el condicionamiento tiene que deshacer después, y solo la deshace hasta donde la aritmética lo permite (ver más abajo).
• Un único frontal de muestreo simultáneo para todas las referencias y el receptor. El condicionamiento se construye enteramente con espectros cruzados, así que un desfase entre canales los inclina todos; un convertidor que secuencia sus canales mete ese desfase de regalo.
• Ganancias fijas durante toda la tirada, y todas las fuentes funcionando a la vez en el punto de operación que se está diagnosticando.
• Promedios de sobra. Cada paso de condicionamiento gasta uno, así que la última entrada del orden lleva $n_{\mathrm{d}}-(q-1)$; apunta a unos cuantos cientos de promedios efectivos antes de leer una coherencia parcial pequeña como un cero físico.
• Anota junto al resultado: modelo, posición y sensibilidad del sensor de cada canal, la velocidad y la carga de la máquina, la longitud de segmento y el orden de condicionamiento realmente aplicado (res.order); sin este último, el reparto no es reproducible.

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 $n_{\mathrm{d}}-(i-1)$ promedios efectivos (Ecs. 9.100/9.101) y la coherencia múltiple de $q$ entradas lleva $n_{\mathrm{d}}-(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. Supón que junto a las dos de arriba se registra una tercera fuente débilmente acoplada, un ventilador de la sala contigua, que no es coherente con ninguna de las dos. Pasa order para saltarte el orden por defecto:

x3 = noise_signal(fs, 32.0, color="white", seed=4)   # la tercera fuente, débil
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

Las coherencias ordinaria y múltiple salen idénticas diga lo que diga order; las parciales y las porciones por fuente no, porque a cada entrada solo se le acredita lo que las anteriores dejaron sin explicar.

Comprueba antes las entradas entre sí

El condicionamiento invierte la matriz de espectros cruzados de entrada, así que solo separa fuentes que sean de verdad separables. Cada paso de Schur divide por la potencia que a la entrada pivote le queda después de eliminar las anteriores, de modo que dos referencias que son casi la misma señal apenas dejan nada en ese pivote y el paso siguiente amplifica sin cota el ruido de estimación. Bendat y Piersol tratan el límite como un caso aparte (su figura 7.3): una coherencia entre entradas igual a uno significa que los dos registros son una sola fuente que llega a la salida por dos caminos, y entonces no existe ningún reparto entre dos entradas.

Esta implementación no lanza ningún error en ese caso. Un pivote condicionado cuya potencia cae por debajo de un suelo relativo se salta, y la entrada aporta exactamente cero, que es idéntico al resultado estrella de esta página, una fuente que solo se correlaciona, incluso cuando ambas fuentes son reales y lo único mal puesto eran las referencias.

Así que haz la comprobación antes de creerte el reparto:

1. Calcula la coherencia ordinaria entre cada par de entradas, un cross_spectral_density por pareja, antes de llamar a miso_coherence.
2. En cualquier banda donde una pareja pase de 0,9 aproximadamente, no leas ni las coherencias parciales ni las porciones por fuente. La coherencia múltiple y $G_{yy}$ siguen siendo válidas ahí; la atribución no.
3. Confírmalo reordenando. Una descomposición bien condicionada apenas se mueve al reordenar las entradas; una mal condicionada se mueve mucho.

Los remedios son físicos, no numéricos: acercar una referencia a su propia fuente para que deje de oír a la otra, hacer funcionar las fuentes de una en una si la máquina lo permite, o aceptar que la pareja es una sola entrada y fundirlas.

Merece la pena nombrar otros dos escollos. 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 eliminación de tendencia), de modo que una coherencia MISO y un power_spectral_density calculados con la misma longitud de segmento son congruentes 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 $L_{iy}$), 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 correlacionadas 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: signals.miso.

Referencias

• Bendat, J. S. y Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4.ª 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.