Medición de sistemas: Golay, barridos, inversión

La guía de acústica de salas recupera respuestas al impulso con los dos caballos de batalla de la ISO 18233: el barrido exponencial y la MLS. Esta página añade la capa de ingeniería de medición a su alrededor, tomada de la literatura de medición de funciones de transferencia y no de una norma: los pares complementarios de Golay (una tercera excitación cuya deconvolución está exactamente libre de ruido de correlación), los barridos conformados cuyo espectro de magnitud sigue cualquier objetivo prescrito conservando el factor de cresta de un seno barrido, y la inversión regularizada de una respuesta medida - la herramienta que convierte la respuesta medida de un altavoz o micrófono en un ecualizador seguro, y el puente entre las dos anteriores: la inversa de una respuesta medida es el objetivo espectral natural para el barrido de la siguiente medición.

1. Pares complementarios de Golay

Un par de Golay son dos secuencias binarias a, b de longitud L = 2**n, construidas con una recursión de dos líneas - anexar b a a, y -b a a (Havelock, Parte I, Cap. 6, Ec. (1)). Su propiedad definitoria es algebraica, no aproximada: las dos autocorrelaciones periódicas suman

$$\sum _j{a}_j{a}_{j+k}+b_j{b}_{j+k}=\{\begin{array}{ll}2L& k=0\\ 0& k\ne 0,\end{array}$$

una delta exacta (Ec. (2)). La autocorrelación de una MLS tiene un residuo -1 en cada retardo no nulo; los lóbulos laterales de Golay se cancelan idénticamente. La medición emite cada código por turnos - reproducir a periódicamente, grabar un período en régimen estacionario, y lo mismo con b - y suma las dos correlaciones cruzadas circulares (Ec. (4)): para un sistema lineal e invariante en el tiempo sin ruido, la respuesta al impulso vuelve con precisión de máquina, una identidad de forma cerrada que los tests y el informe de conformidad fijan en 1e-13.

import numpy as np
from scipy import signal
from phonometry import golay_impulse_response, golay_pair

fs = 48000
pair = golay_pair(14)                       # dos códigos de 16384 muestras

# Medición simulada: reproducir cada código periódicamente y grabar un
# período en régimen estacionario de un sistema paso banda tipo sala.
b, a = signal.butter(2, [80.0, 12000.0], btype="bandpass", fs=fs)
length = pair[0].size
rec_a = signal.lfilter(b, a, np.tile(pair[0], 3))[2 * length:]
rec_b = signal.lfilter(b, a, np.tile(pair[1], 3))[2 * length:]

ir = golay_impulse_response(rec_a, rec_b, pair, fs=fs)
print(ir.method, ir.size)                   # golay 16384
ir.plot(language="es")

La recuperación de Golay de la respuesta al impulso de un sistema paso banda: la suma de las correlaciones complementarias cae sobre la respuesta verdadera con precisión de máquina (error máximo 2·10⁻¹⁴ aquí), la identidad de forma cerrada que fijan los tests.

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from phonometry import golay_impulse_response, golay_pair

fs = 48000
pair = golay_pair(14)                       # dos códigos de 16384 muestras
b, a = signal.butter(2, [200.0, 2000.0], btype="bandpass", fs=fs)
length = pair[0].size
rec_a = signal.lfilter(b, a, np.tile(pair[0], 3))[2 * length:]
rec_b = signal.lfilter(b, a, np.tile(pair[1], 3))[2 * length:]

ir = golay_impulse_response(rec_a, rec_b, pair, fs=fs)

# Una línea: la respuesta al impulso recuperada:
ir.plot(language="es")
plt.show()

# A mano, contra la respuesta verdadera del sistema:
impulse = np.zeros(length)
impulse[0] = 1.0
true_ir = signal.lfilter(b, a, impulse)
err = np.max(np.abs(np.asarray(ir) - true_ir))

t_ms = 1e3 * np.arange(length) / fs
view = t_ms <= 6.0
fig, ax = plt.subplots()
ax.plot(t_ms[view], np.asarray(ir)[view], label="RI recuperada")
ax.plot(t_ms[view], true_ir[view], "r--",
        label="Respuesta verdadera del sistema")
ax.set(xlabel="Tiempo [ms]", ylabel="Amplitud",
       title=f"máx |recuperada - verdadera| = {err:.1e}")
ax.legend()
plt.show()

El resultado es el mismo ImpulseResponseResult que devuelven los frontales de barrido y MLS, así que todo lo que viene después - parámetros de sala, curvas de caída, STI - lo consume sin cambios. Las grabaciones de varios períodos se promedian síncronamente, de modo que el ruido de fondo no correlado cae 3 dB por duplicación mientras la parte determinista permanece exacta. Los compromisos son un espejo de la MLS: la recuperación es circular (el sistema debe decaer dentro de un período del código, o un ImpulseResponseWarning señala la cola con aliasing), los productos de distorsión se reparten por el período en vez de separarse como en un barrido, y los dos regímenes estacionarios secuenciales hacen del par el más expuesto de la familia a la variación temporal (Sec. 2 de Xiang) - el precio de la complementariedad exacta.

2. Barridos con un espectro objetivo arbitrario

La energía de un barrido en una frecuencia dada puede fijarse de dos maneras: con su amplitud, o con cuánto tiempo permanece en ella. La síntesis en el dominio de la frecuencia de Mueller y Massarani (Secs. 4.2-4.3) usa la segunda palanca: definir la magnitud objetivo |H(f)|, hacer crecer el retardo de grupo en proporción a la potencia del objetivo,

$${\tau }_G(f)={\tau }_G(f-df)+C|H(f){|}^2,C=\frac{{\tau }_G(f_{end})-{\tau }_G(f_{start})}{\sum |H(f){|}^2},$$

(Ecs. (11)-(12)), integrar el retardo de grupo en una fase, y aplicar la FFT inversa. El barrido sigue entonces cualquier forma espectral con una envolvente casi constante, manteniendo el factor de cresta cerca del ideal de 3,02 dB de un seno barrido - a diferencia de una señal de ruido con el mismo espectro, que queda ~6 dB por encima. shaped_sweep_signal implementa la construcción con los detalles de limitación en banda y de fase en Nyquist del artículo; target es "pink" (el énfasis clásico de medición de salas, por defecto), "white", o cualquier par (frecuencias_hz, magnitud_db):

from phonometry import shaped_sweep_signal

fs = 48000
sweep = shaped_sweep_signal(fs, 50.0, 5000.0, 2.0, target="pink")
print(round(sweep.crest_factor_db, 1))      # 4.2 (dB; el ideal es 3,02)
sweep.plot(language="es")

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 shaped_sweep_signal

fs = 48000
res = shaped_sweep_signal(fs, 50.0, 5000.0, 2.0, target="pink")
x = np.asarray(res)

nperseg = 8192   # 75 % de solape: con 50 % un barrido riza ~2 dB
freqs, psd = sp_signal.welch(x, fs=fs, nperseg=nperseg,
                             noverlap=3 * nperseg // 4)
welch_db = 10.0 * np.log10(psd)
welch_db -= welch_db[(freqs >= 50.0) & (freqs <= 5000.0)].max()
target_db = 20.0 * np.log10(np.maximum(res.magnitude, 1e-300))
target_db -= target_db[(res.frequencies >= 50.0)
                       & (res.frequencies <= 5000.0)].max()

fig, axes = plt.subplots(2, 1, figsize=(10, 7))
axes[0].plot(np.arange(x.size) / fs, x, lw=0.5)
axes[0].set_xlabel("Tiempo [s]")
axes[0].set_ylabel("Amplitud")
axes[1].semilogx(freqs[1:], welch_db[1:], lw=1.3,
                 label="Espectro de Welch")
axes[1].semilogx(res.frequencies[1:], target_db[1:], "r--",
                 label="Objetivo rosa (-3 dB por octava)")
axes[1].axvspan(50.0, 5000.0, alpha=0.08, label="Banda del barrido")
axes[1].set_xlabel("Frecuencia [Hz]")
axes[1].set_ylabel("Nivel re máximo en banda [dB]")
axes[1].set_ylim(-60.0, 8.0)
axes[1].legend()
plt.tight_layout()
plt.show()

Los metadatos de la síntesis viajan con el resultado: frequencies y magnitude son la magnitud exacta, limitada en banda, impuesta al espectro; group_delay es la trayectoria tiempo-frecuencia del barrido, y crest_factor_db la relación pico/RMS conseguida. El propósito de la conformación es ingeniería de la relación señal-ruido: adaptar el espectro emitido al ruido de fondo del recinto (o a los límites de potencia del altavoz) compra una SNR independiente de la frecuencia que ningún posprocesado puede recuperar (Mueller y Massarani, Sec. 3).

La deconvolución no necesita nada nuevo: el resultado actúa como su propio array de referencia para el método espectral de impulse_response, que divide de nuevo la coloración del barrido:

import numpy as np
from scipy import signal
from phonometry import impulse_response, shaped_sweep_signal

fs = 48000
freqs = np.array([50.0, 200.0, 1000.0, 8000.0])
emphasis = np.array([12.0, 6.0, 0.0, 0.0])   # refuerzo LF contra el retumbe
sweep = shaped_sweep_signal(fs, 50.0, 8000.0, 3.0,
                            target=(freqs, emphasis))

# Medición simulada de un sistema conocido (sustituir por emitir/grabar):
b, a = signal.butter(2, [200.0, 4000.0], btype="bandpass", fs=fs)
excitation = np.concatenate([np.asarray(sweep), np.zeros(fs)])
recorded = signal.lfilter(b, a, excitation)

ir = impulse_response(recorded, excitation, fs, length=fs)
ir.plot(language="es")

El énfasis vuelve a ponderar el suelo de ruido de la medición, no la respuesta recuperada. Una sutileza heredada de la limitación dura en banda: el núcleo limitador de banda de la deconvolución es de fase cero, así que una pequeña cola anticausal queda al final del búfer completo de deconvolución - use return_full=True cuando persiga décimas de dB en los bordes de la banda.

3. Inversión espectral regularizada

Ecualizar con una respuesta medida significa invertirla, y un recíproco puro 1/H(f) estalla allí donde el sistema radia poca energía - un valle en banda, todo lo que queda fuera de banda - convirtiendo el ecualizador en un amplificador de ruido. Mueller y Massarani confinan la inversión con un paso banda (Secs. 3.1, 4.5); la forma general de ese confinamiento es la regularización de Tikhonov dependiente de la frecuencia de Kirkeby y Nelson,

$$H_{inv}(f)=\frac{H^{\ast }(f)}{|H(f){|}^2+\epsilon (f)},$$

con ε(f) pequeña dentro de la banda a ecualizar y grande fuera. En banda, la magnitud ecualizada |H·H_inv| se desvía de la unidad exactamente en ε/(|H|² + ε) - una forma cerrada que la suite de conformidad comprueba bin a bin - y fuera de banda la ganancia del filtro nunca puede superar 1/(2·√ε), el máximo analítico de x/(x² + ε). Un retardo de modelado de medio bloque del filtro hace causal la inversa, en general anticausal, de una respuesta de fase mixta (Kirkeby y Nelson, Sec. 2.4).

import numpy as np
from scipy import signal
from phonometry import regularized_inverse_filter

fs = 48000.0
# Una respuesta medida tipo altavoz (sustituir por una RI medida).
b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)
imp = np.zeros(2048)
imp[0] = 1.0
h = signal.lfilter(b, a, imp)

inv = regularized_inverse_filter(h, fs, f_range=(200.0, 4000.0))
print(round(inv.flatness_db, 5))    # 1e-05 (dB: desviación en banda de 0 dB)
print(round(inv.max_gain_db, 1))    # -6.0 (dB: refuerzo fuera de banda acotado)
inv.plot(language="es")

Mostrar el código de esta figura

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from phonometry import regularized_inverse_filter

fs = 48000.0
b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)
imp = np.zeros(2048)
imp[0] = 1.0
h = signal.lfilter(b, a, imp)

res = regularized_inverse_filter(h, fs, f_range=(200.0, 4000.0))
f = res.frequencies[1:]
h_mag = np.abs(res.response_spectrum)[1:]
inv_mag = np.abs(res.spectrum)[1:]
peak = h_mag.max()

fig, ax = plt.subplots(figsize=(10, 6))
ax.semilogx(f, 20 * np.log10(h_mag / peak), label="Medida $|H|$")
ax.semilogx(f, 20 * np.log10(inv_mag * peak),
            label=r"Inversa $|H_{\mathrm{inv}}|$")
ax.semilogx(f, 20 * np.log10(h_mag * inv_mag), lw=1.8,
            label=r"Ecualizada $|H \cdot H_{\mathrm{inv}}|$")
ax.axvspan(200.0, 4000.0, alpha=0.08, label="Banda ecualizada")
ax.set_ylim(-50.0, 15.0)
ax.set_xlabel("Frecuencia [Hz]")
ax.set_ylabel("Magnitud [dB]")
ax.legend()
plt.show()

Ambos niveles de regularización son fracciones del pico de |H|², generalizando la regularization escalar de impulse_response: regularization_inside (por defecto 1e-6) fija con qué exactitud se aplana la banda, regularization_outside (por defecto 1.0) limita el refuerzo fuera de banda 6 dB por debajo de la unidad en banda, y un fundido geométrico de transition_octaves octavas (por defecto 1/3) conecta ambos suavemente. El apply() del resultado ecualiza cualquier grabación con el retardo de modelado ya eliminado, y la función acepta directamente un ImpulseResponseResult - su frecuencia de muestreo viaja con él:

import numpy as np
from scipy import signal
from phonometry import (impulse_response, regularized_inverse_filter,
                        sweep_signal)

fs = 48000
sweep = sweep_signal(fs, 50.0, 20000.0, 1.0)
excitation = np.concatenate([sweep, np.zeros(fs // 2)])
b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)
recorded = signal.lfilter(b, a, excitation)     # en la práctica: emitir/grabar

ir = impulse_response(recorded, excitation, fs)  # vale cualquier frontal
inv = regularized_inverse_filter(ir, f_range=(200.0, 10000.0))
flat_recording = inv.apply(recorded)             # retardo ya eliminado
inv.plot(language="es")   # magnitudes medida, inversa y ecualizada (figura de arriba)

Cerrando el círculo: preénfasis desde una respuesta medida

Las dos mitades de esta página componen el propio flujo de trabajo de Mueller y Massarani (su Fig. 18): medir el altavoz, invertir su respuesta en la banda de transmisión y entregar la magnitud inversa a la síntesis del barrido como objetivo - la siguiente medición radia entonces un espectro acústico plano, con la ecualización hecha por la excitación en vez de por un posprocesado que amplifica ruido:

import numpy as np
from scipy import signal
from phonometry import regularized_inverse_filter, shaped_sweep_signal

fs = 48000.0
b, a = signal.butter(2, [100.0, 8000.0], btype="bandpass", fs=fs)
imp = np.zeros(2048)
imp[0] = 1.0
h = signal.lfilter(b, a, imp)                    # la respuesta medida

inv = regularized_inverse_filter(h, fs, f_range=(200.0, 4000.0))
target_db = 20.0 * np.log10(
    np.abs(inv.spectrum[1:]) * np.abs(inv.response_spectrum).max()
)
sweep = shaped_sweep_signal(int(fs), 200.0, 4000.0, 3.0,
                            target=(inv.frequencies[1:], target_db))
sweep.plot(language="es")   # forma de onda y espectro contra el objetivo invertido

Relación con las demás herramientas

El par de Golay se une a sweep_signal y mls_signal en la familia de adquisición ISO 18233; elija el barrido por el rechazo de la distorsión, la MLS por hardware heredado, y el par cuando importe la deconvolución exactamente libre de ruido de un sistema invariante en el tiempo (bancos de HRTF, útiles de calibración). El análisis armónico de lo que los barridos descartan vive en distorsión con seno barrido. La maquinaria de Welch usada para verificar el barrido conformado es la página de análisis espectral calibrado, y las respuestas ecualizadas alimentan la misma cadena posterior que cualquier respuesta al impulso.

Qué cubre esta guía

Cubierto. El par complementario de Golay (Golay 1961; Havelock, Kuwano y Vorländer eds., 2008, parte I, cap. 6 de Xiang): la recursión de adjunción, la identidad exacta de autocorrelaciones complementarias y la recuperación por suma de correlaciones, implementadas por golay_pair y golay_impulse_response. La síntesis de barridos por retardo de grupo de Müller y Massarani (2001, secs. 4.2-4.3), implementada por shaped_sweep_signal. La regularización de Tikhonov dependiente de la frecuencia de Kirkeby y Nelson (1999, Ec. 17 y sec. 2.4), implementada por regularized_inverse_filter.

No cubierto. Son métodos de ingeniería de medición tomados de la literatura de funciones de transferencia, no de una norma de certificación, así que no hay ninguna cláusula de cumplimiento contra la que comprobar. La inversión original de Kirkeby y Nelson apunta a la reproducción sonora multialtavoz, con una regularización matricial para la cancelación de diafonía entre varios canales. regularized_inverse_filter implementa solo el caso escalar de un único canal de su Ec. 17.

Referencias

• Golay, M. J. E. (1961). Complementary series. IRE Transactions on Information Theory, 7(2), 82-87. https://doi.org/10.1109/TIT.1961.1057620La construcción original de los pares complementarios del §1.
• Havelock, D., Kuwano, S. y Vorländer, M. (eds.). (2008). Handbook of signal processing in acoustics. Springer. https://doi.org/10.1007/978-0-387-30441-0Parte I, capítulo 6 (Xiang, Digital Sequences): la recursión de Golay del §1, la identidad de autocorrelaciones complementarias de la Ec. (2) y el procedimiento de recuperación en frecuencia de la Ec. (4) y la Fig. 2. ISBN 978-0-387-77698-9.
• Kirkeby, O. y Nelson, P. A. (1999). Digital filter design for inversion problems in sound reproduction. Journal of the Audio Engineering Society, 47(7/8), 583-595. La regularización de Tikhonov dependiente de la frecuencia del §3 (su Ec. (17)) y el retardo de modelado que hace causal la inversa de fase mixta (Sec. 2.4).
• Müller, S. y Massarani, P. (2001). Transfer-function measurement with sweeps. Journal of the Audio Engineering Society, 49(6), 443-471. El capítulo de síntesis de barridos tras el §2: construcción en el dominio de la frecuencia a partir de magnitud y retardo de grupo (Sec. 4.2), la recursión del retardo de grupo para espectros de magnitud arbitrarios (Sec. 4.3, Ecs. (11)-(12)) y la discusión de inversión limitada en banda de las Secs. 3.1 y 4.5. Se consultó la edición extendida "Director's Cut".