Change Log

May 2026

MAD-gated BIC for polynomial degree selection

The degree search now uses BIC (Bayesian Information Criterion) instead of MAD (Mean Absolute Deviation) for selecting the best polynomial degree. BIC penalizes model complexity, preferring simpler models when fit quality is comparable. To prevent BIC from selecting models with substantially worse median fit, a 20% MAD gate is applied: only candidates whose MAD is within 20% of the best MAD are eligible for BIC comparison. This was motivated by KIP broadband data where pure BIC selected M=2,N=4 (A0=39.86, unphysical) over M=2,N=3 (A0=0.86, correct) — a marginal RSS improvement at the cost of low-frequency accuracy.

Config file argument for run.sh

run.sh now accepts an optional config file path: ./run.sh [--build] [--no-gui] [config.py]. Defaults to setup.py if omitted. The regression test harness uses this to pass per-dataset configs directly, eliminating the previous backup/restore logic.

Try-both-pick-better bandwidth comparison

When the coherent bandwidth covers less than 90% of the fitting range at the configured COHE_LIMIT, Stage 3 now runs the degree search at both the original threshold and a bandwidth-expanded threshold, then compares results by BIC. The expanded fit is only adopted if it produces a strictly lower BIC. This can never cause a regression — at worst it produces the same result as without the rule.

80% Nyquist fitting truncation

Stage 3 now truncates the fitting range to 80% of Nyquist before polynomial fitting. The anti-alias filter rolloff near Nyquist is a digitizer artifact (typically a 100+ dB cliff from FIR decimation) that wastes polynomial degrees on spurious high-frequency poles/zeros. At 100 sps, the cap is 40 Hz — well above typical NORMFREQ values but below the anti-alias rolloff. Applied unconditionally to all channel types. Investigation on RS1D v11 (100 sps velocity) showed that removing the anti-alias contamination improves MAD by 5× (2.1% → 0.4%), eliminates unphysical scale/A0 values, and shifts sensitivity by ~2% toward the true value. The constant FIT_NYQUIST_FRAC = 0.80 is hardcoded as a physical constraint. Full-range arrays are preserved for response output so view3/view4 plots still cover the full spectrum.

Quality grading system (A/B/C/F)

Every run now receives an automatic quality grade based on MAD (fit residual), A0 normalization factor, degree symmetry, and recovery actions (pole reflection, scale correction, coherence fallback). The grade appears in Stage 4’s summary, on the pole/zero plot (view3), and on the response comparison plot (view4). Grade A indicates high confidence; Grade F indicates unreliable results. All grading data is saved in out3.diagnostics.json for automation. See Interpreting the results for threshold details.

R²-gated linear phase removal for flat-response instruments

Stage 3 now detects and removes linear phase delays (time offsets) that prevent stable polynomial fitting for flat-response instruments like accelerometers. The phase is fitted with a linear slope and gated on R²: only when R² > 0.95 (phase dominated by a pure delay) is the delay stripped before polynomial fitting. Velocity sensors with rich non-linear phase from natural resonances (R² < 0.8 typically) are left untouched. This fixes 100% fitting failure on RS6D v3 ENZ accelerometer channels and improves fit quality for other accelerometer datasets.

Coherence-sweep delay estimation

The delay estimator now sweeps descending coherence thresholds (0.99, 0.95, 0.90, 0.80) for selecting phase data points, using the highest-coherence subset first. Previously, a fixed 0.95 threshold was used, which could include noisy lower-coherence points that contaminated the linear fit and suppressed R² below the removal gate. With the sweep, the purest phase signal is always tried first. This fixes cases where R² appeared low at the default threshold but was near-perfect (R² > 0.999) when restricted to the highest-quality data. The log now reports which coherence threshold was used for the delay estimate.

Dual timing display in Stage 4 summary

Stage 4’s summary now reports timing from both sources: Stage 0 (integer-sample cross-correlation on raw waveforms) and Stage 3 (sub-sample linear phase fit on the transfer function). When delay removal is applied, this is also clearly indicated. This gives users full visibility into the timing measurements and what action was taken.

March 2026

setup.py pre-flight validation (Stage 0)

Stage 0 now validates setup.py before any processing begins. It checks that KFILE, UFILE, and KRESP point to existing files, and verifies that numeric parameters (LOWFRQ, NORMFREQ, ZERO_LIMIT, COHE_LIMIT, OVERLAP_FRACTION, NUM_OF_WINDOWS, NMAX, and optionally LP_FILTER_HZ) are within their valid ranges. If any check fails, all errors are reported together and the pipeline halts before reading any data. This catches common misconfiguration mistakes (wrong paths, out-of-range values) immediately rather than mid-run.

Physical stability enforcement in Stage 3

Stage 3 now enforces the physical constraint that all poles must have Re(p) < 0 (left half-plane). Poles with Re(p) > 0 are unstable — they imply exponentially growing impulse response, which is unphysical for any real seismometer. When the nonlinear optimizer converges to an unstable solution, the positive-real-part poles are reflected to the left half-plane (p → −Re(p) + Im(p)·j) and the polynomial is rebuilt from the corrected roots. Additionally, if the polynomial fit converges to a sign-flipped solution (negative leading coefficient, which would produce a negative sensitivity in Stage 4), the scale is forced positive. Both corrections are flagged with a ⚠️ warning in the output. Stage 4 retains a secondary guard and reports ❌ if any unstable poles or negative sensitivity survive.

Stage 3 degree search now reports skipped combinations

When DETERMINE_DEGREE = True, Stage 3 searches polynomial degree combinations M,N up to NMAX. Previously, combinations that yielded unstable poles (after the nonlinear fit) were silently ignored. They are now explicitly reported as ⛔ M=x, N=y → skipped (unstable poles), giving the user full visibility into which configurations were considered and rejected.

Config-parameter annotations on all view plots

All view scripts now overlay relevant configured thresholds and reference frequencies directly on their plots:

  • view0 / view1: coherence plot and transfer function plot show COHE_LIMIT as a red dashed horizontal line and LP_FILTER_HZ (if set) as an orange dotted vertical line.

  • view3: response comparison plot shows NORMFREQ as a labeled vertical line; poles-and-zeros plot shows unstable poles (Re(p) > 0) in orange and zeros that will be snapped to the origin (|z| < ZERO_LIMIT) highlighted in orange, with a corner annotation showing the ZERO_LIMIT value.

  • view4: response comparison plots show NORMFREQ as a labeled vertical line in both amplitude and phase subplots.

``–no-gui`` flag for headless runs

run.sh now accepts a --no-gui flag that forces the Matplotlib Agg backend regardless of whether a display is available. Plots are written to files only and no windows are opened. Useful for batch processing or running over SSH without X11 forwarding. Flags can be combined: ./run.sh --build --no-gui.

Low-pass filter control (``LP_FILTER_HZ``)

A new configuration parameter, LP_FILTER_HZ, has been added to setup.py. It applies a zero-phase Butterworth low-pass filter to both input traces before transfer function estimation (stages 1 and 1e). Because the same filter is applied identically to both traces, it cancels in the transfer function ratio and introduces no phase distortion. This removes aliased energy near Nyquist that can otherwise contaminate coherence estimates. Set to None (default) to disable. Recommended value: 0.7–0.8 × Nyquist.

Timing offset check in Stage 0

Stage 0 now cross-correlates the two input waveforms to estimate any timing offset between the sensors. Results are reported in seconds and samples. A zero-sample lag is ideal; 1–2 samples may be negligible; larger offsets indicate a GPS unlock or timing problem that should be corrected before proceeding.

Raw waveform plot in Stage 0 / view0

Stage 0 now generates raw waveform plots of both input traces (via view0.py) alongside the instrument response comparison plots. This allows visual inspection of the input data before committing to a full pipeline run.

Improved logging across all stages

All stages now emit structured, human-readable output with emoji status indicators (✅ success, ❌ error, ⚠️ warning) and indented step labels. Stage 3 highlights the best-fit polynomial candidate with MAD residuals. Stage 4 reports channel ID, A0 normalization factor, and sensitivity before writing the RESP file.

January 2026

Docker-based Python 3 rewrite

The distribution model has changed from a licensed binary package to an open-source Docker container running Python 3. The old sudo respgen / stage?.exe workflow has been replaced with ./run.sh stage?.py. Dependencies are managed via requirements.txt (ObsPy, NumPy, SciPy, Matplotlib). The configuration file is setup.py (copied from setup.py.example); it is bind-mounted into the container at runtime and is not committed to the repository.

Key changes:

  • All stage scripts rewritten in Python 3 (stage0.py through stage4.py, view0.py through view4.py, view_cohe.py)

  • X11 forwarding enabled for interactive matplotlib plots

  • run.sh auto-mounts data and output directories and forwards signals correctly (--init)

  • Output written to output/ directory on the host