seismic-dsp

Procesamiento de señales sismológicas en el navegador: recortar, filtrar, transformar, diezmar, estimar espectros, integrar y detectar disparos sobre trazas, sin dependencias en runtime.

Alcance

La librería trabaja sobre muestras que ya están en memoria. No incluye lectura de formatos (MiniSEED, SAC, StationXML), clientes FDSN, catálogos, taup, ni manejo de fechas y husos horarios. Todo eso encaja mejor en el servidor: si el servidor ya puede enviar un Float32Array, parsear MiniSEED también en el cliente duplica el trabajo.

El remuestreo se limita a factores enteros. decimate(factor) está incluido; pasar de 100 Hz a 40 Hz queda fuera de la v1, junto con la deconvolución con polos y ceros, rotate (ZNE→ZRT) y assertAligned.

Dos particularidades

Los huecos son parte del modelo. Un hueco no se representa con ceros, porque en sismología el cero es un valor perfectamente válido y la confusión entre ambos es una fuente habitual de errores. Cada función documenta qué hace cuando encuentra un hueco, y cada una de esas declaraciones tiene su propia prueba.

La corrección se verifica contra scipy —y contra obspy en la capa sismológica—. Los resultados se comparan con fixtures generados desde scipy.signal, y los de ruido y disparo desde obspy.signal, con la tolerancia declarada dentro del propio fichero. Cada fixture registra la versión de la biblioteca que lo produjo y el sha256 del script generador, que vive en este repositorio. Los dos ficheros cuyo oráculo es obspy declaran además su versión: un fixture que sobrevive a una actualización de su oráculo ha dejado de serlo.

Notas de uso

butter devuelve secciones de segundo orden (sos) por defecto. A los órdenes habituales en sismología, la forma directa {b, a} pierde precisión por mal condicionamiento del polinomio. {b, a} sigue disponible para órdenes bajos y para contrastar con implementaciones que la usen.

compareToNoiseModel pide la banda de validez de tu corrección de respuesta y devuelve NaN fuera de ella. La razón: con una corrección de sensibilidad escalar, una traza de banda ancha puede quedar por debajo del modelo de ruido bajo alrededor de los 100 s de periodo, porque la respuesta del sensor ya ha caído y el escalar no lo refleja. Recortar el overlay donde deja de ser válido evita mostrar un resultado físicamente imposible.

La FFT acepta cualquier longitud y no la redondea. rfft(x, { n }) con n = 200 o n = 60 001 transforma exactamente esa longitud, como numpy: radix mixto (2, 3, 4, 5), la FFT real como una compleja de media longitud, y Bluestein solo cuando la longitud tiene primos grandes. Lo que no depende de los datos se guarda en un plan por longitud, en una caché de como mucho 16 planes y 32 MiB que no altera ningún resultado. Desde la 0.2.0 un espectrograma de segmentos de 200 muestras es unas 20 veces más rápido que en la 0.1.0, y los resultados coinciden con los de entonces hasta el redondeo (≤ 3·10⁻¹⁵ del pico).

Signal guarda la cantidad física como estructura, no como cadena: { kind, scale, exponent, perHz }. integrate y differentiate la actualizan y lanzan cuando no saben cómo hacerlo. Así, una comprobación que necesite distinguir aceleración de velocidad no depende de que la cadena se haya escrito 'm/s^2' o 'm/s²'.

Estado

La API está implementada y probada de punta a punta, de L0 a L6: 1260 pruebas en 29 ficheros, veintidós fixtures de oráculo y cero dependencias en runtime. Las diecisiete filas del contrato de huecos tienen su prueba — el recuento vive en docs/contrato-de-huecos.md.

| Fase | Contenido | Oráculo | Estado | | ---- | ----------------------------------------- | -------------- | ------ | | L0 | Signal, Span, huecos, recorte, origen | — | ✅ | | L1 | FFT y ventanas | numpy · scipy | ✅ | | L2 | Filtros | scipy | ✅ | | L3 | Tendencia y huecos | scipy | ✅ | | L4 | Espectral | scipy | ✅ | | L5 | Diezmado | scipy · propio | ✅ | | L6 | Ruido y sismología | obspy | ✅ |

Uso

npm install @kuronosu/seismic-dsp

import {
  signal,
  slice,
  effectiveOrigin,
  presentSegments,
  span,
} from '@kuronosu/seismic-dsp'

const s = signal(muestras, 100, {
  gaps: [span(2000, 2600)],
  originOffset: 1735689600,
})

const ventana = slice(s, 6000, 10_000)
effectiveOrigin(ventana) // 1735689660 — aritmética entera, sin error acumulado
presentSegments(ventana) // los tramos con datos

Desarrollo

npm install
npm run check       # formato + typecheck + tests
npm run format      # aplica Prettier
npm run build       # ESM + CJS + .d.ts
npm run fixtures    # regenera los oráculos
npm run readme:tolerances   # reescribe la tabla de tolerancias de este fichero

El formato lo fija Prettier (.prettierrc.json) y lo aplica VS Code al guardar (.vscode/settings.json); npm run check lo verifica, así que el editor y la línea de comandos no pueden discrepar. Los ficheros de fixtures/ quedan fuera: los escribe generate.py y reformatearlos produciría un diff que nadie ha pedido.

Los fixtures se regeneran solo a mano. Actualizar scipy es un paso deliberado: se regeneran los ficheros y el diff se revisa como cualquier otro cambio. Regenerarlos automáticamente en CI haría que dejaran de funcionar como oráculo, porque pasarían a reflejar siempre el comportamiento actual. El cargador comprueba el sha256 del generador y avisa si el fichero ha quedado desfasado respecto al script.

La tolerancia se define en el fixture, no en el test, y su parte absoluta se deriva de la escala de los datos (atol = atolRel · scale). Es una diferencia práctica: al actualizar de scipy 1.17.1 a 1.18.1, sosfiltfilt varió hasta 2,4e-14 · escala, y un atol fijo de 1e-14 habría hecho fallar tres de las ocho combinaciones de referencia sin que hubiera ningún error en la librería.

Versiones del generador: ver fixtures/requirements.lock.

Tolerancias declaradas

Lo que cada fixture promete, tal y como el fichero lo declara. La escribe npm run readme:tolerances leyendo fixtures/*.json, y npm run check se pone rojo si deja de coincidir con ellos: escrita a mano, divergiría.

Aquí están las tolerancias, no los márgenes. Cuánto tendría que empeorar la implementación antes de que un test se ponga rojo se mide contra cada fichero y se justifica en su campo tolerance.note, que es donde vive esa cuenta.

Las pruebas comparan |a − b| <= atol + rtol·|b|, como numpy.allclose, con la tolerancia que sale del fichero y no del test. La perilla es atolRel: el suelo absoluto se DERIVA de la escala de los datos, atol = atolRel · escala, donde escala es max|valores| del propio fichero.

| Fixture | Oráculo | atolRel | escala | atol | rtol | | ------------------ | -------------------------------------------------------------------------------------------------------------- | --------- | ------ | ---------- | ------- | | decimate.json | el cuerpo de scipy.signal.decimate: sosfiltfilt(butter, x)[::q] | 1e-12 | 1.41 | 1.41e-12 | 1e-12 | | design.json | scipy.signal.butter (output='sos' y output='ba') | 1e-13 | 2690 | 2.69e-10 | 1e-12 | | fft.json | numpy.fft.fft / ifft / rfft / irfft, numpy.abs / angle / unwrap | 1e-13 | 147 | 1.47e-11 | 1e-12 | | fftfreq.json | numpy.fft.fftfreq / numpy.fft.rfftfreq | 1e-15 | 125 | 1.25e-13 | 1e-12 | | filtering.json | scipy.signal.lfilter / sosfilt / lfilter_zi / sosfilt_zi | 1e-12 | 1.57 | 1.57e-12 | 1e-12 | | freqz.json | scipy.signal.freqz / scipy.signal.sosfreqz | 1e-10 | 1 | 1e-10 | 1e-12 | | minmax.json | min/max por cubo (NUESTRO: no hay scipy.signal.minmax) | 1e-12 | 62.9 | 6.29e-11 | 1e-12 | | noise.json | obspy.signal.spectral_estimation.get_nlnm/get_nhnm (Peterson 1993, dominio publico) | 1e-12 | 188 | 1.88e-10 | 1e-12 | | padding.json | scipy.signal.filtfilt / sosfiltfilt (padtype y padlen) | 1e-12 | 1.65 | 1.65e-12 | 1e-12 | | periodogram.json | scipy.signal.periodogram (density, spectrum) + scipy.signal.spectrogram(mode='magnitude') para la amplitud | 1e-12 | 1.28 | 1.28e-12 | 1e-12 | | segments.json | scipy.signal.sosfiltfilt / sosfilt, por SEGMENTO | 1e-12 | 1.57 | 1.57e-12 | 1e-12 | | seismic.json | scipy.integrate.cumulative_trapezoid + numpy.gradient/diff, aplicados POR SEGMENTO PRESENTE | 1e-12 | 85.9 | 8.59e-11 | 1e-12 | | sosfiltfilt.json | scipy.signal.sosfiltfilt | 1e-12 | 1.57 | 1.57e-12 | 1e-12 | | spectrogram.json | scipy.signal.spectrogram(mode='psd') | 1e-12 | 0.631 | 6.31e-13 | 1e-12 | | trend.json | scipy.signal.detrend, numpy.median | 1e-12 | 7.88 | 7.88e-12 | 1e-12 | | trigger.json | obspy.signal.trigger.classic_sta_lta_py / recursive_sta_lta_py / trigger_onset, POR SEGMENTO PRESENTE | 1e-12 | 10 | 1e-11 | 1e-12 | | welch.json | scipy.signal.welch (average='mean' y 'median') + scipy.signal.spectrogram(mode='psd') para el promedio parcial | 1e-12 | 1.09 | 1.09e-12 | 1e-12 | | windows.json | scipy.signal.get_window | 1e-14 | 1 | 1e-14 | 1e-12 |

Las cifras van a tres cifras significativas; las exactas están en cada fichero.

Los fixtures sin suelo absoluto derivado

Con atolRel = 0 el suelo vale 0 sea cual sea la escala, así que estos no declaran una escala con sentido. El 0 es deliberado en todos, y es además el valor seguro: pasado a allClose exige igualdad exacta, no la afloja.

• bessel_i0.json (scipy.special.i0) — atolRel = 0, rtol = 1e-13. Solo relativa, y a propósito: I0 crece exponencialmente y el mayor valor de sus 16 casos es 2.93e20, así que un suelo derivado de esa escala taparía cualquier error en vez de acotarlo.
• conditioning.json (scipy.signal.filtfilt vs scipy.signal.sosfiltfilt) — atolRel = 0, rtol = 0. No compara muestras contra muestras, sino una magnitud —la divergencia entre las rutas ba y sos— contra la banda que el generador midió moviendo los coeficientes un ulp. Su perilla real es params.ulpBandFactor = 4, y vive en params porque ese factor en un campo llamado rtol sería una trampa cargada.
• median_bias.json (despejado de scipy.signal.welch(average='median') vs scipy.signal.spectrogram(mode='psd')) — atolRel = 0, rtol = 1e-12. Solo relativa: el factor es un cociente de dos cantidades que scipy calcula en coma flotante, y rtol cubre ese redondeo y nada más.
• padlen.json (scipy.signal.sosfiltfilt / scipy.signal.filtfilt (padlen por defecto)) — atolRel = 0, rtol = 0. Lo que compara son enteros medidos —padLenSos, padLenBa, shortestFilterable—: en un entero no hay nada que redondear, así que la igualdad es exacta.

Los 22 ficheros salen de scipy 1.18.1 · numpy 2.4.4 · Python 3.14.3. noise.json y trigger.json, además, de obspy 1.5.0. Las versiones exactas están en fixtures/requirements.lock.