Correlación, retardo y envolvente

Donde los estimadores espectrales calibrados describen una señal en frecuencia, esta página cubre sus equivalentes en el dominio del tiempo dentro de phonometry.signals: estimaciones de autocorrelación y correlación cruzada con las tres normalizaciones estándar y sus errores aleatorios de Bendat y Piersol; estimación del retardo (TDE) por el correlador directo, la pendiente de fase del espectro cruzado y la correlación cruzada generalizada (GCC) de Knapp y Carter con las ponderaciones Roth, SCOT, PHAT y de máxima verosimilitud; localización submuestral del pico para retardos y alineación de respuestas al impulso; y la envolvente de Hilbert con fase y frecuencia instantáneas. Los estimadores GCC corren sobre el mismo núcleo de Welch que las densidades espectrales, así que ambas vistas de un par de señales son mutuamente congruentes bin a bin.

1. Estimaciones de correlación

correlation calcula la auto- o la correlación cruzada vía FFT con relleno de ceros para que el producto circular nunca se enrolle (Bendat y Piersol, sección 11.4.2), con el convenio de signos del modelo de retardo del libro: para $y(t)=\alpha x(t-{\tau }_0)+n(t)$ la estimación alcanza su pico en $\tau =+{\tau }_0$ (ecuación 5.21). Hay tres normalizaciones disponibles:

• 'biased': las sumas por retardo divididas por $N$; se atenúa hacia los extremos del registro y queda acotada por $[R_{xx}(0)R_{yy}(0){]}^{1/2}$;
• 'unbiased': divididas por $N-|r|$ (ecuación 11.96), una estimación insesgada de $R_{xy}(\tau )$ cuya varianza crece hacia los extremos;
• 'coefficient': la función coeficiente de correlación ${\rho }_{xy}(\tau )=C_{xy}(\tau )/({\sigma }_x{\sigma }_y)$ en $[-1,1]$ sobre los registros sin media (ecuación 5.16).

import numpy as np
from phonometry import correlation

res = correlation(x, y, fs, normalization="coefficient", max_lag=0.05)
peak = np.argmax(res.values)
print(res.lags[peak], res.values[peak])   # retardo y su coeficiente
res.plot()

El modelo de dos sensores $y(t)=0,8x(t-{\tau }_0)+n(t)$ bajo las tres normalizaciones: la función coeficiente está acotada y tiene su pico en $+{\tau }_0$ (arriba); sobre el rango completo de retardos la estimación sesgada se atenúa hacia los extremos mientras que la insesgada paga su insesgadez con una varianza que crece allí (abajo).

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from phonometry import correlation, noise_signal

fs = 8192.0
delay = 102                                  # 12,45 ms
x = noise_signal(fs, 2.0, seed=4)
interference = noise_signal(fs, 2.0, rms=0.5, seed=5)
y = 0.8 * np.concatenate([np.zeros(delay), x[:-delay]]) + interference

res = correlation(x, y, fs, normalization="coefficient", max_lag=0.05)

# Una línea: la estimación de correlación frente al retardo en segundos:
res.plot(language="es")
plt.show()

# Las tres normalizaciones a mano, con los mismos registros:
biased = correlation(x, y, fs, normalization="biased")
unbiased = correlation(x, y, fs, normalization="unbiased")

fig, (ax_c, ax_n) = plt.subplots(2, 1)
ax_c.plot(1e3 * res.lags, res.values, label="coeficiente")
ax_c.axvline(1e3 * delay / fs, color="r", linestyle="--",
             label="retardo verdadero")
ax_c.set(xlabel="Retardo [ms]", ylabel="Correlación")
ax_n.plot(unbiased.lags, unbiased.values, "r", lw=0.5, label="insesgada")
ax_n.plot(biased.lags, biased.values, lw=0.5, label="sesgada")
ax_n.set(xlabel="Retardo [s]", ylabel="Correlación")
for ax in (ax_c, ax_n):
    ax.legend()
plt.show()

El resultado siempre lleva la función coeficiente junto a la normalización solicitada, porque el coeficiente es lo que necesitan las fórmulas de error. Para datos gaussianos de banda limitada, de ancho $B$, observados durante $T$ segundos (ecuaciones 8.109/8.112, válidas para $T\ge 10|\tau |$ y $BT\ge 5$):

$$\epsilon \left[{\hat{R}}_{xy}(\tau )\right]=\frac{{\left[1+{\rho }_{xy}^{-2}(\tau )\right]}^{1/2}}{\sqrt{2BT}},\epsilon \left[{\hat{R}}_{xx}(0)\right]=\frac{1}{\sqrt{BT}}.$$

res.random_error(signal_bandwidth) la evalúa por retardo con el coeficiente medido, y la función independiente correlation_random_error acepta un coeficiente explícito. En el ejemplo 8.5 del libro, una señal común de potencia $S$ llega a dos sensores que llevan ruidos independientes de potencia $M$ y $N_2$, así que el coeficiente de correlación en el retardo es $\rho =S/\sqrt{(S+M)(S+N_2)}$; con ambos ruidos diez veces la señal eso vale $1/11$, y para $B=100$ Hz durante $T=5$ s la fórmula devuelve $\epsilon \approx 0,35$, uno de los anclajes fijados de conformidad: correlation_random_error(1/11, 100.0, 5.0) imprime 0,349, que es el ejemplo 8.5 entero en una sola llamada.

(La $N$ de más arriba es la longitud del registro en muestras, tal como la usan las viñetas de normalización; la potencia de ruido del segundo sensor se escribe aquí $N_2$ para no confundirlas.) Dos formas cerradas anclan el propio estimador en los tests: la autocorrelación de un seno, $(A^2/2)\cos (2\pi f_0\tau )$, y la autocorrelación $\sin (2\pi B\tau )/(2\pi B\tau )$ del ruido blanco de banda limitada (ecuación 8.120).

2. Estimación del retardo

time_delay estima el retardo de $y$ respecto a $x$ en el modelo de dos sensores $y(t)=\alpha x(t-{\tau }_0)+n(t)$ (B&P, sección 5.1.4) por tres vías:

• 'direct': el pico de la función coeficiente de correlación sobre el registro completo;
• 'phase': la pendiente por mínimos cuadrados, ponderada por $|G_{xy}|$, de la fase del espectro cruzado (ecuación 5.101b): un retardo puro tiene fase exactamente lineal, así que este estimador resuelve retardos fraccionarios a mejor de 1e-3 muestras sin ninguna interpolación de pico, siempre que la fase desenrollada no sea ambigua (retardos limpios y moderados);
• 'gcc': la correlación cruzada generalizada de Knapp y Carter (1976): el espectro cruzado promediado por Welch se pondera con $\psi (f)$ antes de la transformada inversa, afilando el pico que la autocorrelación de la propia señal ensancharía (su ecuación 9).

La imagen física son dos micrófonos y un frente de onda: el camino extra hasta el micrófono lejano es $c$ por el retardo, y todo el problema de estimación es localizar un pico en el eje de retardos.

Del retardo a la geometría

Un retardo casi nunca es la respuesta: es el paso intermedio. Multiplicado por la velocidad del sonido se convierte en una diferencia de camino, y para una onda plana que llega con un ángulo $\theta$ respecto de la normal a una pareja separada $d$,

$$\sin \theta =\frac{c{\tau }_0}{d},$$

de modo que los 2,44 ms del diagrama a 343 m/s son 0,84 m de camino extra, un número que solo tiene sentido para una pareja separada al menos esa distancia. Con ello vienen tres límites, y los tres son comprobables:

• $|c{\tau }_0|>d$ no es físico. Si la estimación supera la separación, el estimador se ha enganchado al pico equivocado o a una reflexión. Poner max_delay=d/c mantiene la búsqueda dentro del margen que permite la geometría, que es la salvaguarda más barata de la página.
• La forma de onda plana pide que la fuente esté en el campo lejano de la pareja. Más cerca, la curvatura del frente de onda sesga $\theta$, y los dos micrófonos ven ángulos distintos del mismo frente.
• La resolución angular es la peor en el eje de la pareja. Derivando, $\mathrm{d}\theta =c\mathrm{d}\tau /(d\cos \theta )$, así que una incertidumbre temporal fija se convierte en un error angular cada vez mayor conforme $\theta$ se acerca a ±90°, donde $\cos \theta \to 0$. Una pareja es precisa por el través y ciega por el eje, que es la razón de que la localización de dirección use más de dos sensores.

La separación $d$ fija por tanto a la vez el margen no ambiguo de la pareja y la resolución angular que puede ofrecer, y es el primer número que hay que elegir cuando el retardo va a convertirse en una dirección.

Las ponderaciones de la Tabla I de Knapp y Carter, con las condiciones que el artículo asocia a cada una:

weighting	$\psi (f)$	Comportamiento y condiciones
'none'	$1$	El correlador simple: la delta en el retardo queda convolucionada con la autocorrelación de la señal; pico ancho con señales coloreadas.
'roth'	$1/G_{xx}$	Suprime las bandas donde el primer sensor es ruidoso; aun así ensancha salvo que ese ruido sea espectralmente similar a la señal.
'scot'	$1/\sqrt{G_{xx}{G}_{yy}}$	Preblanquea ambos canales simétricamente; coincide con Roth cuando los sensores son iguales.
'phat'	$1/|G_{xy}|$	Idealmente una delta en el retardo para ruidos incorrelacionados (su ecuación 23), pero la ponderación ignora la relación señal-ruido, así que las bandas sin señal aportan fase aleatoria de magnitud unidad. Necesita potencia de señal en toda la banda de análisis.
'ml'	${\gamma }^2/(|G_{xy}|(1-{\gamma }^2))$	El procesador de máxima verosimilitud de Hannan-Thomson: un PHAT atenuado por la varianza de fase que cada banda realmente soporta. Alcanza la cota de Cramér-Rao; la opción segura cuando la señal no llena la banda.

Las condiciones de Knapp y Carter están todas escritas para ruido de sensor: una pareja en anecoico con ruido independiente añadido en cada micrófono. La degradación con la que se topa un acústico en interiores es de otra naturaleza, porque una reflexión es una copia coherente de la señal y no ruido, y eso cambia la elección:

• Domina el ruido de sensor (una sala silenciosa, una pareja al aire libre, una base larga): 'ml' es la óptima, porque atenúa cada banda según la varianza de fase que esa banda realmente soporta.
• Domina la reverberación (mala relación directo-reverberante): 'phat' es la opción habitual precisamente porque descarta la magnitud, así que una reflexión fuerte no puede dominar el pico por ser intensa.
• No converge ninguna de las dos: enventana el registro alrededor de la llegada directa y correla solo los primeros milisegundos, que es la única maniobra que elimina una reflexión en lugar de reponderarla.

Conviene memorizar la firma del fallo, porque ningún número devuelto la delata: un correlograma cuyo pico principal no sobresale mucho de sus vecinos, o un retardo que salta entre segmentos del mismo registro, significa que el camino directo no es dominante y que el retardo que tienes es el de una reflexión. delay_std no te va a avisar: la ecuación 8.129 modela la dispersión de un único pico, no la elección entre dos.

Los dos estimadores caen en el retardo de 2,441 ms (2,4413 ms el directo, 2,4418 ms el PHAT), así que la discusión no es sobre la exactitud con un par limpio, sino sobre cuánto se puede fiar uno del pico. El preblanqueo estrecha el lóbulo principal de 0,73 ms a media altura a 0,24 ms, dos muestras a 8192 Hz, porque el correlador directo convoluciona la delta con la autocorrelación de una señal limitada en banda a 800 Hz mientras que el PHAT tira esa magnitud. Lee la anchura, no la posición: es la que dice si podría haber una reflexión escondida bajo el pico.

Mostrar el código de esta figura

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

fs = 8192.0
delay = 20  # muestras
b, a = sp_signal.butter(2, 800.0 / (fs / 2.0))   # señal común coloreada
s = sp_signal.lfilter(b, a, noise_signal(fs, 4.0, color="white", seed=10))
x = s + noise_signal(fs, 4.0, color="white", rms=0.02, seed=11)
y = np.roll(s, delay) + noise_signal(fs, 4.0, color="white", rms=0.02, seed=12)

direct = time_delay(x, y, fs, method="direct", max_delay=0.01)
phat = time_delay(x, y, fs, method="gcc", weighting="phat",
                  nperseg=2048, max_delay=0.01)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(1e3 * direct.lags,
        direct.correlation / np.max(np.abs(direct.correlation)),
        label="Correlación cruzada directa")
ax.plot(1e3 * phat.lags, phat.correlation, label="GCC-PHAT")
ax.axvline(1e3 * delay / fs, ls="--", color="k", label="Retardo verdadero")
ax.set_xlabel("Retardo [ms]")
ax.set_ylabel("Correlación normalizada")
ax.legend()
plt.show()

Los métodos de pico de correlación refinan el pico muestral con interpolación parabólica de tres puntos, opcionalmente tras un sobremuestreo local de banda limitada (upsample=16 remuestrea una ventana alrededor del pico dieciséis veces antes de la parábola). La precisión submuestral presupone que el pico está sobremuestreado, es decir, que las señales están limitadas en banda por debajo de Nyquist; en un par limitado a $0,4f_{\mathrm{s}}$ los tests fijan el error alcanzable en ≲0,1 muestras con la parábola sola y ≲2e-3 muestras con upsample=16. Para la GCC el retardo debe caber en medio segmento de Welch; aumenta nperseg para retardos largos.

Con signal_bandwidth dado, el resultado también lleva la incertidumbre de localización del pico de la ecuación 8.129 de B&P y su intervalo $\pm 2\sigma$ (ecuación 8.130):

$$\sigma ({\hat{\tau }}_0)\approx {\left(\frac{3}{4}\right)}^{1/4}\frac{\sqrt{\epsilon [{\hat{R}}_{xy}({\tau }_0)]}}{\pi B}.$$

from phonometry import time_delay

res = time_delay(x, y, fs, method="gcc", weighting="ml",
                 nperseg=2048, upsample=16, signal_bandwidth=1000.0)
print(res.delay, res.delay_samples)     # segundos y muestras fraccionarias
print(res.delay_std, res.delay_interval)  # sigma de la ec. 8.129, +/-2 sigma
res.plot()                              # correlación con el retardo marcado

La fórmula modela el pico de la función de correlación continua, así que trata el intervalo como una cota conservadora de orden de magnitud: el Monte Carlo con semilla de la batería de tests observa la dispersión real por debajo de la predicción.

3. Retardo y alineación de respuestas al impulso

La correlación cruzada de una respuesta al impulso con un impulso unidad ideal es la propia RI, así que la localización submuestral del pico de su magnitud es su tiempo de llegada. impulse_response_delay aplica exactamente el mismo refinamiento que el pico de la TDE (sobremuestreo local de banda limitada, ×8 por defecto, más la parábola), y con una RI de reference mide el retardo entre el par a partir de su correlación cruzada sobre el registro completo (los transitorios de un solo disparo no son registros estacionarios, así que se usa el correlador directo en lugar de la GCC promediada por Welch):

from phonometry import align_impulse_responses, impulse_response_delay

t_arrival = impulse_response_delay(ir, fs)              # segundos desde t = 0
dt = impulse_response_delay(ir_b, fs, reference=ir_a)   # retardo del par

res = align_impulse_responses(ir_b, ir_a, fs)  # elimina el retardo estimado
res.plot()                       # superposición referencia frente a alineada

Una RI medida con un retardo fraccionario de 7,37 muestras (trazos) se alinea de vuelta sobre su referencia: el desplazamiento de banda limitada elimina las 7,36 muestras estimadas y la traza alineada (puntos) cae sobre la referencia.

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal as sp_signal
from phonometry import align_impulse_responses, fractional_delay

fs = 48000.0
t = np.arange(int(0.03 * fs)) / fs
rng = np.random.default_rng(6)
# Pulso de referencia de banda limitada: ráfaga tonal gaussiana de 2 kHz en 5 ms.
ir_a = sp_signal.gausspulse(t - 0.005, fc=2000.0, bw=0.5)
ir_b = fractional_delay(ir_a, 7.37)[: ir_a.size]
ir_b += 0.005 * rng.standard_normal(ir_a.size)

res = align_impulse_responses(ir_b, ir_a, fs)

# Una línea: referencia y RI alineada superpuestas:
res.plot(language="es")
plt.show()

# A mano, desde los campos que lleva el resultado:
t_ms = 1e3 * t
fig, ax = plt.subplots()
ax.plot(t_ms, res.reference, lw=1.6, label="RI de referencia")
ax.plot(t_ms, ir_b, "0.5", lw=1.0, linestyle="--", label="RI medida")
ax.plot(t_ms, res.aligned[: t.size], "r:", label="RI alineada")
ax.set(xlabel="Tiempo [ms]", ylabel="Amplitud",
       title=f"retardo eliminado: {res.delay_samples:.2f} muestras")
ax.legend()
plt.show()

align_impulse_responses elimina el retardo estimado con un desplazamiento fraccionario exacto de banda limitada (una rampa de fase en el dominio de la frecuencia sobre un registro con relleno de ceros, así que nada se enrolla): la herramienta para promediar conjuntos de RI o comparar mediciones tomadas a distancias ligeramente distintas. Los tests sintéticos de retardo fraccionario documentan la precisión alcanzable sobre un pulso suave de banda limitada: en torno a 1e-2 muestras con la parábola sola, 1e-3 con el upsample=8 por defecto y por debajo de 1e-5 a ×32.

4. Envolvente de Hilbert y frecuencia instantánea

envelope construye la señal analítica $z(t)=x(t)+j\tilde{x}(t)$ mediante la construcción de espectro unilateral que recomiendan Bendat y Piersol (ecuación 13.25) y devuelve las tres magnitudes del capítulo 13 sobre un mismo eje temporal:

$$A(t)={\left[x^2(t)+{\tilde{x}}^2(t)\right]}^{1/2},\theta (t)=\arctan \frac{\tilde{x}(t)}{x(t)},f(t)=\frac{1}{2\pi }\frac{d\theta }{dt}.$$

Estas tres tienen sentido físico solo si el registro es de verdad una portadora con una modulación lenta: el ancho de banda de la modulación tiene que quedar por debajo de la frecuencia de portadora, para que los dos espectros no se solapen. Esa es la condición de banda estrecha (de Bedrosian), y es exactamente por lo que se cumple la identidad AM exacta de más abajo. Nada en la llamada la impone: envelope devuelve números para cualquier vector que le entregues.

Una violación es fácil de reconocer una vez sabes qué buscar. La envolvente deja de ser suave y empieza a seguir los semiciclos individuales de la forma de onda, y la frecuencia instantánea da picos, se sale por completo de la banda de la señal y se vuelve negativa, algo que ocurre cada vez que la señal analítica pasa cerca del origen del plano complejo, y que es un artefacto geométrico y no una frecuencia física. Aplicadas a voz, a una respuesta al impulso de sala o a ruido de banda ancha, las tres magnitudes hacen esto. El arreglo es limitar antes el registro a banda estrecha: filtrar en paso banda alrededor de la componente que interesa, exactamente como hace el espectro de la envolvente con su parámetro band=, o leer la cresta de un espectrograma cuando coexisten varias componentes, porque la construcción de Hilbert no tiene forma de representar más de una a la vez.

Para una portadora modulada en amplitud $u(t)\cos (2\pi f_{0}t)$ la envolvente recupera $u(t)$ exactamente (ecuación 13.27): la batería de conformidad fija la envolvente AM recuperada y el par $\cos \to \sin$ de la Tabla 13.1 al nivel de 1e-9, y la frecuencia instantánea de un barrido sigue su rampa.

from phonometry import envelope

res = envelope(x, fs)
print(res.envelope, res.instantaneous_frequency)
res.plot()               # señal + envolvente, frecuencia instantánea

slow = envelope(x, fs, decimation_factor=32)   # con antialias, fs/32

Las magnitudes de Hilbert de un modo de 250 Hz golpeado: la envolvente traza el decaimiento exponencial (arriba), y la frecuencia instantánea se asienta sobre la portadora mientras el modo domina, temblando a medida que la señal se hunde en el suelo de ruido (abajo, diezmado ×8 con antialias).

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from phonometry import envelope

fs = 8192.0
t = np.arange(int(0.4 * fs)) / fs
rng = np.random.default_rng(7)
x = np.exp(-t / 0.1) * np.sin(2 * np.pi * 250.0 * t)   # un modo golpeado
x += 0.001 * rng.standard_normal(t.size)

res = envelope(x, fs, decimation_factor=8)

# Una línea: señal + envolvente y la frecuencia instantánea:
res.plot(language="es")
plt.show()

# A mano, desde los campos que lleva el resultado:
fig, (ax_e, ax_f) = plt.subplots(2, 1, sharex=True)
ax_e.plot(t, res.signal, lw=0.6, label="Señal")
ax_e.plot(res.times, res.envelope, "r", lw=1.8, label="Envolvente A(t)")
ax_e.plot(res.times, -res.envelope, "r", lw=1.8)
ax_e.legend()
ax_f.plot(res.times, res.instantaneous_frequency, lw=0.9,
          label="Frecuencia instantánea f(t)")
ax_f.axhline(250.0, color="g", linestyle="--", label="portadora de 250 Hz")
ax_f.set(xlabel="Tiempo [s]", ylabel="Frecuencia [Hz]", ylim=(230, 270))
ax_f.legend()
plt.show()

La envolvente de una señal de banda limitada es a su vez de baja frecuencia, así que el resultado ofrece diezmado opcional: un filtro antialias FIR de fase cero por defecto, o submuestreo simple con antialias=False, exactamente el convenio que la cadena de sonoridad/aspereza ECMA-418-2 de phonometry.psychoacoustics aplica internamente tras su paso de banda auditivo (fórmulas 65/119 de la norma), apropiado cuando la entrada ya es de banda estrecha.

Relación con los estimadores espectrales

time_delay (métodos GCC y de fase) corre sobre el mismo núcleo de Welch (ventana, política de solape, calibración sin eliminación de tendencia, valores por defecto de segmento) que cross_spectral_density y los estimadores de respuesta en frecuencia H1/H2, así que una GCC, una coherencia y un espectro cruzado calculados con la misma longitud de segmento coinciden bin a bin; el estimador 'phase' es literalmente la pendiente de la fase del CrossSpectralDensityResult, ponderada como prescribe la ecuación 5.101b.

La envolvente de Hilbert de aquí es la compañera en el tiempo del espectro de la envolvente, que lee las mismas modulaciones como líneas discretas en frecuencia; y la alineación submuestral comparte su núcleo de banda limitada con las herramientas públicas de retardo fraccionario y remuestreo.

Qué cubre esta guía

• Cubierto

  El capítulo 5, la sección 8.4 y el capítulo 13 de Bendat y Piersol: correlation con las normalizaciones sesgada, insesgada y coeficiente y sus fórmulas de error aleatorio (ecuaciones 8.109/8.112); time_delay con los métodos directo, de pendiente de fase y la correlación cruzada generalizada de Knapp y Carter (1976), incluidas las ponderaciones Roth, SCOT, PHAT y de máxima verosimilitud de su Tabla I; la incertidumbre de localización del pico de la ecuación 8.129; impulse_response_delay y align_impulse_responses; y la envelope de Hilbert con fase y frecuencia instantáneas (capítulo 13).

• No cubierto

  correlation y time_delay modelan un único retardo de camino común entre exactamente dos sensores (el modelo $y(t)=\alpha x(t-{\tau }_0)+n(t)$ de la ecuación 5.21). Un registro con varias llegadas (multitrayecto, un camino directo más reflexiones) no se separa automáticamente: time_delay y el correlador directo solo informan del mayor pico, la misma limitación que echo_detection señala en el cepstrum dentro de Cepstrum, ecos y espectro de la envolvente. Estimar retardos entre más de dos sensores, como haría un array o un beamformer, implica llamar tú mismo a estas funciones por pares: no hay ninguna resolución de TDOA multisensor incorporada.

Véase también

• Cepstrum, ecos y espectro de la envolvente: la vía sin referencia hasta un retardo, y la envolvente que construye esta página.
• Análisis espectral calibrado: el espectro cruzado con el que se construyen las ponderaciones GCC, y las reglas de adquisición de dos canales.
• Señales de prueba: el núcleo de retardo fraccionario de banda limitada que hay detrás de la alineación submuestral.
• Promediado síncrono en el tiempo: la misma alineación aplicada por revolución.
• Referencia de la API: signals.correlation y signals.envelope.

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.1.4 y 5.2.6-5.2.7 (retardo por correlación y espectro cruzado), 8.4 (errores aleatorios de las estimaciones de correlación y de la localización del pico), 11.4 (cómputo FFT con relleno de ceros) y capítulo 13 (transformadas de Hilbert, envolvente y fase instantánea). ISBN 978-0-470-24877-5.
• Knapp, C. H. y Carter, G. C. (1976). The generalized correlation method for estimation of time delay. IEEE Transactions on Acoustics, Speech, and Signal Processing, 24(4), 320-327. https://doi.org/10.1109/TASSP.1976.1162830El marco GCC, las ponderaciones de la Tabla I con sus condiciones y el procesador de máxima verosimilitud (Hannan-Thomson).