2.25.3.21. pycsamt.ai.processing.tsdenoise#

TimeSeriesDenoiser – MMF-SVM-K-SVD raw field-time-series denoising.

Every other estimator in pycsamt.ai.processing works on the frequency-domain impedance tensor (EMDenoiser) or on station-level tables derived from it. This module works one stage further upstream: on the raw TSData sample stream (Ex, Ey, Hx, Hy, Hz), the same object pycsamt.ts.readers.read_ts() produces and pycsamt.ts.process.ts_to_spectra() FFTs into cross-power spectra. It is meant to run before pycsamt.ts.process.preprocess(), as a strong-interference-suppression pass on the field record itself, complementing rather than duplicating that function’s light gap-fill and the Huber-robust segment weighting already built into cross_spectra().

The method is a direct implementation of Gui et al. (2024) [1] – mathematical-morphological filtering (MMF), a support-vector-machine (SVM) signal/noise classifier, and K-SVD dictionary-learning denoising – extending the authors’ earlier MMF-K-SVD method (Gui et al., 2021) [2] with the SVM triage stage that prevents high-quality low-amplitude MT signal from being treated as noise.

Pipeline#

mmf_split() separates each channel into a low-frequency part (kept untouched – MT signal below ~1 Hz is exactly what a mathematical morphological filter with a sub-second structuring element preserves) and a high-frequency residual that may carry strong interference.
The high-frequency residual is cut into win_seconds windows; SignalQualityClassifier (a linear-kernel SVM trained on four complexity features – sample entropy, fuzzy entropy, approximate entropy, box-counting fractal dimension, see compute_entropy_features()) labels each window clean or noisy, mirroring the paper’s Fig. 4-6.
Contiguous runs of noisy windows are handed to KSVDDenoiser, which learns a small dictionary from overlapping patches of the run itself (self-supervised – no pre-training needed) and subtracts the atoms’ sparse reconstruction, which – because transient interference is far easier to sparsely represent than the genuinely random background signal – approximates the noise profile, not the signal.
Low-frequency part, untouched clean windows, and denoised noisy runs are merged back into one channel.

TimeSeriesDenoiser orchestrates all three stages and accepts / returns a TSData record directly via TimeSeriesDenoiser.apply(), the same TSData-in / TSData-out convention apply() uses for site collections.

References

Functions

`approx_entropy`(x, *[, m, r])	Approximate entropy (Pincus, 1995).
`box_dimension`(x, *[, n_scales])	Box-counting fractal dimension of the `(t, x(t))` graph.
`compute_entropy_features`(segments, *[, m, ...])	Batch feature extraction: (SE, FE, AE, BD) for each row of segments.
`fuzzy_entropy`(x, *[, m, r, n_exp])	Fuzzy entropy (Chen et al., 2007).
`generate_synthetic_library`(win_len, *[, ...])	Build a synthetic clean/noisy training library, mirroring Fig.
`ksvd_dictionary`(Y, n_atoms, sparsity, *[, ...])	K-SVD dictionary learning (Aharon, Elad & Bruckstein, 2006).
`mmf_split`(x, *[, size, dt])	Split a time series into a low-frequency part and a high-frequency residual using mathematical morphological filtering (MMF).
`omp`(D, y, sparsity)	Orthogonal matching pursuit for one signal.
`omp_batch`(D, Y, sparsity)	Orthogonal matching pursuit for a batch of signals (patches).
`sample_entropy`(x, *[, m, r])	Sample entropy (Richman & Moorman, 2000).
`snr_db`(clean, denoised)	Signal-to-noise ratio, Eq.
`time_domain_ncc`(clean, denoised)	Normalized cross-correlation, Eq.

Classes

`KSVDDenoiser`([patch_len, n_atoms, sparsity, ...])	K-SVD dictionary-learning denoiser for a single time-series excursion (typically the noisy windows `SignalQualityClassifier` flagged, concatenated over one contiguous run).
`SignalQualityClassifier`([kernel, C, gamma, ...])	Linear-kernel SVM classifying 10-s time-series windows as high-quality or noisy, from four complexity features.
`TimeSeriesDenoiser`([mmf_size, win_seconds, ...])	MMF-SVM-K-SVD field time-series denoiser (Gui et al., 2024).

pycsamt.ai.processing.tsdenoise.mmf_split(x, *, size=None, dt=None)[source]

Split a time series into a low-frequency part and a high-frequency residual using mathematical morphological filtering (MMF).

Follows Tang et al. (2012) / Li et al. (2020): the low-frequency envelope is the average of the “open-then-close” and “close-then-open” morphological filters with a flat 1-D structuring element of length size,

\[x_{low} = \tfrac{1}{2}\left[ \mathrm{OC}(x) + \mathrm{CO}(x) \right], \qquad \mathrm{OC}(x) = \delta_g(\epsilon_g(x)),\;\; \mathrm{CO}(x) = \epsilon_g(\delta_g(x))\]

where \(\epsilon_g\)/\(\delta_g\) are grayscale erosion/dilation. Averaging the two orderings removes both positive and negative transients symmetrically – a spike that survives opening-then-closing is suppressed by closing-then-opening and vice versa – which is what makes MMF, unlike CEEMD, VMD, or wavelet decomposition, remain reliable at extracting the low-frequency signal even when strong interference is present (see the module docstring and Fig. 7 of Gui et al., 2024).

Parameters:

x (ndarray, shape (n,)) – Input channel. NaN samples are linearly interpolated before filtering and restored as NaN in both outputs.
size (int or None) –
Structuring-element length in samples (odd, >= 3). When None, defaults to the number of samples spanning approximately 2 s (rounded to the nearest odd integer, using dt) – a spike/transient shorter than this is treated as “high frequency”; MT signal below roughly 1 Hz survives as “low frequency”.

This is the parameter that matters most, and the paper itself treats it as data-dependent (Li et al., 2020’s own “detailed principles … on selecting the corresponding structural element sizes”). size must exceed the widest interference pulse the record contains, or that pulse’s plateau survives erosion/dilation and leaks straight into low instead of being rejected into high – e.g. on a real 15 Hz field record with ~3-4 s wide square-wave bursts, the 2 s default left visible spikes in low and only cut noisy-window SNR by a few dB; raising size to ~8 s (121 samples) rejected the bursts cleanly and raised the same window’s SNR by ~17 dB. Inspect low against a plot of x and widen size until the plateaus/steps disappear from it.
dt (float or None) – Sampling interval (s), used only to derive the default size. Required when size is None.

Returns:

low, high – low + high == x (up to floating point). low is the MMF envelope; high is the residual carrying any strong interference plus genuine high-frequency signal.

Return type:

ndarray, shape (n,)

See also

TimeSeriesDenoiser: Uses this as its first stage.

pycsamt.ai.processing.tsdenoise.sample_entropy(x, *, m=2, r=None)[source]

Sample entropy (Richman & Moorman, 2000).

Larger for more irregular, less self-similar signals – the high-quality MT background – and small (near 0) for the stereotyped, self-similar shapes of strong-interference transients (Fig. 4a of Gui et al., 2024).

Parameters:

x (ndarray, shape (n,))
m (int, default 2) – Embedding dimension.
r (float or None) – Tolerance. Defaults to 0.2 * std(x) (Pincus’s convention), the value used throughout the paper.

Returns:

np.nan when the window is too short or degenerate (std(x) == 0) to estimate.

Return type:

float

pycsamt.ai.processing.tsdenoise.fuzzy_entropy(x, *, m=2, r=None, n_exp=2.0)[source]

Fuzzy entropy (Chen et al., 2007).

Replaces sample entropy’s hard d <= r match with a smooth exponential membership \(\exp(-(d/r)^{n})\), making it less sensitive to the tolerance choice; same clean-vs-noise separation as sample_entropy() (Fig. 4b).

Parameters:

x (ndarray, shape (n,))
m (int, default 2)
r (float or None) – Defaults to 0.2 * std(x).
n_exp (float, default 2.0) – Membership-function steepness.

Return type:

float

pycsamt.ai.processing.tsdenoise.approx_entropy(x, *, m=2, r=None)[source]

Approximate entropy (Pincus, 1995).

Same interpretation as sample_entropy() (larger for high-quality signal, Fig. 4c), but counts self-matches, which biases it slightly and makes it less consistent for short windows – kept for parity with the paper’s four-feature vector.

Parameters:

x (ndarray, shape (n,))
m (int, default 2)
r (float or None) – Defaults to 0.2 * std(x).

Return type:

float

pycsamt.ai.processing.tsdenoise.box_dimension(x, *, n_scales=8)[source]

Box-counting fractal dimension of the (t, x(t)) graph.

The time axis and amplitude axis are both normalised to [0, 1] before counting, then a log(box count) vs. log(1/box size) regression slope is fit over dyadic box sizes (\(1/2, 1/4, \dots\)). Larger for the ragged, high-quality MT background than for a smooth or blocky interference shape (Fig. 4d).

Parameters:

x (ndarray, shape (n,))
n_scales (int, default 8) – Number of dyadic box sizes to regress over (capped so the finest box size still spans >= 2 samples).

Return type:

float

pycsamt.ai.processing.tsdenoise.compute_entropy_features(segments, *, m=2, r_frac=0.2, fuzzy_n=2.0)[source]

Batch feature extraction: (SE, FE, AE, BD) for each row of segments.

Parameters:

segments (ndarray, shape (n_segments, win_len))
m (int, default 2) – Embedding dimension shared by SE/FE/AE.
r_frac (float, default 0.2) – Tolerance as a fraction of each segment’s own std.
fuzzy_n (float, default 2.0) – Fuzzy-membership steepness (see fuzzy_entropy()).

Returns:

F – Columns [sample_entropy, fuzzy_entropy, approx_entropy, box_dimension]. Degenerate/too-short rows are NaN.

Return type:

ndarray, shape (n_segments, 4)

pycsamt.ai.processing.tsdenoise.generate_synthetic_library(win_len, *, n_per_class=500, clean_amp=1.0, noise_amp=8.0, seed=None)[source]

Build a synthetic clean/noisy training library, mirroring Fig. 3 of Gui et al. (2024).

Each noisy sample is a weak-amplitude colored-noise background (the same generator as the clean class) with one strong, dominant-amplitude transient superimposed – charge/discharge, square-wave, pulse, or triangle-wave, chosen uniformly at random – matching the “strong interference signal samples have very obvious shapes and amplitudes compared to high-quality signal samples” description in the paper. Used by SignalQualityClassifier.from_synthetic() when no measured labelled library is available.

Parameters:

win_len (int) – Segment length in samples (e.g. round(win_seconds / dt)).
n_per_class (int, default 500) – Number of clean and noisy samples each (paper’s own default).
clean_amp (float, default 1.0) – Background amplitude scale (arbitrary units – only the clean/noisy amplitude ratio matters to the entropy features).
noise_amp (float, default 8.0) – Transient amplitude scale, relative to clean_amp.
seed (int or None)

Returns:

segments (ndarray, shape (2 * n_per_class, win_len))
labels (int ndarray, shape (2 * n_per_class,)) – 1 = high-quality/clean, 0 = noisy.

Return type:

tuple[ndarray, ndarray]

class pycsamt.ai.processing.tsdenoise.SignalQualityClassifier(kernel='linear', C=1.0, gamma='scale', m=2, r_frac=0.2, fuzzy_n=2.0, random_state=None)[source]

Bases: BaseEMProcessor

Linear-kernel SVM classifying 10-s time-series windows as high-quality or noisy, from four complexity features.

Faithful to section 2.2 / 3.1 of Gui et al. (2024): a small, high-margin classifier on [sample_entropy, fuzzy_entropy, approx_entropy, box_dimension] (see compute_entropy_features()) is enough to separate strong interference (low complexity, very consistent shape) from genuine MT signal (high complexity, weak amplitude, strong randomness) – no deep network is needed for a 4-feature, high-margin problem like this one.

Parameters:

kernel (str, default "linear") – Passed to sklearn.svm.SVC – "linear" matches the paper’s own choice (Fig. 5).
C (float, default 1.0)
gamma (str or float, default "scale")
m (int, default 2) – Embedding dimension shared by SE/FE/AE (see compute_entropy_features()).
r_frac (float, default 0.2)
fuzzy_n (float, default 2.0)
random_state (int or None)

Examples

>>> from pycsamt.ai.processing.tsdenoise import SignalQualityClassifier
>>> clf = SignalQualityClassifier.from_synthetic(
...     win_len=150, random_state=0,
... )
>>> clf.predict(segments)

classmethod from_synthetic(win_len, *, n_per_class=500, clean_amp=1.0, noise_amp=8.0, random_state=None, **kwargs)[source]

Build and train from generate_synthetic_library().

Parameters:

win_len (int)
n_per_class (int) – Passed to generate_synthetic_library().
clean_amp (float) – Passed to generate_synthetic_library().
noise_amp (float) – Passed to generate_synthetic_library().
random_state (int or None)
**kwargs – Passed to the constructor.

Return type:

SignalQualityClassifier

fit(X, y, **kwargs)[source]

Train on raw segments.

Parameters:

X (ndarray, shape (n_segments, win_len)) – Raw time-series windows (not pre-computed features).
y (int ndarray, shape (n_segments,)) – 1 = high-quality, 0 = noisy.

Return type:

self

transform(X)[source]

Quality score for each segment.

Parameters:: X (ndarray, shape (n_segments, win_len))
Returns:: scores – Probability of being high-quality, in [0, 1]. NaN rows (degenerate segments) score 0.0 (treated as noisy/untrustworthy).
Return type:: ndarray, shape (n_segments,)

predict(X, *, threshold=0.5)[source]

Predict labels (1=high-quality, 0=noisy).

Parameters:

X (ndarray)
threshold (float)

Return type:

ndarray

pycsamt.ai.processing.tsdenoise.omp(D, y, sparsity)[source]

Orthogonal matching pursuit for one signal.

Parameters:

D (ndarray, shape (L, K)) – Column-normalized dictionary.
y (ndarray, shape (L,))
sparsity (int) – Maximum number of active atoms.

Returns:

code

Return type:

ndarray, shape (K,)

pycsamt.ai.processing.tsdenoise.omp_batch(D, Y, sparsity)[source]

Orthogonal matching pursuit for a batch of signals (patches).

Parameters:

D (ndarray, shape (L, K))
Y (ndarray, shape (L, N))
sparsity (int)

Returns:

codes

Return type:

ndarray, shape (K, N)

pycsamt.ai.processing.tsdenoise.ksvd_dictionary(Y, n_atoms, sparsity, *, n_iter=10, seed=None)[source]

K-SVD dictionary learning (Aharon, Elad & Bruckstein, 2006).

Alternates OMP sparse coding with a per-atom SVD update – for each atom, the rank-1 SVD of the residual restricted to the patches currently using that atom simultaneously updates the atom and its coefficients, which is what distinguishes K-SVD from a plain alternating-least-squares dictionary update (Fig. 2 of Gui et al., 2024). Atoms left unused after a coding pass are replaced by the training patch with the largest current reconstruction error, the standard K-SVD safeguard against dictionary collapse when there are few training patches relative to n_atoms.

Parameters:

Y (ndarray, shape (L, N)) – Training patches, one per column.
n_atoms (int) – Dictionary size K (capped to N if larger).
sparsity (int) – Maximum active atoms per patch (OMP budget).
n_iter (int, default 10)
seed (int or None)

Returns:

D (ndarray, shape (L, K)) – Learned, column-normalized dictionary.
X (ndarray, shape (K, N)) – Final sparse codes.

Return type:

tuple[ndarray, ndarray]

class pycsamt.ai.processing.tsdenoise.KSVDDenoiser(patch_len=32, n_atoms=64, sparsity=4, n_iter=10, overlap=0.5, random_state=None)[source]

Bases: BaseEMProcessor

K-SVD dictionary-learning denoiser for a single time-series excursion (typically the noisy windows SignalQualityClassifier flagged, concatenated over one contiguous run).

Self-supervised: fit() learns the dictionary directly from overlapping patches of the (noisy) input itself – no pre-training or labelled data needed, matching section 2.3 of Gui et al. (2024). Because a short, low-sparsity dictionary can represent a stereotyped transient far better than the genuinely random background signal, the sparse reconstruction of the noisy patches approximates the noise, not the signal, and transform() returns the noisy input minus that reconstruction.

Parameters:

patch_len (int, default 32) – Patch length in samples.
n_atoms (int, default 64) – Dictionary size (capped to the number of extracted patches).
sparsity (int, default 4) – OMP budget per patch.
n_iter (int, default 10) – K-SVD alternations.
overlap (float, default 0.5) – Fractional patch overlap in [0, 1).
random_state (int or None)

Examples

>>> from pycsamt.ai.processing.tsdenoise import KSVDDenoiser
>>> den = KSVDDenoiser(random_state=0).fit(noisy_run)
>>> clean_run = den.transform(noisy_run)

fit(X, **kwargs)[source]

Learn the dictionary from patches of X itself.

Parameters:: X (ndarray, shape (n,))
Return type:: self

noise_profile(X)[source]

Sparse reconstruction of X against the learned dictionary – the extracted noise contour (Fig. 13c of Gui et al., 2024).

Parameters:: X (ndarray)
Return type:: ndarray

transform(X)[source]

Denoise X: X - noise_profile(X).

Parameters:: X (ndarray, shape (n,))
Return type:: ndarray, shape (n,)

class pycsamt.ai.processing.tsdenoise.TimeSeriesDenoiser(mmf_size=None, win_seconds=10.0, quality_clf=None, noise_amp_ratio=8.0, patch_len=32, n_atoms=64, sparsity=4, ksvd_iter=10, ksvd_overlap=0.5, quality_threshold=0.5, random_state=None)[source]

Bases: BaseEMProcessor

MMF-SVM-K-SVD field time-series denoiser (Gui et al., 2024).

Orchestrates mmf_split(), SignalQualityClassifier, and KSVDDenoiser into the three-stage pipeline described in the module docstring, and works directly on a TSData record via apply() – the same place pycsamt.ts.process.preprocess() sits, just upstream of it.

Parameters:

mmf_size (int or None) – Structuring-element length (samples) for mmf_split(). None derives it from each channel’s own dt (~2 s) – the single most important knob to tune: it must exceed the widest interference pulse in the record, see mmf_split()’s docstring for a worked before/after.
win_seconds (float, default 10.0) – SVM classification window length – the paper’s own choice.
quality_clf (SignalQualityClassifier or None) – Pre-trained classifier to reuse (e.g. across stations sharing the same sampling rate and amplitude scale). When None, one is trained automatically the first time transform() needs it, with generate_synthetic_library()’s clean_amp set to a robust (MAD-based) estimate of that channel’s own high-frequency background level – required because raw channel amplitude varies enormously across acquisition systems (raw ADC counts, mV/km, nT, …); training the classifier at a fixed, unrelated amplitude scale leaves every real feature vector far outside the training envelope, a regime where a linear SVM’s decision boundary no longer reflects “more extreme = more noisy” and can misclassify confidently.
noise_amp_ratio (float, default 8.0) – Synthetic interference amplitude as a multiple of the auto-detected (or supplied) clean_amp – matches the roughly 8x contrast Fig. 3 of Gui et al. (2024) shows between their strong-interference and high-quality signal samples.
patch_len (int) – Passed to each per-run KSVDDenoiser.
n_atoms (int) – Passed to each per-run KSVDDenoiser.
sparsity (int) – Passed to each per-run KSVDDenoiser.
ksvd_iter (int) – Passed to each per-run KSVDDenoiser.
ksvd_overlap (float) – Passed to each per-run KSVDDenoiser.
quality_threshold (float, default 0.5) – Windows scoring below this are treated as noisy.
random_state (int or None)

Examples

>>> from pycsamt.ts import read_ts
>>> from pycsamt.ai.processing.tsdenoise import TimeSeriesDenoiser
>>> ts = read_ts("data/MT/TS/kap103as.ts/kap103as.ts")
>>> clean_ts = TimeSeriesDenoiser(random_state=0).apply(ts)

fit(X=None, y=None, *, dt=None, clean_amp=1.0, n_per_class=500, **kwargs)[source]

Train (or attach) the internal SignalQualityClassifier.

Optional – transform()/apply() train one automatically on first use if none is attached, with clean_amp set from the channel’s own data (see quality_clf’s docstring entry above). Call explicitly to control dt / clean_amp / n_per_class directly, or to reuse the same classifier for several apply() calls without retraining each time.

Parameters:

X (ignored) – Present for BaseEMProcessor interface compatibility; training data always comes from generate_synthetic_library().
y (ignored) – Present for BaseEMProcessor interface compatibility; training data always comes from generate_synthetic_library().
dt (float) – Sampling interval (s), needed to size the classification window (round(win_seconds / dt) samples).
clean_amp (float, default 1.0) – Background amplitude scale for the synthetic training library – must be in the same units, and roughly the same order of magnitude, as the channel’s actual high-frequency background level (see quality_clf’s docstring entry). Left at the default 1.0 only makes sense for data already normalised to that scale.
n_per_class (int, default 500) – Passed to generate_synthetic_library().

Return type:

self

transform(X, *, dt, channel='')[source]

Denoise a single 1-D channel array.

Parameters:

X (ndarray, shape (n,))
dt (float) – Sampling interval (s).
channel (str) – Label recorded in diagnostics_ (purely informational).

Return type:

ndarray, shape (n,)

apply(ts, *, channels=None, inplace=False)[source]

Denoise a TSData record.

Parameters:

ts (TSData)
channels (list of str or None) – Channels to denoise. None (default) processes every channel present.
inplace (bool, default False) – When False (the default), ts is left untouched and a denoised copy is returned. When True, ts is modified in place and returned.

Returns:

Same station/metadata, with each requested channel replaced by its denoised reconstruction.

Return type:

TSData

Examples

>>> clean_ts = TimeSeriesDenoiser(
...     random_state=0,
... ).apply(ts)

property diagnostics_: DataFrame[source]

Per-window classification diagnostics from the last apply()/transform() call.

Returns:: Columns: channel, win_start, win_stop (s), score (probability of high-quality), label ("clean"/"noisy").
Return type:: DataFrame

pycsamt.ai.processing.tsdenoise.snr_db(clean, denoised)[source]

Signal-to-noise ratio, Eq. (7) of Gui et al. (2024).

\[SNR = 20 \log_{10} \frac{\|m(s)\|_2}{\|m(s) - g(s)\|_2}\]

Parameters:

clean (ndarray) – Reference (original high-quality) signal \(m(s)\).
denoised (ndarray) – Reconstructed signal \(g(s)\).

Returns:

Decibels; higher is better.

Return type:

float

pycsamt.ai.processing.tsdenoise.time_domain_ncc(clean, denoised)[source]

Normalized cross-correlation, Eq. (8) of Gui et al. (2024).

\[NCC = \frac{\sum_s m(s) g(s)} {\sqrt{\sum_s m(s)^2 \sum_s g(s)^2}}\]

Parameters:

clean (ndarray)
denoised (ndarray)

Returns:

In [-1, 1]; closer to 1 means higher similarity.

Return type:

float