Primeros pasos

phonometry es una biblioteca de Python para la medida acústica, desde los bancos de filtros de fracción de octava y la ponderación frecuencial hasta los niveles sonoros y las métricas normalizadas. Cada métrica se verifica frente a la norma que la rige: 536 comprobaciones numéricas en total. Instálala con pip install phonometry; el primer ejemplo de abajo filtra una señal en tercios de octava y devuelve el nivel de presión acústica de cada banda, y el apartado siguiente ancla esos niveles a un tono de calibrador, que es lo que los convierte en una medición.

Instalación

Desde PyPI (recomendada)

pip install phonometry

Extras opcionales

pip install phonometry[plot]    # matplotlib, para gráficas de respuesta y métodos .plot() de resultados
pip install phonometry[perf]    # numba, ponderación temporal 'impulse' más rápida
pip install phonometry[report]  # reportlab + svglib, para que los métodos .report() generen fichas PDF normativas (su panel de gráfica necesita también [plot])
pip install phonometry[full]    # todos los anteriores (recomendado)

Recomiendo pip install phonometry[full]: instala de una vez matplotlib, numba, reportlab y svglib, con lo que quedan habilitadas todas las funcionalidades de la biblioteca. La instalación base calcula todas las métricas solo con NumPy y SciPy; lo único que no queda disponible son las figuras (.plot() y las gráficas de respuesta de los filtros), las fichas PDF normativas (.report()) y el núcleo compilado que acelera la ponderación temporal impulse.

Un matiz sobre [full]: numba declara numpy<2.5, así que phonometry[full] (igual que phonometry[perf]) resuelve NumPy por debajo de 2.5 mientras que una instalación simple obtiene la última versión. numba solo acelera la ponderación temporal impulse, de modo que si prefieres mantener NumPy al día, instala phonometry[plot,report] y deja fuera [perf].

Desde un clon

git clone https://github.com/jmrplens/phonometry.git
cd phonometry
pip install .

Como submódulo de git

git submodule add https://github.com/jmrplens/phonometry.git
# Después, instálalo en modo editable para usarlo desde tu proyecto
pip install -e ./phonometry

La cadena de procesado de un vistazo

Todo análisis con phonometry es un subconjunto de una misma cadena: tomar la señal en bruto, convertirla a unidades físicas, ponderarla en frecuencia, dividirla en bandas normalizadas, suavizarla en el tiempo y reducirla a métricas:

Cadena de procesado de phonometry: señal, calibración, ponderación frecuencial, banco de filtros de octava, ponderación temporal y métricas, con la norma verificada en cada etapa

Léela como una secuencia de decisiones y no como una receta fija. La calibración va primero porque todo lo que viene después es un cociente respecto a 20 µPa, así que un factor de escala equivocado desplaza todos los números posteriores en la misma cantidad y nada de lo que viene detrás puede detectarlo. La ponderación frecuencial la define IEC 61672-1 como un único filtro aplicado a la señal de banda ancha; aplicar en su lugar correcciones de ponderación banda a banda después del banco es un atajo que solo se sostiene para una señal estacionaria. El banco decide la resolución en frecuencia, y la balística y las métricas deciden el comportamiento temporal, actuando las dos sobre la señal ponderada y en ese orden, porque la ponderación temporal exponencial eleva la señal al cuadrado y por tanto no conmuta con nada de lo anterior.

El ejemplo del apartado siguiente recorre Señal → Octava y nada más: sin calibración, así que sus niveles son pascales por hipótesis y no por medición; sin ponderación, así que son Z (sin ponderar); sin balística, así que cada banda es un único valor eficaz sobre todo el registro. «Dar sentido físico a las muestras», más abajo, devuelve la primera etapa a su sitio, y cada etapa es una función o clase independiente que puedes usar por separado; las guías las recorren de izquierda a derecha (Calibración → Ponderación frecuencial → Bancos de filtros → Ponderación temporal → Niveles).

Un primer análisis: bandas de tercio de octava (sin calibrar)

Divide una señal en bandas de tercio de octava y lee el nivel de cada una.

import numpy as np
from phonometry import filters

fs = 48000
t = np.linspace(0, 1, fs, endpoint=False)
# Señal compuesta: 100 Hz + 1000 Hz
signal = np.sin(2 * np.pi * 100 * t) + np.sin(2 * np.pi * 1000 * t)

# Aplicar el banco de filtros de 1/3 de octava
spl, freq = filters.octave_filter(signal, fs=fs, fraction=3)

print(f"Bandas: {freq}")
# Bandas: [12.589254117941678, 15.848931924611138, ..., 19952.623149688785]  (33 bandas)
print(f"SPL [dB]: {spl}")
# SPL [dB]: [46.88395351 47.96774897 49.04991279 ...]  — 90,71 dB a 100 Hz y 90,95 dB a 1 kHz

Las frecuencias de banda salen como números, no como las etiquetas de un analizador. IEC 61260-1 sitúa las frecuencias centrales sobre una rejilla en base 10, $f_{m} = 1000 \cdot 1 0^{x /10 b}$ con $b = 3$ para los tercios de octava, que es por lo que la más baja es 12,589 Hz y la más alta 19 952,6 Hz. Los números redondos que imprime un sonómetro de mano son las designaciones nominales que la misma norma da a esas mismas bandas: pasa nominal=True y freq vuelve como ['12.5', '16', …, '20k'], que es lo que quieres en un eje o en un informe, mientras que los valores exactos son lo que quieres para calcular. Hay 33 bandas porque el rango por defecto va de 12 Hz a 20 kHz; limits=[f_min, f_max] lo estrecha, y estrecharlo suele merecer la pena: consulta Bancos de filtros para la rejilla.

Lee el espectro antes de fiarte de él. Dos de las 33 bandas contienen la señal: 90,71 dB a 100 Hz y 90,95 dB a 1 kHz, ambas a menos de 0,3 dB de 90,97 dB, el nivel de una sinusoide de amplitud unidad leída como pascales ( $20 lo g_{10} (0, 7071/20 μ Pa)$ ). Esa coincidencia es la comprobación de que el banco hace una aritmética defendible; el pequeño déficit en 100 Hz es la respuesta del filtro de banda diezmado en su propio centro, no energía perdida. Todo lo demás de la gráfica es el banco de filtros mirándose a sí mismo: los 63,79 dB a 80 Hz y los 65,13 dB a 125 Hz, y los 53,18 dB a 800 Hz y los 62,10 dB a 1,25 kHz, son los dos tonos colándose por las faldas de los filtros contiguos, de 27 a 29 dB por debajo; la subida lenta de 46,88 dB a 12,5 Hz a 54,71 dB a 50 Hz es la banda atenuada lejana de esos filtros más el suelo numérico. La regla práctica que se deduce: en un espectro de bandas, trata como artefacto del análisis todo lo que quede más de unos 40 dB por debajo de la banda más alta, hasta que una medida de fondo diga lo contrario. Verificación de clase de filtro es la máscara que fija cuánto tienen que caer esas faldas.

Todavía no son niveles físicos. octave_filter devuelve siempre $20 lo g_{10} (p_{rms} /20 μ Pa)$ y, con la calibración por defecto, toma el valor numérico de cada muestra como si fuera la presión en pascales. Una sinusoide de amplitud unidad se lee, por tanto, como 0,7071 Pa eficaces, es decir 90,97 dB: exactamente correcto para una señal sintética definida en pascales, y exactamente incorrecto para una grabación, cuyas muestras valen lo que hayan dado el convertidor y la cadena de ganancia. Hay dos formas honestas de obtener un nivel defendible. Pasa el factor de pascales por unidad digital obtenido con un tono de calibrador, calibration=filters.LevelCalibration(factor=S) (consulta Calibración y dBFS), o, cuando la señal nunca tuvo escala física, pide LevelCalibration(dbfs=True) y lee niveles referidos a fondo de escala digital, que salen negativos y no se pueden confundir con un SPL.

Análisis espectral en tercios de octava de una señal de seis tonos con la PSD en bruto de fondo

Un ejemplo más rico que el fragmento de arriba: seis tonos (20, 100, 500, 2000, 4000 y 15000 Hz) de amplitud 100, así que los niveles de banda llegan a unos 131 dB. La curva gris es la PSD de la señal en bruto, desplazada hacia abajo para compartir el eje.

Mostrar el código de esta figura

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

fs = 48000
t = np.linspace(0, 5, 5 * fs, endpoint=False)
# Seis tonos de amplitud 100, así que las bandas quedan en torno a 130 dB
tones = [20, 100, 500, 2000, 4000, 15000]
y = 100 * np.sum([np.sin(2 * np.pi * f * t) for f in tones], axis=0)
# Aplicar el banco de filtros de 1/3 de octava
spl, freq = filters.octave_filter(y, fs=fs, fraction=3, limits=[12.0, 20000.0])

# Fondo gris: la PSD de la señal en bruto (Welch), desplazada justo por debajo
# de los SPL de banda para comparar ambas formas espectrales en un mismo eje.
f_psd, psd = scipy.signal.welch(y, fs, nperseg=8192)
psd_db = 10 * np.log10(psd + 1e-12)
psd_db += np.max(spl) - np.max(psd_db) - 5

fig, ax = plt.subplots()
ax.semilogx(f_psd, psd_db, color="gray", alpha=0.6, label="PSD de la señal en bruto")
ax.semilogx(freq, spl, marker="o", markerfacecolor="white",
            label="Bandas de 1/3 de octava")
ax.set_xlabel("Frecuencia [Hz]")
ax.set_ylabel("SPL [dB]")
ax.set_xlim(11, 25000)
ax.legend(loc="lower right")
plt.show()

Esa misma llamada puede dibujar además el banco que acaba de diseñar, que es para lo que sirve el extra [plot]: pasa response_plot=filters.ResponsePlot(show=True) y los 33 filtros de banda aparecen en un mismo eje antes de que vuelvan los niveles.

Respuesta en frecuencia del banco de filtros Butterworth de tercio de octava, 33 bandas de 12,5 Hz a 20 kHz con orden 6

El banco por defecto: Butterworth, orden 6, tercio de octava. Cada curva cruza −3 dB en sus propios bordes de banda, que es de donde salen las faldas del párrafo anterior, y la máscara de clase 1 de IEC 61260-1 es la que fija con qué rapidez tienen que caer.

Mostrar el código de esta figura

spl, freq = filters.octave_filter(
    signal, fs=fs, fraction=3, response_plot=filters.ResponsePlot(show=True)
)

La Galería de arquitecturas de filtro tiene esta misma lámina para las otras cuatro arquitecturas.

Dar sentido físico a las muestras

Los niveles de arriba no son una medición, porque nada conectó los números del array con una presión. Esa conexión es un único factor $S$ , la sensibilidad en pascales por unidad digital, y se obtiene grabando un calibrador de nivel conocido a través del mismo micrófono, preamplificador, interfaz y ajuste de ganancia con los que vas a medir. Un calibrador acústico de clase 1 produce 94 dB re 20 µPa a 1 kHz, que son $2 \times 1 0^{- 5} \cdot 1 0^{94/20} = 1, 0024$ Pa eficaces, y no 1 Pa: una diferencia de 0,02 dB que vale la pena tomarse en serio, porque se propaga a todos los niveles posteriores.

# Lo que la cadena de adquisición le hace a la presión. Este micrófono, este
# preamplificador y esta interfaz entregan juntos 0.005 unidades digitales por
# pascal; en una medición real nunca conoces este número, que es justo para lo
# que sirve el calibrador.
chain = 0.005

# Lo que el calibrador escribe en disco: 3 s de 1 kHz a 94 dB re 20 uPa.
p94 = 2e-5 * 10 ** (94.0 / 20)                       # 1.0024 Pa eficaces
calibrator = chain * np.sqrt(2) * p94 * np.sin(
    2 * np.pi * 1000 * np.arange(3 * fs) / fs
)

# ...y la medición: los mismos dos tonos por la misma cadena.
recording = chain * signal

# La sensibilidad es la presión conocida entre el valor eficaz del tono grabado.
cal = p94 / np.sqrt(np.mean(calibrator ** 2))
print(f"{cal:.1f} Pa por unidad digital")
# 200.0 Pa por unidad digital

spl, bands = filters.octave_filter(
    recording, fs=fs, fraction=3, nominal=True,
    calibration=filters.LevelCalibration(factor=cal),
)
print(f"{bands[9]} Hz: {spl[9]:.2f} dB SPL, {bands[19]}: {spl[19]:.2f} dB SPL")
# 100 Hz: 90.71 dB SPL, 1k: 90.95 dB SPL

Hay dos cosas que merece la pena notar. Los niveles calibrados son los mismos números que en el primer análisis, porque allí la señal sintética estaba definida en pascales y esta grabación es ese mismo campo visto a través de una cadena que el calibrador acaba de deshacer; pero esta vez son decibelios re 20 µPa porque los ha puesto ahí una medición, no porque lo hayamos supuesto. Y esa misma grabación leída sin el factor da 44,93 dB en la banda de 1 kHz, 46,0 dB por debajo: un nivel sin calibrar no es aproximadamente correcto, está desviado en lo que haya dado la cadena de ganancia.

En la práctica no se divide a mano. metrology.sensitivity(calibrator, target_spl=94.0, fs=fs) devuelve ese mismo factor y, de paso, comprueba la estabilidad a corto plazo del tono grabado igual que IEC 60942 cualifica al propio calibrador, de modo que un micrófono mal acoplado se detecta aquí en vez de corromper en silencio todos los niveles posteriores. Calibración y dBFS lleva esa llamada, la regla de deriva antes/después y la alternativa digital en dBFS para señales que nunca tuvieron escala física.

Un número en lugar de treinta y tres

Un sonómetro muestra un nivel, no un espectro. Sobre una señal calibrada ese nivel es una línea, porque tanto la ponderación A como el promedio energético están definidos sobre la presión de banda ancha:

weighted = filters.weighting_filter(cal * recording, fs, curve="A")
laeq = 10 * np.log10(np.mean(weighted ** 2) / (2e-5) ** 2)
print(f"LAeq = {laeq:.2f} dB(A)")
# LAeq = 91.03 dB(A)

El nivel sin ponderar de ese mismo segundo es 93,98 dB: la ponderación A quita casi 3 dB, prácticamente todos del tono de 100 Hz, al que la curva atenúa unos 19 dB. signals.laeq(recording, fs, calibration_factor=cal) es la forma de una línea de las dos anteriores y aplica la ponderación internamente, así que no apliques nunca la ponderación A a la señal antes de pasarla a una función de nivel: eso la pondera dos veces. Las envolventes Fast y Slow, $L_{A E}$ , $L_{C p e ak}$ y los niveles percentiles son llamadas de la misma clase: Ponderación temporal y Niveles integrados y estadísticos las tienen, y Construye un sonómetro recorre la cadena entera, del tono de calibrador a todas ellas, en una sola página.

Analizar un archivo de audio

Una medición real son dos grabaciones, no una, y las condiciones en que se hacen deciden lo que valen los niveles:

La misma cadena, sin tocar nada por medio. Graba primero el tono de calibrador y después la medición, con el mismo micrófono, preamplificador, interfaz y ganancia. Fija la ganancia antes del tono de calibrador y no la vuelvas a mover: el factor solo vale mientras la cadena siga como se calibró, y un mando de ganancia tocado es un error mayor que cualquier deriva.
Deja los picos de 10 a 12 dB por debajo del fondo de escala y desactiva todo lo automático. El control automático de ganancia, la supresión de ruido, la limitación y la ecualización de la interfaz cambian todos el nivel en función del nivel, así que una muestra procesada ya no se puede convertir en un nivel de presión acústica. Una recortada tampoco.
Graba a 48 kHz o más. A 44,1 kHz se descartan las bandas por encima de $f_{s} /2$ : el banco devuelve 32 bandas que terminan en 15 848,9 Hz en lugar de las 33 que terminan en 19 952,6 Hz de arriba, y lo avisa con un PhonometryWarning.
Repite la comprobación del calibrador al final. La diferencia entre antes y después es la cota de deriva de todo lo capturado en medio, y un criterio habitual invalida la serie cuando supera 0,5 dB.

Cadena de calibración: calibrador acústico acoplado en el micrófono a 94,0 dB y 1 kHz, preamplificador, interfaz de audio y la sensibilidad que convierte unidades digitales en pascales

La primera de las dos grabaciones: el calibrador acoplado sobre la cápsula a 94,0 dB / 1 kHz, por el mismo preamplificador e interfaz que la medición que viene después, produciendo el factor de pascales por unidad digital del que depende todo lo que sigue.

import numpy as np
from scipy.io import wavfile
from phonometry import filters

# Los dos archivos salen de la misma cadena, en este orden y sin tocar nada.
fs, calibrator = wavfile.read("calibrator.wav")
fs, signal = wavfile.read("measurement.wav")

# wavfile.read devuelve (muestras, canales); el banco espera (canales, muestras)
if calibrator.ndim > 1:
    calibrator = calibrator[:, 0]
if signal.ndim > 1:
    signal = signal[:, 0]        # micrófono de medida en el canal 1

# Sensibilidad en pascales por unidad digital, a partir del tono de calibrador
# de 94 dB. metrology.sensitivity(calibrator, target_spl=94.0, fs=fs) es la
# llamada de la biblioteca, y valida además la estabilidad del tono.
cal = 2e-5 * 10 ** (94.0 / 20) / np.sqrt(np.mean(calibrator.astype(float) ** 2))

spl, freq = filters.octave_filter(
    signal, fs=fs, fraction=3, nominal=True,
    calibration=filters.LevelCalibration(factor=cal),
)

Reduce un archivo multicanal al canal con el que mediste, o traspónlo. La guarda de arriba no es adorno: wavfile.read devuelve un array (muestras, canales) y el banco lee un array 2D como (canales, muestras). Pasa un archivo estéreo tal cual y no salta nada: un array estéreo de 2400 muestras devuelve un spl de forma (2400, 33), un espectro por cada par de muestras, y un archivo de un segundo tarda minutos en producir números sin sentido. Consulta Multicanal y rendimiento para el convenio, y para el caso en que se analizan varios canales calibrados a la vez, cada uno con su propia sensibilidad.

El audio en enteros se convierte de tipo, no se reescala. La salida de wavfile.read se ejecuta sin error porque los enteros se convierten internamente a float64, pero una muestra int16 de +32767 entra en el banco como 32767.0 y no como 1.0. Todos los niveles salen, por tanto, $20 lo g_{10} (32768) = 90, 3$ dB por encima de los de la misma señal leída como float: medido, una sinusoide de 1 kHz de amplitud unidad tiene su máximo en 90,95 dB y esa misma sinusoide escrita en int16 lo tiene en 181,26 dB. Divide por el valor de fondo de escala (signal / 32768.0 para int16) antes de analizar, o deja que lo absorba el paso de calibración: una sensibilidad determinada a partir de una grabación de calibrador en el mismo formato entero ya lleva el factor incorporado. Mezclar los dos, un calibrador leído como float con una medida leída como entero, o al revés, se equivoca en 90 dB sin ningún síntoma. Las reglas están en Calibración y dBFS.

Un número por banda, sobre el archivo entero. octave_filter devuelve un único nivel eficaz por banda para todo el array, quitando antes la media (detrend=True); responde a «cuánta energía hubo en cada banda durante esta grabación», que es la pregunta correcta para una fuente estacionaria y la equivocada para cualquier cosa que cambie. Un paso de vehículo, una máquina que cicla o una sala por la que alguien pasa vuelven como un promedio energético que nadie había pedido. Si el nivel se mueve durante el registro, trocea la señal y llama al banco por trozos, o toma la vía normalizada y usa la balística exponencial Fast/Slow que muestra un sonómetro; mode="peak" da el pico de banda en lugar del promedio, y un pico de banda no es el $L_{C p e ak}$ . Hay un suelo por el otro lado: un filtro de banda necesita varios ciclos de su propia frecuencia central para establecerse, así que un registro de un segundo dice muy poco sobre las bandas de 12,5 o 16 Hz por muchos dígitos que devuelva.

Siguientes pasos

El análisis en octavas de arriba usa phonometry.filters, uno de los dieciocho espacios de nombres de importación. La documentación los agrupa en diez áreas, una por tema de la barra lateral, desde la psicoacústica y la acústica de salas, edificación y vibraciones hasta el ruido ambiental, aeronáutico y submarino, la electroacústica y la simulación de ondas FDTD. Cada objeto de resultado expone una figura .plot(language="en"|"es") de una línea y, cuando una norma define un formato de informe, un método .report() que genera la ficha PDF normativa.

Construir un sonómetro: la misma cadena de principio a fin, del tono de calibrador a $L_{A F} (t)$ , $L_{A e q}$ , $L_{A E}$ , $L_{C p e ak}$ , los niveles percentiles y la comprobación de clase de cada etapa
Todas las guías: el mapa, todas las guías agrupadas por tema con una línea sobre cada una
Glosario: un símbolo, su unidad, el apartado que lo define y la guía que lo calcula
Galería de arquitecturas de filtro: elige una arquitectura e inspecciona sus respuestas
Calibración y dBFS: obtén valores SPL reales
Por qué phonometry: la filosofía de diseño basada en conformidad
Informe de conformidad: el valor esperado y el calculado de las 536 comprobaciones
Referencia de la API: todos los parámetros de cada función
Bibliografía: los libros y artículos que sustentan cada guía, cada uno con un enlace verificado

Creado y mantenido por José Manuel Requena Plens ORCID Todos los proyectos