Make the filters zero-phase again and use padlen

[jsh9]

Fixes #77. One of the four P0 fixes from #114.

Please decide before merging: this changes the GoF scores in Demo_08_GoF_Scores.ipynb, mostly S1 and S2 (see "Effect on GoF scores" below). MATLAB's GoF code filters with zero phase, as this PR does. The S1 and S2 changes come almost entirely from a branch in the Python S_() that MATLAB's S_() does not have.

Problem

_filter_kernel() filtered with scipy.signal.sosfilt(sos, x). That runs the filter once, forward in time, so the output is delayed, and padlen was never used.

• On a symmetric pulse that peaks at 10.00 s, the filtered peak moved to 10.23 s (low-pass), 10.21 s (high-pass), 10.30 s (band-pass) and 10.94 s (band-stop).
• padlen=10 and padlen=1000 gave identical results.
• The docstrings promise zero phase distortion, and describe padlen as the filtfilt() pad length.
• This came in with 976a63a (Goodness-of-Fit Scoring, QOL, and Functional Updates #22), which switched butter() to SOS output and replaced filtfilt(b, a, x, padlen=padlen) with sosfilt(sos, x).

Changes

• PySeismoSoil/helper_signal_processing.py: _filter_kernel() calls scipy.signal.sosfiltfilt(sos, x, padlen=padlen). This keeps the SOS form from Goodness-of-Fit Scoring, QOL, and Functional Updates #22 and restores zero phase and padlen. It affects lowpass(), highpass(), bandpass(), bandstop() and baseline(), and through them Ground_Motion's filter methods, Ground_Motion.baseline_correct() and the GoF scores.
• The 8 padlen docstring links (here and in class_ground_motion.py) now point to sosfiltfilt instead of filtfilt. The meaning of padlen is unchanged.
• padlen=None gives the same pad length as the old filtfilt(b, a) default: 15 for 4th-order low/high-pass, 27 for band-pass/band-stop. As before Goodness-of-Fit Scoring, QOL, and Functional Updates #22, the signal must be longer than the pad length. Ground_Motion's filter methods default to padlen=150. Everything filtered in the tests and notebooks is long enough.
• CHANGELOG.md: a "Fixed" entry.
• The example notebooks were re-run, as CONTRIBUTING.md requires.

Tests

12 new tests, all of which fail on the old code (e.g. assert np.float64(10.23) == 10.0 ± 1.0e-05):

• tests/test_helper_signal_processing.py:
  • test_filter__zero_phase[lowpass|highpass|bandpass|bandstop]: a symmetric pulse keeps its peak time.
  • test_filter__padlen[...]: padlen=10 and padlen=1000 give different outputs. It uses a cosine, which is nonzero at the ends.
• tests/test_class_ground_motion.py::test_filter__zero_phase[...]: the same pulse check through the Ground_Motion methods, with their default padlen=150.

No existing test value changed. The GoF tests use analytic expectations (identical signals, a scaled copy), so the phase does not affect them.

Effect on GoF scores (needs a decision)

The re-run Demo_08_GoF_Scores.ipynb gives different scores:

Record	S1 before → after	S2 before → after	Average before → after
0	-0.044 → 0.165	-0.026 → 0.232	-0.115 → -0.093
1	-0.014 → 0.357	-0.033 → 1.959	-0.069 → 0.190
2	-0.010 → 0.309	-0.012 → 1.952	-0.488 → -0.245
3	-0.084 → 0.264	-0.056 → 1.878	-0.340 → -0.098
4	-0.057 → 0.285	-0.057 → 1.788	-0.333 → -0.137

S3–S7 move by about 0.01–0.07, and S8 of record 0 goes from -0.172 to -0.353.

Why S1 and S2 change: almost all of it comes from the quiet pre-event part of each record, where the measured record's normalized integrals are close to 0.

• The Python S_() takes the absolute difference where the measured value is below eps2 = 1e-12, and the relative difference elsewhere.
• With the causal filter, the measured values there were exactly 0 (clamped to 2.2e-16). So they all got the absolute difference, and scored about 0.
• The backward pass of the zero-phase filter leaks about 1e-11 into them. Values above 1e-12 now get the relative difference, which saturates (about 10) against the simulation's ~1e-6.

MATLAB's S_() has no absolute-difference branch. This is the one in d_1234.m, which computes S1–S4:

function score = S_(meas,simu)
    meas(meas == 0) = eps;
    simu(simu == 0) = eps;
    rela_diff = (simu-meas)./meas;
    score = erf(rela_diff*1)*10;
end

d_567.m and d_89.m have the same function without the two eps lines. The eps2 branch was added to the Python port in 829d39b ("Gof working update", 2024-07-15, on the #22 branch), four months after the MATLAB comparison in 715ce1f.

With MATLAB's S_() in place of the Python one, S1 and S2 are about 2 for every record, with either filter:

Record	S1 before → after	S1 with MATLAB's S_(): causal / zero-phase	S2 before → after	S2 with MATLAB's S_(): causal / zero-phase
0	-0.044 → 0.165	1.975 / 1.973	-0.026 → 0.232	1.994 / 2.005
1	-0.014 → 0.357	2.002 / 1.990	-0.033 → 1.959	1.984 / 1.993
2	-0.010 → 0.309	2.007 / 1.980	-0.012 → 1.952	2.005 / 2.022
3	-0.084 → 0.264	1.932 / 1.923	-0.056 → 1.878	1.961 / 1.939
4	-0.057 → 0.285	1.960 / 1.936	-0.057 → 1.788	1.961 / 1.927

So with MATLAB's formula, the filter moves S1 and S2 by less than 0.04. Neither the old nor the new Python S1 and S2 match MATLAB's formula on these records. The difference is the eps2 branch, which this PR does not change.

MATLAB reference: the MATLAB source of SeismoSoil's GoF functions filters both records with zero phase. SeismoSoil ships these functions as compiled .p files, and the source is not public. gofScores.m computes the scores with d_1234.m, d_567.m and d_89.m. d_1234.m filters like this, and d_567.m the same way:

a_m = baselineCorrHighPass(a_m);
a_s = baselineCorrHighPass(a_s);
...
        a_m = zeroPhaseBandPass(a_m,fmin,fmax,0,filter_order);
        a_s = zeroPhaseBandPass(a_s,fmin,fmax,0,filter_order);
...
    a_m = zeroPhaseHighPass(a_m,fmin);
    a_s = zeroPhaseHighPass(a_s,fmin);

zeroPhaseBandPass.m wraps filtfilt, as do the high-pass, low-pass and band-stop versions:

    [b,a] = butter(filter_order,[freq_lo_norm,freq_hi_norm]);
    filtered_x = filtfilt(b,a,x);

baselineCorrHighPass.m, which d_89.m also calls, filters with the zero-phase high-pass:

a_new = zeroPhaseHighPass([t_cut,a_cut],cutoff_freq);

So this PR makes the Python filters match MATLAB. When the port was compared with MATLAB in 715ce1f ("Edits to make Python gof results match Matlabs within floating pt. error", 2024-03-04), _filter_kernel() already used sosfilt(). It came from a4f7856 ("Update filter kernel to use SOS representation", 2024-02-26), which replaced filtfilt(b, a, x, padlen=padlen) with sosfilt(sos, x). The data for that comparison is not in the repo.

Testing

On the current head, which has main merged in (now with #74, #123 and #115), run locally on Python 3.13:

• python -m pytest tests: 253 passed.
• pre-commit run -a and pydoclint PySeismoSoil: pass.
• No notebook has an error output, and Demo_08_GoF_Scores.ipynb shows the "after" scores above.
• The table with MATLAB's S_() comes from running Demo_08's five records through GOF_Scores.calc_scores(fmin=1.0, fmax=15.0), with helper_gof_scores.S_ replaced by a port of the MATLAB function above, once with each filter. With the Python S_(), the same script reproduces the "before" and "after" columns.

CI passes on this head: tests on Linux, macOS and Windows with Python 3.11–3.13, lint, and the notebook run.

Notes for review

Found while fixing this, and not changed here:

• baseline() delays its output by one sample: the realignment slice a_new[nr_zpad - cross_bound_left - 1:] is off by one. With this fix a pulse peak lands at 10.01 s instead of 10.00 s. That is why baseline_correct() has no zero-phase test here. MATLAB's baselineCorrHighPass.m has the same slice (1-based), so the port matches MATLAB here:

  a_new = a_new(nr_zpad-cross_bound_left+1:end);  % adjusting left end

• The eps2 branch in S_(), described above.

#115 (#60) is merged, and main is merged into this branch. #116 (#79) and #118 (#78) change other library code; #118 changes amplify_motion(), which this PR does not touch. All three re-run the notebooks and add an entry to the "Fixed" group in CHANGELOG.md, so after one of them merges, the others need main merged in again (keep both entries, then re-run the notebooks).

🤖 Generated with Claude Code

https://claude.ai/code/session_01W7fK8T2pPLf3tPCnMCKL1T