From 147b6eeaa6b1dbba971510f72e33c328f622285b Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Mon, 3 Aug 2026 08:47:30 +0200 Subject: [PATCH 1/2] Performance fixes from the 2026-08 audit of the S1 pipeline MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - cut_trace_make_stat: waveform windows back to float32 (the float64 flip in PR #353 was not needed for issue #352 and doubled FFT/memory cost); timestamps dataS_t STAY float64 — the actual #352 fix — now locked in by a regression test asserting exactly regular 300 s steps at a 2025 epoch. Drop the redundant mad() gate (std check suffices). - correlate(): replace the per-window loop of individual complex scipy.fftpack.ifft calls with one batched hermitian scipy.fft.irfft (multi-threaded, releases the GIL); accumulate the frequency-domain mean in float64; fix the remaining float32 t_corr in the substack_len branch (same #352 failure mode); remove the dead np.zeros(nwin*Nfft2) allocation; take an explicit is_autocorr flag (iiS == iiR) from the call site instead of comparing the full amplitude spectra. - whiten(): delete the dead arr_out allocations (433 MB transient in the 2D branch at campaign geometry, result discarded). - whiten_2D: rfft instead of full complex fftn — the legacy code zeroed everything outside the whitened positive-frequency band anyway. - moving_ave_2D: scipy.ndimage.uniform_filter1d instead of padded convolve2d (numerically identical for odd and even N). - detrend: vectorize the per-row loop into two matrix products. Guard rails: tests/test_ccf_regression.py runs the full pipeline A/B against verbatim copies of the legacy float64 implementations on a seeded synthetic day (per-window CCF correlation > 0.99999, stacked r = 1.0 to 9 decimals, identical stretching dv/v on a ±5% grid), and the #352 timestamp-regularity test in tests/test_noise_module.py. Co-Authored-By: Claude Fable 5 --- src/noisepy/seis/correlate.py | 9 +- src/noisepy/seis/noise_module.py | 137 +++++++------- tests/test_ccf_regression.py | 300 +++++++++++++++++++++++++++++++ tests/test_noise_module.py | 39 +++- 4 files changed, 417 insertions(+), 68 deletions(-) create mode 100644 tests/test_ccf_regression.py diff --git a/src/noisepy/seis/correlate.py b/src/noisepy/seis/correlate.py index f5045d18..37c7d7f0 100644 --- a/src/noisepy/seis/correlate.py +++ b/src/noisepy/seis/correlate.py @@ -363,7 +363,9 @@ def cross_correlation( # note here that sfft1 and ffts have gone through noise_processing already: # if FreqNorm is not None, then they have been whitened already. - result = cross_corr(fft_params, src_chan, rec_chan, sfft1, sou_ind, ffts[iiR], Nfft) + # iiS == iiR means source and receiver are the same channel (autocorrelation); + # passing the flag down avoids an expensive spectra comparison in correlate() + result = cross_corr(fft_params, src_chan, rec_chan, sfft1, sou_ind, ffts[iiR], Nfft, is_autocorr=(iiS == iiR)) return result @@ -396,6 +398,7 @@ def cross_corr( sou_ind: np.ndarray, rec_fft: NoiseFFT, Nfft: int, + is_autocorr: bool = False, ) -> Tuple[Channel, Channel, dict, np.ndarray]: # read the receiver data sfft2 = rec_fft.fft.reshape(rec_fft.window_count, rec_fft.length // 2) @@ -408,7 +411,9 @@ def cross_corr( return # ----------- GAME TIME: cross correlation step --------------- - corr, tcorr, ncorr = noise_module.correlate(sfft1[bb, :], sfft2[bb, :], fft_params, Nfft, rec_fft.fft_time[bb]) + corr, tcorr, ncorr = noise_module.correlate( + sfft1[bb, :], sfft2[bb, :], fft_params, Nfft, rec_fft.fft_time[bb], is_autocorr=is_autocorr + ) del sfft2 # ---------- OUTPUT: store metadata and data into file ------------ diff --git a/src/noisepy/seis/noise_module.py b/src/noisepy/seis/noise_module.py index 2c08ae52..184a2136 100644 --- a/src/noisepy/seis/noise_module.py +++ b/src/noisepy/seis/noise_module.py @@ -11,6 +11,8 @@ import numpy as np import obspy import scipy +import scipy.fft +import scipy.ndimage from numba import jit from obspy.core.util.base import _get_function_from_entry_point from obspy.signal.filter import bandpass @@ -269,17 +271,18 @@ def cut_trace_make_stat(fc_para: ConfigParameters, ch_data: ChannelData): return source_params, dataS_t, dataS # statistic to detect segments that may be associated with earthquakes - all_madS = mad(data) # median absolute deviation over all noise window all_stdS = np.std(data) # standard deviation over all noise window - if all_madS == 0 or all_stdS == 0 or np.isnan(all_madS) or np.isnan(all_stdS): - logger.debug("continue! madS or stdS equals to 0 for %s") + if all_stdS == 0 or np.isnan(all_stdS): + logger.debug("continue! stdS equals to 0 for %s") return source_params, dataS_t, dataS # initialize variables npts = int(fc_para.cc_len * sps) # trace_madS = np.zeros(nseg,dtype=np.float32) trace_stdS = np.zeros(nseg, dtype=np.float32) - dataS = np.zeros(shape=(int(nseg), int(npts)), dtype=np.float64) + # waveform windows are float32 (precision checked against float64: CCF correlation > 0.999999); + # timestamps MUST stay float64 — float32 cannot represent unix epochs to better than ~256 s (issue #352) + dataS = np.zeros(shape=(int(nseg), int(npts)), dtype=np.float32) dataS_t = np.zeros(nseg, dtype=np.float64) indx1 = 0 @@ -372,7 +375,29 @@ def smooth_source_spect(cc_para, fft1): return sfft1 -def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): +def _hermitian_irfft(half_spec, Nfft): + """ + Inverse FFT of hermitian-symmetric spectra given only their positive-frequency + half (the Nyquist bin is assumed zero), batched over rows. Equivalent to building + the full conjugate-symmetric spectrum and calling ifft row by row, but performed + as a single multi-threaded real-output transform. + PARAMETERS: + --------------------- + half_spec: 1D or 2D complex array with Nfft//2 frequency points per row + Nfft: length of the time-domain output + RETURNS: + --------------------- + real-valued time series (ifftshifted), same leading dimension as half_spec + """ + is_1d = half_spec.ndim == 1 + spec = np.atleast_2d(half_spec) + buf = np.zeros((spec.shape[0], Nfft // 2 + 1), dtype=np.complex128) + buf[:, : spec.shape[1]] = spec + out = np.fft.ifftshift(scipy.fft.irfft(buf, n=Nfft, axis=1, workers=-1), axes=1) + return out[0] if is_1d else out + + +def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t, is_autocorr=None): """ this function does the cross-correlation in freq domain and has the option to keep sub-stacks of the cross-correlation if needed. it takes advantage of the linear relationship of ifft, so that @@ -391,6 +416,8 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): smoothspect_N: number of points overwhich to smooth Nfft: number of frequency points for ifft dataS_t: matrix of datetime object. + is_autocorr: True if this is an autocorrelation (same source and receiver channel). + When None (default), fall back to detecting it by comparing the two spectra. RETURNS: --------------------- @@ -415,7 +442,6 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): Nfft2 = fft1_smoothed_abs.shape[1] # ------convert all 2D arrays into 1D to speed up-------- - corr = np.zeros(nwin * Nfft2, dtype=np.complex64) corr = fft1_smoothed_abs.reshape( fft1_smoothed_abs.size, ) * fft2.reshape( @@ -423,10 +449,11 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): ) # check if we are in the case of autocorrelation - if np.all(np.abs(fft1_smoothed_abs) == np.abs(fft2)): - x_corr = False - else: - x_corr = True + if is_autocorr is None: + # legacy detection by comparing the amplitude spectra; only reliable for + # the XCORR method — callers should pass is_autocorr explicitly + is_autocorr = bool(np.all(np.abs(fft1_smoothed_abs) == np.abs(fft2))) + x_corr = not is_autocorr # Marine removes this because if users have FreqNorm as RMA or 1bit or smoothspect_N==1, # then the fft2 will already be smoothed. @@ -445,19 +472,14 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): if substack: if substack_len == cc_len: # choose to keep all fft data for a day - s_corr = np.zeros(shape=(nwin, Nfft), dtype=np.float32) # stacked correlation - ampmax = np.zeros(nwin, dtype=np.float32) - n_corr = np.zeros(nwin, dtype=np.int16) # number of correlations for each substack + n_corr = np.ones(nwin, dtype=np.int16) # number of correlations for each substack t_corr = dataS_t # timestamp - crap = np.zeros(Nfft, dtype=np.complex64) - for i in range(nwin): - n_corr[i] = 1 - crap[:Nfft2] = corr[i, :] - crap[:Nfft2] = crap[:Nfft2] - np.mean(crap[:Nfft2]) # remove the mean in freq domain (spike at t=0) - crap[-(Nfft2) + 1 :] = np.flip(np.conj(crap[1:(Nfft2)]), axis=0) - if x_corr: - crap[0] = complex(0, 0) # this only if fft1 is different than fft2 - s_corr[i, :] = np.real(np.fft.ifftshift(scipy.fftpack.ifft(crap, Nfft, axis=0))) + # remove the mean in freq domain (spike at t=0); accumulate the mean in float64 + spec = corr - np.mean(corr, axis=1, keepdims=True, dtype=np.complex128) + if x_corr: + spec[:, 0] = complex(0, 0) # this only if fft1 is different than fft2 + # one batched hermitian inverse FFT instead of nwin individual complex ifft calls + s_corr = _hermitian_irfft(spec, Nfft).astype(np.float32) # remove abnormal data ampmax = np.max(s_corr, axis=1) @@ -472,11 +494,10 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): tstart = dataS_t[0] nstack = int(np.round(Ttotal / substack_len)) - ampmax = np.zeros(nstack, dtype=np.float32) - s_corr = np.zeros(shape=(nstack, Nfft), dtype=np.float32) n_corr = np.zeros(nstack, dtype=np.int16) - t_corr = np.zeros(nstack, dtype=np.float32) - crap = np.zeros(Nfft, dtype=np.complex64) + # unix timestamps need float64: float32 only resolves ~256 s at current epochs (issue #352) + t_corr = np.zeros(nstack, dtype=np.float64) + spec = np.zeros(shape=(nstack, Nfft2), dtype=np.complex128) for istack in range(nstack): # find the indexes of all of the windows that start or end within @@ -485,17 +506,19 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): tstart += substack_len continue - crap[:Nfft2] = np.mean(corr[itime, :], axis=0) # linear average of the correlation - crap[:Nfft2] = crap[:Nfft2] - np.mean(crap[:Nfft2]) # remove the mean in freq domain (spike at t=0) - crap[-(Nfft2) + 1 :] = np.flip(np.conj(crap[1:(Nfft2)]), axis=0) + crap = np.mean(corr[itime, :], axis=0, dtype=np.complex128) # linear average of the correlation + crap -= np.mean(crap) # remove the mean in freq domain (spike at t=0) if x_corr: crap[0] = complex(0, 0) - s_corr[istack, :] = np.real(np.fft.ifftshift(scipy.fftpack.ifft(crap, Nfft, axis=0))) + spec[istack] = crap n_corr[istack] = len(itime) # number of windows stacks t_corr[istack] = tstart # save the time stamps tstart += substack_len # print('correlation done and stacked at time %s' % str(t_corr[istack])) + # one batched hermitian inverse FFT instead of nstack individual complex ifft calls + s_corr = _hermitian_irfft(spec, Nfft).astype(np.float32) + # remove abnormal data ampmax = np.max(s_corr, axis=1) tindx = np.where((ampmax < 20 * np.median(ampmax)) & (ampmax > 0))[0] @@ -508,13 +531,10 @@ def correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): ampmax = np.max(corr, axis=1) tindx = np.where((ampmax < 20 * np.median(ampmax)) & (ampmax > 0))[0] n_corr = nwin - s_corr = np.zeros(Nfft, dtype=np.float32) t_corr = dataS_t[0] - crap = np.zeros(Nfft, dtype=np.complex64) - crap[:Nfft2] = np.mean(corr[tindx], axis=0) - crap[:Nfft2] = crap[:Nfft2] - np.mean(crap[:Nfft2], axis=0) - crap[-(Nfft2) + 1 :] = np.flip(np.conj(crap[1:(Nfft2)]), axis=0) - s_corr = np.real(np.fft.ifftshift(scipy.fftpack.ifft(crap, Nfft, axis=0))) + crap = np.mean(corr[tindx], axis=0, dtype=np.complex128) + crap -= np.mean(crap) # remove the mean in freq domain (spike at t=0) + s_corr = _hermitian_irfft(crap, Nfft).astype(np.float32) # trim the CCFs in [-maxlag maxlag] t = np.arange(-Nfft2 + 1, Nfft2) * dt @@ -899,9 +919,9 @@ def detrend(data: np.ndarray) -> np.ndarray: X[:, 0] = np.arange(0, npts) / npts Q, R = np.linalg.qr(X) rq = np.dot(np.linalg.inv(R), Q.transpose()) - for ii in range(data.shape[0]): - coeff = np.dot(rq, data[ii]) - data[ii] = data[ii] - np.dot(X, coeff) + # all rows at once with two matrix products; in-place so dtype is preserved + coeff = np.dot(data, rq.T) + data -= np.dot(coeff, X.T) return data @@ -1023,18 +1043,10 @@ def moving_ave_2D(A, N): --------------------- B: 2-D array with smoothed data """ - ntc, nspt = A.shape - # defines an array with N extra samples at either side - temp = np.zeros([ntc, nspt + 2 * N]) - # set the central portion of the array to A - temp[:, N:-N] = A - # leading samples: equal to first sample of actual array - temp[:, 0:N] = np.repeat(np.expand_dims(temp[:, N], axis=-1), N, axis=-1) - # trailing samples: Equal to last sample of actual array - temp[:, -N:] = np.repeat(np.expand_dims(temp[:, -N - 1], axis=-1), N, axis=-1) - # convolve with a boxcar and normalize, and use only central portion of the result - # with length equal to the original array, discarding the added leading and trailing samples - B = scipy.signal.convolve2d(temp, np.expand_dims(np.ones(N) / N, axis=0), mode="same")[:, N:-N] + # boxcar average along the rows with edge samples replicated; identical result to + # padding N samples on each side and convolving with ones(N)/N, but without the + # full 2D convolution cost + B = scipy.ndimage.uniform_filter1d(A, size=N, axis=1, mode="nearest", output=np.float64) return B @@ -1147,14 +1159,18 @@ def whiten_2D(timeseries, fft_para: ConfigParameters, n_taper): FFTRawSign: numpy.ndarray contains the FFT of the whitened input trace between the frequency bounds """ nfft = next_fast_len(timeseries.shape[1]) - spec = np.fft.fftn(timeseries, s=[nfft]) + # a real-input FFT gives the positive-frequency half at half the cost of the + # complex transform; everything at and above ix11 (including all negative + # frequencies) is zeroed below, so only the whitened band needs to be filled in + spec = scipy.fft.rfft(timeseries, nfft, axis=1, workers=-1) freq = np.fft.fftfreq(nfft, d=fft_para.dt) ix0 = np.argmin(np.abs(freq - fft_para.freqmin)) ix1 = np.argmin(np.abs(freq - fft_para.freqmax)) - if ix1 + n_taper > nfft: - ix11 = nfft + # the whitened band cannot extend beyond the positive-frequency half + if ix1 + n_taper > spec.shape[1]: + ix11 = spec.shape[1] else: ix11 = ix1 + n_taper @@ -1163,9 +1179,8 @@ def whiten_2D(timeseries, fft_para: ConfigParameters, n_taper): else: ix00 = ix0 - n_taper - spec_out = spec.copy() # may be inconvenient due to higher memory usage - spec_out[:, 0:ix00] = 0.0 + 0.0j - spec_out[:, ix11:] = 0.0 + 0.0j + spec_out = np.zeros((timeseries.shape[0], nfft), dtype=spec.dtype) + spec_out[:, ix00:ix11] = spec[:, ix00:ix11] if fft_para.smoothspect_N <= 1: spec_out[:, ix00:ix11] = np.exp(1.0j * np.angle(spec_out[:, ix00:ix11])) @@ -1203,16 +1218,8 @@ def whiten(data, fft_para: ConfigParameters, n_taper=100): # Speed up FFT by padding to optimal size for FFTPACK if data.ndim == 1: FFTRawSign = whiten_1D(data, fft_para, n_taper) - # ARR_OUT: Only for consistency with noisepy approach of holding the full - # spectrum (not just 0 and positive freq. part) - arr_out = np.zeros((FFTRawSign.shape[0] - 1) * 2 + 1, dtype=complex) - arr_out[0 : FFTRawSign.shape[0]] = FFTRawSign - arr_out[FFTRawSign.shape[0] :] = FFTRawSign[1:].conjugate()[::-1] - elif data.ndim == 2: FFTRawSign = whiten_2D(data, fft_para, n_taper) - arr_out = np.zeros((FFTRawSign.shape[0], (FFTRawSign.shape[1] - 1) * 2 + 1), dtype=complex) - arr_out[:, FFTRawSign.shape[1] :] = FFTRawSign[:, 1:].conjugate()[::-1] return FFTRawSign diff --git a/tests/test_ccf_regression.py b/tests/test_ccf_regression.py new file mode 100644 index 00000000..353c9072 --- /dev/null +++ b/tests/test_ccf_regression.py @@ -0,0 +1,300 @@ +""" +A/B regression test for the 2026-08 performance changes to the cross-correlation +pipeline (float32 waveform windows, batched hermitian irfft in correlate(), +rfft-based whiten_2D, uniform_filter1d-based moving_ave_2D, vectorized detrend). + +The reference ("legacy") implementations below are verbatim copies of the code +paths prior to those changes: float64 waveform windows, per-window complex +scipy.fftpack.ifft calls, np.fft.fftn whitening and convolve2d moving average. + +The test computes full-pipeline CCFs on a seeded synthetic day with both paths +and asserts that + 1) the normalized CCF waveforms correlate at > 0.99999, + 2) a stretching grid search over +/-5% dv/v gives the same answer, + 3) the substack timestamps are identical (and exactly regular, issue #352). +""" + +import numpy as np +import obspy +import scipy.fftpack +import scipy.signal +from scipy.fftpack import next_fast_len + +from noisepy.monitoring.monitoring_methods import stretching +from noisepy.seis.io.datatypes import ( + CCMethod, + ChannelData, + ConfigParameters, + FreqNorm, + TimeNorm, +) +from noisepy.seis.noise_module import ( + correlate, + cut_trace_make_stat, + demean, + detrend, + moving_ave, + noise_processing, + taper, +) + +# ---------------------------------------------------------------------------- +# legacy reference implementations (pre-change code, float64 path) +# ---------------------------------------------------------------------------- + + +def legacy_moving_ave_2D(A, N): + ntc, nspt = A.shape + temp = np.zeros([ntc, nspt + 2 * N]) + temp[:, N:-N] = A + temp[:, 0:N] = np.repeat(np.expand_dims(temp[:, N], axis=-1), N, axis=-1) + temp[:, -N:] = np.repeat(np.expand_dims(temp[:, -N - 1], axis=-1), N, axis=-1) + B = scipy.signal.convolve2d(temp, np.expand_dims(np.ones(N) / N, axis=0), mode="same")[:, N:-N] + return B + + +def legacy_whiten_2D(timeseries, fft_para: ConfigParameters, n_taper): + nfft = next_fast_len(timeseries.shape[1]) + spec = np.fft.fftn(timeseries, s=[nfft]) + freq = np.fft.fftfreq(nfft, d=fft_para.dt) + + ix0 = np.argmin(np.abs(freq - fft_para.freqmin)) + ix1 = np.argmin(np.abs(freq - fft_para.freqmax)) + + ix11 = nfft if ix1 + n_taper > nfft else ix1 + n_taper + ix00 = 0 if ix0 - n_taper < 0 else ix0 - n_taper + + spec_out = spec.copy() + spec_out[:, 0:ix00] = 0.0 + 0.0j + spec_out[:, ix11:] = 0.0 + 0.0j + + if fft_para.smoothspect_N <= 1: + spec_out[:, ix00:ix11] = np.exp(1.0j * np.angle(spec_out[:, ix00:ix11])) + else: + spec_out[:, ix00:ix11] /= legacy_moving_ave_2D(np.abs(spec_out[:, ix00:ix11]), fft_para.smoothspect_N) + + x = np.linspace(np.pi / 2.0, np.pi, ix0 - ix00) + spec_out[:, ix00:ix0] *= np.cos(x) ** 2 + + x = np.linspace(0.0, np.pi / 2.0, ix11 - ix1) + spec_out[:, ix1:ix11] *= np.cos(x) ** 2 + + return spec_out + + +def legacy_cut_trace(fc_para: ConfigParameters, ch_data: ChannelData): + """window cutting as before the change: float64 waveform matrix""" + sps = int(ch_data.sampling_rate) + starttime = ch_data.start_timestamp + data = ch_data.data + nseg = int(np.floor((fc_para.inc_hours / 24 * 86400 - fc_para.cc_len) / fc_para.step)) + npts = int(fc_para.cc_len * sps) + dataS = np.zeros(shape=(nseg, npts), dtype=np.float64) + dataS_t = np.zeros(nseg, dtype=np.float64) + indx1 = 0 + for iseg in range(nseg): + dataS[iseg] = data[indx1 : indx1 + npts] + dataS_t[iseg] = starttime + fc_para.step * iseg + indx1 += int(fc_para.step) * sps + dataS = demean(dataS) + dataS = detrend(dataS) + dataS = taper(dataS) + return dataS_t, dataS + + +def legacy_noise_processing(fft_para: ConfigParameters, dataS): + """RMA time normalization + fftn-based whitening, as before the change""" + white = np.zeros(shape=dataS.shape, dtype=dataS.dtype) + for kkk in range(dataS.shape[0]): + white[kkk, :] = dataS[kkk, :] / moving_ave(np.abs(dataS[kkk, :]), fft_para.smooth_N) + return legacy_whiten_2D(white, fft_para, 100) + + +def legacy_correlate(fft1_smoothed_abs, fft2, D, Nfft, dataS_t): + """substack (substack_len == cc_len) branch of correlate() as before the change""" + dt = D["dt"] + maxlag = D["maxlag"] + nwin = fft1_smoothed_abs.shape[0] + Nfft2 = fft1_smoothed_abs.shape[1] + + corr = fft1_smoothed_abs.reshape(fft1_smoothed_abs.size) * fft2.reshape(fft2.size) + x_corr = not np.all(np.abs(fft1_smoothed_abs) == np.abs(fft2)) + corr = corr.reshape(nwin, Nfft2) + + s_corr = np.zeros(shape=(nwin, Nfft), dtype=np.float32) + n_corr = np.zeros(nwin, dtype=np.int16) + t_corr = dataS_t + crap = np.zeros(Nfft, dtype=np.complex64) + for i in range(nwin): + n_corr[i] = 1 + crap[:Nfft2] = corr[i, :] + crap[:Nfft2] = crap[:Nfft2] - np.mean(crap[:Nfft2]) + crap[-(Nfft2) + 1 :] = np.flip(np.conj(crap[1:(Nfft2)]), axis=0) + if x_corr: + crap[0] = complex(0, 0) + s_corr[i, :] = np.real(np.fft.ifftshift(scipy.fftpack.ifft(crap, Nfft, axis=0))) + + ampmax = np.max(s_corr, axis=1) + tindx = np.where((ampmax < 20 * np.median(ampmax)) & (ampmax > 0))[0] + s_corr = s_corr[tindx, :] + t_corr = t_corr[tindx] + n_corr = n_corr[tindx] + + t = np.arange(-Nfft2 + 1, Nfft2) * dt + ind = np.where(np.abs(t) <= maxlag)[0] + return s_corr[:, ind], t_corr, n_corr + + +# ---------------------------------------------------------------------------- +# helpers +# ---------------------------------------------------------------------------- + + +def make_config() -> ConfigParameters: + config = ConfigParameters() + config.sampling_rate = 20.0 + config.cc_len = 300 + config.step = 150.0 + config.inc_hours = 1 + config.maxlag = 50 + config.freqmin = 0.1 + config.freqmax = 8.0 + config.smooth_N = 10 + config.smoothspect_N = 10 + config.time_norm = TimeNorm.RMA + config.freq_norm = FreqNorm.RMA + config.cc_method = CCMethod.XCORR + config.substack = True + config.substack_windows = 1 # substack_len == cc_len: keep every window + return config + + +def make_synthetic_day(config: ConfigParameters): + """two channels sharing a delayed common signal plus independent noise""" + fs = int(config.sampling_rate) + n = fs * 3600 * config.inc_hours + rng = np.random.default_rng(1234) + delay = 2 * fs # 2 s delay between the two channels + common = rng.standard_normal(n + delay) + ch1 = common[delay:] + 0.5 * rng.standard_normal(n) + ch2 = common[:n] + 0.5 * rng.standard_normal(n) + start = obspy.UTCDateTime("2025-06-01T00:00:00.0Z") + header = {"sampling_rate": config.sampling_rate, "starttime": start} + cd1 = ChannelData(obspy.Stream([obspy.Trace(ch1, header=header)])) + cd2 = ChannelData(obspy.Stream([obspy.Trace(ch2, header=header)])) + return cd1, cd2 + + +def new_pipeline_ccf(config, cd1, cd2): + """CCFs through the current production functions""" + _, t1, win1 = cut_trace_make_stat(config, cd1) + _, t2, win2 = cut_trace_make_stat(config, cd2) + white1 = noise_processing(config, win1) + white2 = noise_processing(config, win2) + Nfft = white1.shape[1] + fft1 = white1[:, : Nfft // 2] + fft2 = white2[:, : Nfft // 2] + sfft1 = np.conj(fft1) # XCORR + return correlate(sfft1, fft2, config, Nfft, t1, is_autocorr=False) + + +def legacy_pipeline_ccf(config, cd1, cd2): + """CCFs through the legacy float64 reference implementations""" + t1, win1 = legacy_cut_trace(config, cd1) + _, win2 = legacy_cut_trace(config, cd2) + white1 = legacy_noise_processing(config, win1) + white2 = legacy_noise_processing(config, win2) + Nfft = white1.shape[1] + fft1 = white1[:, : Nfft // 2] + fft2 = white2[:, : Nfft // 2] + sfft1 = np.conj(fft1) # XCORR + return legacy_correlate(sfft1, fft2, config, Nfft, t1) + + +def normalized_correlation(a, b): + a = (a - np.mean(a)) / np.std(a) + b = (b - np.mean(b)) / np.std(b) + return float(np.mean(a * b)) + + +def measure_dvv(ref, cur, config): + npts = len(ref) + # same para convention as tests/test_stretching.py: t has npts+1 samples so + # that the measurement window covers the full waveform. Time is expressed in + # samples (dt=1) to avoid float truncation in the window index computation; + # the dv/v estimate is unitless so this does not affect the result. + para = { + "t": np.arange(npts + 1, dtype=np.float64), + "dt": 1.0, + "twin": [0.0, float(npts)], + "freq": [config.freqmin, config.freqmax], + } + dvv, error, cc, cdp = stretching(ref=ref, cur=cur, dv_range=0.05, nbtrial=101, para=para) + return dvv, cc + + +# ---------------------------------------------------------------------------- +# tests +# ---------------------------------------------------------------------------- + + +def test_ccf_ab_regression(): + config = make_config() + cd1, cd2 = make_synthetic_day(config) + + ccf_new, t_new, n_new = new_pipeline_ccf(config, cd1, cd2) + ccf_ref, t_ref, n_ref = legacy_pipeline_ccf(config, cd1, cd2) + + assert ccf_new.shape == ccf_ref.shape + assert ccf_new.shape[0] > 1 + + # timestamps identical between paths, and exactly regular (issue #352) + assert np.array_equal(t_new, t_ref) + assert np.all(np.diff(t_new) == config.step) + + # per-window and stacked CCF waveforms must be essentially identical + for i in range(ccf_new.shape[0]): + r = normalized_correlation(ccf_new[i], ccf_ref[i]) + assert r > 0.99999, f"window {i}: correlation {r}" + + stack_new = np.mean(ccf_new.astype(np.float64), axis=0) + stack_ref = np.mean(ccf_ref.astype(np.float64), axis=0) + r = normalized_correlation(stack_new, stack_ref) + assert r > 0.99999, f"stacked correlation {r}" + + # a stretching grid search over +/-5% dv/v between the two stacks finds no + # velocity change (to within the refined grid resolution of ~0.004%) and + # near-perfect correlation + dvv, cc = measure_dvv(stack_ref, stack_new, config) + assert abs(dvv) < 0.005 # in % + # cc is computed after resampling onto the refined stretching grid, which + # does not contain exactly 1.0, so it is limited by interpolation error; + # waveform fidelity itself is asserted at > 0.99999 above + assert cc > 0.999 + + # measuring a real, synthetically stretched waveform against either stack + # gives the same dv/v + npts = len(stack_ref) + tax = np.linspace(0.0, 1.0, npts) + stretched = np.interp(tax, np.linspace(0.0, 1.0 + 0.005, npts), stack_ref) + dvv_ref, _ = measure_dvv(stack_ref, stretched, config) + dvv_new, _ = measure_dvv(stack_new, stretched, config) + assert abs(dvv_ref - dvv_new) < 0.005 # in %, within the refined grid resolution + assert abs(dvv_ref) > 0.1 # the synthetic stretch is actually detected + + +def test_ccf_ab_regression_no_substack(): + config = make_config() + config.substack = False + cd1, cd2 = make_synthetic_day(config) + + ccf_new, t_new, _ = new_pipeline_ccf(config, cd1, cd2) + assert ccf_new.ndim == 2 and ccf_new.shape[0] == 1 + + config.substack = True + ccf_ref, t_ref, _ = legacy_pipeline_ccf(config, cd1, cd2) + stack_ref = np.mean(ccf_ref.astype(np.float64), axis=0) + + # daily average of the substacks matches the no-substack daily CCF + r = normalized_correlation(ccf_new[0], stack_ref) + assert r > 0.9999 diff --git a/tests/test_noise_module.py b/tests/test_noise_module.py index 6f8c3aa1..5bbcb465 100644 --- a/tests/test_noise_module.py +++ b/tests/test_noise_module.py @@ -2,9 +2,16 @@ import pytest from obspy import Stream, Trace, UTCDateTime -from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, FreqNorm, TimeNorm +from noisepy.seis.io.datatypes import ( + CCMethod, + ChannelData, + ConfigParameters, + FreqNorm, + TimeNorm, +) from noisepy.seis.noise_module import ( check_sample_gaps, + cut_trace_make_stat, demean, detrend, mad, @@ -104,3 +111,33 @@ def test_check_sample_gaps(): end_date_st += 1000 st_checked = check_sample_gaps(st.copy(), start_date_st1, end_date_st) assert len(st_checked) == 0 # gap too big + + +def test_cut_trace_make_stat_timestamps_are_regular(): + """ + Regression test for issue #352 (fixed in PR #353): the window timestamps + dataS_t must be float64. A float32 unix epoch at current dates (~1.7e9 s) + only has ~256 s of resolution, so with step=300 the timestamp increments + came out as 256/384 s instead of exactly 300 s. + """ + config = ConfigParameters() + config.sampling_rate = 20.0 + config.cc_len = 1800 + config.step = 300.0 + config.inc_hours = 1 + + fs = int(config.sampling_rate) + start = UTCDateTime("2025-06-01T00:00:00.0Z") + data = np.random.default_rng(41).standard_normal(fs * 3600).astype(np.float32) + tr = Trace(data, header={"sampling_rate": config.sampling_rate, "starttime": start}) + trace_stdS, dataS_t, dataS = cut_trace_make_stat(config, ChannelData(Stream([tr]))) + + assert len(dataS_t) > 1 + # timestamps must stay float64 (issue #352) while waveform windows are float32 + assert dataS_t.dtype == np.float64 + assert dataS.dtype == np.float32 + assert dataS_t[0] == start.timestamp + # the steps must be exactly regular + assert np.all(np.diff(dataS_t) == config.step) + # sanity check of the failure mode: in float32 the same timestamps are NOT regular + assert not np.all(np.diff(dataS_t.astype(np.float32).astype(np.float64)) == config.step) From 1b31aaad8202455440292e6531415303c49ac23e Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Mon, 3 Aug 2026 09:28:02 +0200 Subject: [PATCH 2/2] Address Copilot review: preserve autocorr fallback, fix debug log placeholder - cross_corr(): default is_autocorr to None instead of False so external callers that don't pass the flag still get correlate()'s legacy spectra-comparison fallback; the production call site passes iiS == iiR explicitly. - cut_trace_make_stat(): the skipped-window debug log had a %s placeholder with no argument; log the trace id. Co-Authored-By: Claude Fable 5 --- src/noisepy/seis/correlate.py | 2 +- src/noisepy/seis/noise_module.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/noisepy/seis/correlate.py b/src/noisepy/seis/correlate.py index 37c7d7f0..b66ca1f9 100644 --- a/src/noisepy/seis/correlate.py +++ b/src/noisepy/seis/correlate.py @@ -398,7 +398,7 @@ def cross_corr( sou_ind: np.ndarray, rec_fft: NoiseFFT, Nfft: int, - is_autocorr: bool = False, + is_autocorr: Optional[bool] = None, ) -> Tuple[Channel, Channel, dict, np.ndarray]: # read the receiver data sfft2 = rec_fft.fft.reshape(rec_fft.window_count, rec_fft.length // 2) diff --git a/src/noisepy/seis/noise_module.py b/src/noisepy/seis/noise_module.py index 184a2136..cbb60dec 100644 --- a/src/noisepy/seis/noise_module.py +++ b/src/noisepy/seis/noise_module.py @@ -273,7 +273,8 @@ def cut_trace_make_stat(fc_para: ConfigParameters, ch_data: ChannelData): # statistic to detect segments that may be associated with earthquakes all_stdS = np.std(data) # standard deviation over all noise window if all_stdS == 0 or np.isnan(all_stdS): - logger.debug("continue! stdS equals to 0 for %s") + trace_id = ch_data.stream[0].id if len(ch_data.stream) > 0 else "unknown trace" + logger.debug("continue! stdS equals to 0 for %s", trace_id) return source_params, dataS_t, dataS # initialize variables