From ab61e97495bbd5e1ef6fa328d8bc4ba02051662a Mon Sep 17 00:00:00 2001 From: Riviss Date: Fri, 19 Sep 2025 14:30:21 -0700 Subject: [PATCH 1/6] Add pre-smoothing spectral SNR mask option --- sourcespec/ssp_build_spectra.py | 111 +++++++++++++++++++++++++++++--- 1 file changed, 101 insertions(+), 10 deletions(-) diff --git a/sourcespec/ssp_build_spectra.py b/sourcespec/ssp_build_spectra.py index 6eeb591b..2656a9f8 100644 --- a/sourcespec/ssp_build_spectra.py +++ b/sourcespec/ssp_build_spectra.py @@ -31,6 +31,7 @@ from sourcespec.ssp_setup import ssp_exit from sourcespec.ssp_util import ( smooth, cosine_taper, moment_to_mag, MediumProperties) +# NOTE: 'smooth' is a time-domain window convolution; it assumes finite values. from sourcespec.ssp_geom_spreading import ( geom_spread_r_power_n, geom_spread_r_power_n_segmented, geom_spread_boatwright, geom_spread_teleseismic) @@ -472,6 +473,81 @@ def _displacement_to_moment(stats, config): return 4 * math.pi * v3 * rho / (fsa * stats.radiation_pattern) +def _pre_smoothing_snr_mask(config, spec, specnoise): + """ + Suppress low-frequency bins dominated by noise *before* log-frequency smoothing. + This prevents the smoothing kernel from bleeding high noise amplitudes + into the usable band. + Policy: + - Compute per-bin SNR on linear spectra (moment units) after + geometrical spreading & conversion (same stage as in _build_spectrum, + just before smoothing). + - If 'spectral_snr_mask_threshold' <= 0 or missing -> no-op. + - mode 'left_of_first' (default): find the first frequency f_cross where + SNR >= threshold and clamp all lower-frequency amplitudes to at most the + amplitude at f_cross. Optionally apply a gentle ramp (half-cosine) over + 'spectral_snr_mask_ramp_decades' to avoid a hard step. + - mode 'binary': clamp every bin with SNR < threshold to at most the + minimum unmasked amplitude (conservative). + Notes: + - We *do not* insert NaNs: the smoother needs finite values (ssp_util.smooth). + - This runs before _smooth_spectrum(), i.e., before 'make_freq_logspaced()'. + """ + th = getattr(config, 'spectral_snr_mask_threshold', None) + if not th or th <= 0: + return + mode = getattr(config, 'spectral_snr_mask_mode', 'left_of_first') + ramp_dec = float(getattr(config, 'spectral_snr_mask_ramp_decades', 0.0) or 0.0) + + # spec and specnoise share the same linear frequency grid at this stage + freq = spec.freq + s = spec.data + n = specnoise.data + if s.size == 0 or n.size == 0 or freq.size == 0: + return + eps = np.finfo(float).eps + sn = s / np.maximum(n, eps) + + if mode == 'left_of_first': + # Index of first crossing where SNR >= th; if none, do nothing. + idx = np.argmax(sn >= th) + if sn[idx] < th: + # never crosses threshold + spec_id = spec.get_id() + logger.info(f'{spec_id}: pre-smoothing SNR mask enabled ' + f'(threshold={th:g}), but no crossing found: skipped') + return + A0 = s[idx] + f0 = freq[idx] + # Clamp everything left of f0 to at most A0 + left = np.where(freq < f0)[0] + if left.size: + s[left] = np.minimum(s[left], A0) + # Optional soft ramp (half-cosine) across a band of width 'ramp_dec' decades + if ramp_dec > 0.0: + f1 = f0 / (10.0 ** ramp_dec) + band = np.where((freq >= f1) & (freq < f0))[0] + if band.size: + # w goes 0->1 across [f1, f0]; blend clamped (A0) with original + x = (freq[band] - f1) / max(f0 - f1, eps) + w = 0.5 * (1.0 - np.cos(np.pi * x)) + s[band] = np.minimum(s[band], A0 * w + s[band] * (1.0 - w)) + spec.data = s + spec_id = spec.get_id() + logger.info(f'{spec_id}: pre-smoothing SNR mask applied ' + f'(threshold={th:g}, mode={mode}, f_cross={f0:.4f} Hz, ' + f'ramp_decades={ramp_dec:.3f})') + else: + # 'binary' fallback: clamp any low-SNR bin. + mask = sn < th + if np.any(mask): + # Use the minimum amplitude among unmasked bins as a conservative cap. + unmasked = np.where(~mask)[0] + A0 = np.min(s[unmasked]) if unmasked.size else np.min(s) + s[mask] = np.minimum(s[mask], A0) + spec.data = s + + def _smooth_spectrum(spec, smooth_width_decades=0.2): """ Smooth spectrum in a log10-freq space. @@ -497,7 +573,7 @@ def _smooth_spectrum(spec, smooth_width_decades=0.2): spec.make_freq_logspaced(log_df) -def _build_spectrum(config, trace): +def _build_spectrum(config, trace, do_smooth=True): try: spec = Spectrum(obspy_trace=trace) except ValueError as e: @@ -545,13 +621,14 @@ def _build_spectrum(config, trace): f'skipping spectrum\n{str(e)}' ) from e # smooth spectrum. This also creates log-spaced frequencies and data - try: - _smooth_spectrum(spec, config.spectral_smooth_width_decades) - except ValueError as e: - raise RuntimeError( - f'{spec.id}: Error smoothing spectrum: ' - f'skipping spectrum\n{str(e)}' - ) from e + if do_smooth: + try: + _smooth_spectrum(spec, config.spectral_smooth_width_decades) + except ValueError as e: + raise RuntimeError( + f'{spec.id}: Error smoothing spectrum: ' + f'skipping spectrum\n{str(e)}' + ) from e return spec @@ -869,8 +946,22 @@ def _build_signal_and_noise_spectral_streams( for trace_signal in sorted(signal_st, key=lambda tr: tr.id): trace_noise = noise_st.select(id=trace_signal.id)[0] try: - spec = _build_spectrum(config, trace_signal) - specnoise = _build_spectrum(config, trace_noise) + # Optional low-frequency SNR mask *prior to smoothing*: + th = getattr(config, 'spectral_snr_mask_threshold', None) + use_mask = th is not None and th > 0 + if use_mask: + # Build spectra up to (but not including) smoothing + spec = _build_spectrum(config, trace_signal, do_smooth=False) + specnoise = _build_spectrum(config, trace_noise, do_smooth=False) + # Apply the mask before smoothing to avoid bleed into log-smoothing + _pre_smoothing_snr_mask(config, spec, specnoise) + # Now smooth (this creates log-spaced arrays used downstream) + _smooth_spectrum(spec, config.spectral_smooth_width_decades) + _smooth_spectrum(specnoise, config.spectral_smooth_width_decades) + else: + # Default behavior (no mask) + spec = _build_spectrum(config, trace_signal) + specnoise = _build_spectrum(config, trace_noise) _check_spectral_sn_ratio(config, spec, specnoise) except RuntimeError as msg: # RuntimeError is for skipped spectra From 46de9110ad01a284d5c84f58d64d133aa6409199 Mon Sep 17 00:00:00 2001 From: Riviss Date: Mon, 22 Sep 2025 09:16:15 -0700 Subject: [PATCH 2/6] Include Python modules in sdist --- MANIFEST.in | 2 ++ 1 file changed, 2 insertions(+) diff --git a/MANIFEST.in b/MANIFEST.in index 5476c70b..481ead12 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -5,6 +5,8 @@ include CONTRIBUTORS.txt # include only files in docs, not subdirs like _build include docs/* include imgs/* +# Explicitly include all Python modules so new files are picked up +recursive-include sourcespec *.py # exclude configuration_file.rst, which is generated when # building docs exclude docs/configuration_file.rst From 2caa778351d9714e12c23bd07811dcf7dbf1280c Mon Sep 17 00:00:00 2001 From: Riviss Date: Mon, 22 Sep 2025 09:30:06 -0700 Subject: [PATCH 3/6] Handle string spectral SNR mask threshold --- sourcespec/ssp_build_spectra.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/sourcespec/ssp_build_spectra.py b/sourcespec/ssp_build_spectra.py index 2656a9f8..e8b37f0d 100644 --- a/sourcespec/ssp_build_spectra.py +++ b/sourcespec/ssp_build_spectra.py @@ -948,6 +948,16 @@ def _build_signal_and_noise_spectral_streams( try: # Optional low-frequency SNR mask *prior to smoothing*: th = getattr(config, 'spectral_snr_mask_threshold', None) + if isinstance(th, str): + try: + th = float(th) + except ValueError: + logger.warning( + 'Invalid spectral SNR mask threshold %r: disabling mask', + th) + th = None + else: + setattr(config, 'spectral_snr_mask_threshold', th) use_mask = th is not None and th > 0 if use_mask: # Build spectra up to (but not including) smoothing From 9892a56f9f83f4bd8c47c0dd461cd003e6cfea2b Mon Sep 17 00:00:00 2001 From: Riviss Date: Mon, 22 Sep 2025 09:51:15 -0700 Subject: [PATCH 4/6] Visualize SNR mask thresholds on spectra plots --- sourcespec/ssp_build_spectra.py | 22 ++++++++++++ sourcespec/ssp_plot_spectra.py | 63 +++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+) diff --git a/sourcespec/ssp_build_spectra.py b/sourcespec/ssp_build_spectra.py index e8b37f0d..553c2192 100644 --- a/sourcespec/ssp_build_spectra.py +++ b/sourcespec/ssp_build_spectra.py @@ -508,11 +508,19 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): eps = np.finfo(float).eps sn = s / np.maximum(n, eps) + mask_info = { + 'threshold': float(th), + 'mode': mode, + 'ramp_decades': float(ramp_dec), + 'applied': False + } + if mode == 'left_of_first': # Index of first crossing where SNR >= th; if none, do nothing. idx = np.argmax(sn >= th) if sn[idx] < th: # never crosses threshold + spec.stats.snr_mask_info = mask_info spec_id = spec.get_id() logger.info(f'{spec_id}: pre-smoothing SNR mask enabled ' f'(threshold={th:g}), but no crossing found: skipped') @@ -533,6 +541,13 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): w = 0.5 * (1.0 - np.cos(np.pi * x)) s[band] = np.minimum(s[band], A0 * w + s[band] * (1.0 - w)) spec.data = s + mask_info.update({ + 'applied': True, + 'f_cross': float(f0), + }) + if ramp_dec > 0.0: + mask_info['f_ramp_start'] = float(f1) + spec.stats.snr_mask_info = mask_info spec_id = spec.get_id() logger.info(f'{spec_id}: pre-smoothing SNR mask applied ' f'(threshold={th:g}, mode={mode}, f_cross={f0:.4f} Hz, ' @@ -546,6 +561,13 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): A0 = np.min(s[unmasked]) if unmasked.size else np.min(s) s[mask] = np.minimum(s[mask], A0) spec.data = s + mask_info.update({ + 'applied': True, + 'freq_min': float(np.min(freq[mask])), + 'freq_max': float(np.max(freq[mask])), + 'masked_bins': int(np.count_nonzero(mask)), + }) + spec.stats.snr_mask_info = mask_info def _smooth_spectrum(spec, smooth_width_decades=0.2): diff --git a/sourcespec/ssp_plot_spectra.py b/sourcespec/ssp_plot_spectra.py index 389c0208..1eb29875 100644 --- a/sourcespec/ssp_plot_spectra.py +++ b/sourcespec/ssp_plot_spectra.py @@ -472,6 +472,46 @@ def _snratio_text(spec, ax, color, path_effects): SNRATIO_TEXT_YPOS -= 0.05 +def _snr_mask_text(spec, ax, path_effects): + global SNRATIO_TEXT_YPOS # pylint: disable=global-statement + mask_info = getattr(spec.stats, 'snr_mask_info', None) + if not mask_info: + return + threshold = mask_info.get('threshold') + if threshold is None: + return + mode = mask_info.get('mode', 'left_of_first') + applied = mask_info.get('applied', False) + text = f'Mask ≥{threshold:g}' + if mode == 'left_of_first': + if applied and mask_info.get('f_cross'): + f_cross = mask_info['f_cross'] + text += f' @ {f_cross:.3g} Hz' + ramp_dec = mask_info.get('ramp_decades', 0.0) + f_ramp_start = mask_info.get('f_ramp_start') + if ramp_dec and f_ramp_start: + text += f' (ramp ≥{f_ramp_start:.3g} Hz)' + else: + text += ' (no crossing)' + else: + text += ' (binary)' + if applied: + f_min = mask_info.get('freq_min') + f_max = mask_info.get('freq_max') + if f_min is not None and f_max is not None: + if math.isclose(f_min, f_max): + text += f' ≤{f_max:.3g} Hz' + else: + text += f' {f_min:.3g}-{f_max:.3g} Hz' + else: + text += ', none masked' + ax.text( + 0.95, SNRATIO_TEXT_YPOS, text, ha='right', va='top', + fontsize=8, color='dimgray', path_effects=path_effects, + transform=ax.transAxes, zorder=20) + SNRATIO_TEXT_YPOS -= 0.05 + + def _station_text(spec, ax, color, path_effects, stack_plots): station_text = f'{spec.id[:-1]} {spec.stats.instrtype}' if not stack_plots: @@ -595,6 +635,25 @@ def _plot_fc_and_mw(spec, ax, ax2): ax.axvline(fc, color='#999999', linewidth=2., zorder=1) +def _plot_snr_mask_indicator(spec, ax): + mask_info = getattr(spec.stats, 'snr_mask_info', None) + if not mask_info: + return + if mask_info.get('mode', 'left_of_first') != 'left_of_first': + return + if not mask_info.get('applied'): + return + f_cross = mask_info.get('f_cross') + if not f_cross: + return + ax.axvline( + f_cross, color='#CC9933', linewidth=1.5, linestyle='--', zorder=5) + f_ramp_start = mask_info.get('f_ramp_start') + if f_ramp_start and f_ramp_start < f_cross: + ax.axvspan( + f_ramp_start, f_cross, color='#CC9933', alpha=0.12, zorder=2) + + def _plot_spec(config, plot_params, spec, spec_noise): """Plot one spectrum (and its associated noise).""" plotn = plot_params.plotn @@ -623,9 +682,13 @@ def _plot_spec(config, plot_params, spec, spec_noise): # Write spectral S/N for regular Z,N,E components if orientation not in special_orientations: _snratio_text(spec, ax, color, path_effects) + _snr_mask_text(spec, ax, path_effects) + else: + _snr_mask_text(spec, ax, path_effects) # Plot fc and Mw if a synthetic spectrum S is available if orientation == 'S': _plot_fc_and_mw(spec, ax, ax2) + _plot_snr_mask_indicator(spec, ax) elif plot_type == 'weight': ax.semilogx( spec.freq, spec.data, color=color, alpha=alpha, zorder=20) From 476bd97dbe01f299b5c8af3ef089eb4c3eebe40c Mon Sep 17 00:00:00 2001 From: Riviss Date: Mon, 22 Sep 2025 10:20:23 -0700 Subject: [PATCH 5/6] Apply SNR mask using RSS spectrum and plot event moment --- sourcespec/ssp_build_spectra.py | 168 ++++++++++++++++++++++++++------ sourcespec/ssp_plot_spectra.py | 25 +++++ 2 files changed, 161 insertions(+), 32 deletions(-) diff --git a/sourcespec/ssp_build_spectra.py b/sourcespec/ssp_build_spectra.py index 553c2192..056f237b 100644 --- a/sourcespec/ssp_build_spectra.py +++ b/sourcespec/ssp_build_spectra.py @@ -570,6 +570,102 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): spec.stats.snr_mask_info = mask_info +def _apply_mask_info_to_component(spec, mask_info): + """Apply RSS-derived mask information to a single component spectrum.""" + if not mask_info: + return + # Store a copy of the mask info so later consumers (plots) can reuse it. + spec.stats.snr_mask_info = dict(mask_info) + if not mask_info.get('applied'): + return + mode = mask_info.get('mode', 'left_of_first') + freq = spec.freq + data = spec.data + if mode == 'left_of_first': + f_cross = mask_info.get('f_cross') + if not f_cross: + return + ramp_dec = float(mask_info.get('ramp_decades', 0.0) or 0.0) + eps = np.finfo(float).eps + idx_candidates = np.where(freq >= f_cross)[0] + if not idx_candidates.size: + return + idx = idx_candidates[0] + A0 = data[idx] + left = np.where(freq < f_cross)[0] + if left.size: + data[left] = np.minimum(data[left], A0) + if ramp_dec > 0.0: + f_ramp_start = mask_info.get('f_ramp_start') + if f_ramp_start is None: + f_ramp_start = f_cross / (10.0 ** ramp_dec) + if f_ramp_start < f_cross: + band = np.where((freq >= f_ramp_start) & (freq < f_cross))[0] + if band.size: + x = (freq[band] - f_ramp_start) / max(f_cross - f_ramp_start, eps) + w = 0.5 * (1.0 - np.cos(np.pi * x)) + data[band] = np.minimum(data[band], A0 * w + data[band] * (1.0 - w)) + spec.data = data + else: + freq_min = mask_info.get('freq_min') + freq_max = mask_info.get('freq_max') + if freq_min is None or freq_max is None: + return + mask = (freq >= freq_min) & (freq <= freq_max) + if np.any(mask): + unmasked = np.where(~mask)[0] + A0 = np.min(data[unmasked]) if unmasked.size else np.min(data) + data[mask] = np.minimum(data[mask], A0) + spec.data = data + + +def _apply_rss_snr_mask(config, spec_st, specnoise_st): + """Apply the spectral SNR mask using the root sum of squares spectrum.""" + if not spec_st: + return + spec_ids = {sp.id[:-1] for sp in spec_st} + vertical_codes = config.vertical_channel_codes + wave_type = config.wave_type + for specid in spec_ids: + spec_sel = _select_spectra(spec_st, specid) + if specnoise_st: + specnoise_sel = _select_spectra(specnoise_st, specid) + else: + specnoise_sel = SpectrumStream() + codes = {sp.stats.channel[:-1] for sp in spec_sel} + for code in codes: + components = SpectrumStream( + sp for sp in spec_sel + if sp.stats.channel[:-1] == code and + sp.stats.channel[-1] not in ['H', 'h', 'S', 's', 't'] + ) + if not components: + continue + noises = SpectrumStream( + sp for sp in specnoise_sel + if sp.stats.channel[:-1] == code and + sp.stats.channel[-1] not in ['H', 'h', 'S', 's', 't'] + ) + spec_h = _compute_spec_h( + components, code, vertical_codes, wave_type) + if spec_h is None or getattr(spec_h.stats, 'ignore', False): + continue + specnoise_h = _compute_spec_h( + noises, code, vertical_codes, wave_type) + if specnoise_h is None: + continue + _pre_smoothing_snr_mask(config, spec_h, specnoise_h) + mask_info = getattr(spec_h.stats, 'snr_mask_info', None) + if not mask_info: + continue + mask_info = dict(mask_info) + mask_info['source_channel'] = spec_h.stats.channel + for component in components: + if getattr(component.stats, 'ignore', False): + continue + _apply_mask_info_to_component(component, mask_info) + + def _smooth_spectrum(spec, smooth_width_decades=0.2): """ Smooth spectrum in a log10-freq space. @@ -965,49 +1061,57 @@ def _build_signal_and_noise_spectral_streams( """ spec_st = SpectrumStream() specnoise_st = SpectrumStream() + spec_items = [] + + th = getattr(config, 'spectral_snr_mask_threshold', None) + if isinstance(th, str): + try: + th = float(th) + except ValueError: + logger.warning( + 'Invalid spectral SNR mask threshold %r: disabling mask', + th) + th = None + else: + setattr(config, 'spectral_snr_mask_threshold', th) + use_mask = th is not None and th > 0 + for trace_signal in sorted(signal_st, key=lambda tr: tr.id): trace_noise = noise_st.select(id=trace_signal.id)[0] try: - # Optional low-frequency SNR mask *prior to smoothing*: - th = getattr(config, 'spectral_snr_mask_threshold', None) - if isinstance(th, str): - try: - th = float(th) - except ValueError: - logger.warning( - 'Invalid spectral SNR mask threshold %r: disabling mask', - th) - th = None - else: - setattr(config, 'spectral_snr_mask_threshold', th) - use_mask = th is not None and th > 0 - if use_mask: - # Build spectra up to (but not including) smoothing - spec = _build_spectrum(config, trace_signal, do_smooth=False) - specnoise = _build_spectrum(config, trace_noise, do_smooth=False) - # Apply the mask before smoothing to avoid bleed into log-smoothing - _pre_smoothing_snr_mask(config, spec, specnoise) - # Now smooth (this creates log-spaced arrays used downstream) - _smooth_spectrum(spec, config.spectral_smooth_width_decades) - _smooth_spectrum(specnoise, config.spectral_smooth_width_decades) - else: - # Default behavior (no mask) - spec = _build_spectrum(config, trace_signal) - specnoise = _build_spectrum(config, trace_noise) - _check_spectral_sn_ratio(config, spec, specnoise) + spec = _build_spectrum( + config, trace_signal, do_smooth=not use_mask) + specnoise = _build_spectrum( + config, trace_noise, do_smooth=not use_mask) except RuntimeError as msg: - # RuntimeError is for skipped spectra logger.warning(msg) continue - except SpectrumIgnored as msg: - _ignore_spectrum(msg, spec, specnoise) - trace_original = original_st.select(id=trace_signal.id)[0] - _ignore_trace(msg, trace_original) + spec_items.append((spec, specnoise, trace_signal)) spec_st.append(spec) specnoise_st.append(specnoise) + if not spec_st: logger.error('No spectra left! Exiting.') ssp_exit() + + if use_mask: + # Apply the mask on root-sum-of-squares spectra before smoothing. + _apply_rss_snr_mask(config, spec_st, specnoise_st) + for spec in spec_st: + _smooth_spectrum(spec, config.spectral_smooth_width_decades) + for specnoise in specnoise_st: + _smooth_spectrum(specnoise, config.spectral_smooth_width_decades) + + for spec, specnoise, trace_signal in spec_items: + try: + _check_spectral_sn_ratio(config, spec, specnoise) + except SpectrumIgnored as msg: + _ignore_spectrum(msg, spec, specnoise) + trace_original = original_st.select(id=trace_signal.id)[0] + _ignore_trace(msg, trace_original) + except RuntimeError as msg: + logger.warning(msg) + # build H component _build_H( spec_st, specnoise_st, config.vertical_channel_codes, config.wave_type) diff --git a/sourcespec/ssp_plot_spectra.py b/sourcespec/ssp_plot_spectra.py index 1eb29875..e929307c 100644 --- a/sourcespec/ssp_plot_spectra.py +++ b/sourcespec/ssp_plot_spectra.py @@ -654,6 +654,30 @@ def _plot_snr_mask_indicator(spec, ax): f_ramp_start, f_cross, color='#CC9933', alpha=0.12, zorder=2) +def _plot_event_moment_line(spec, ax, ax2): + """Plot the catalog/event seismic moment as a reference line.""" + if getattr(ax, '_event_moment_drawn', False): + return + event = getattr(spec.stats, 'event', None) + if not event: + return + scalar_moment = getattr(getattr(event, 'scalar_moment', None), 'value', None) + if scalar_moment is None or scalar_moment <= 0: + return + color = '#555555' + ax.axhline( + scalar_moment, color=color, linestyle='--', linewidth=1.5, zorder=4) + ax._event_moment_drawn = True # pylint: disable=protected-access + if ax2 and not getattr(ax2, '_event_moment_drawn', False): + try: + Mw = float(moment_to_mag(scalar_moment)) + except (TypeError, ValueError): + return + ax2.axhline( + Mw, color=color, linestyle='--', linewidth=1.0, zorder=4) + ax2._event_moment_drawn = True # pylint: disable=protected-access + + def _plot_spec(config, plot_params, spec, spec_noise): """Plot one spectrum (and its associated noise).""" plotn = plot_params.plotn @@ -689,6 +713,7 @@ def _plot_spec(config, plot_params, spec, spec_noise): if orientation == 'S': _plot_fc_and_mw(spec, ax, ax2) _plot_snr_mask_indicator(spec, ax) + _plot_event_moment_line(spec, ax, ax2) elif plot_type == 'weight': ax.semilogx( spec.freq, spec.data, color=color, alpha=alpha, zorder=20) From 85f67ee8a7e22632ca01df1250eb47fc38a382a7 Mon Sep 17 00:00:00 2001 From: Riviss Date: Mon, 22 Sep 2025 11:33:09 -0700 Subject: [PATCH 6/6] Add configurable RSS mask frequency range --- sourcespec/config_files/configspec.conf | 6 ++ sourcespec/ssp_build_spectra.py | 81 ++++++++++++++++++++++--- 2 files changed, 78 insertions(+), 9 deletions(-) diff --git a/sourcespec/config_files/configspec.conf b/sourcespec/config_files/configspec.conf index 839347ea..d557d847 100644 --- a/sourcespec/config_files/configspec.conf +++ b/sourcespec/config_files/configspec.conf @@ -271,6 +271,12 @@ spectral_win_length = float(min=1e-99, default=None) # Default value of 0.2 is generally a good choice spectral_smooth_width_decades = float(min=1e-99, default=0.2) +# Frequency range (Hz) to search for the spectral S/N mask threshold. +# Provide two comma-separated values (min_freq, max_freq). The maximum +# value is clipped to 1 Hz inside the code. Leave this parameter unset to +# use the lowest available frequency up to 1 Hz. +spectral_snr_mask_freq_range = float_list(min=0, default=None) + # Residuals file path # An HDF5 file with the mean residuals per station, used for station # correction. This file is generally created using the command diff --git a/sourcespec/ssp_build_spectra.py b/sourcespec/ssp_build_spectra.py index 056f237b..07dc7aff 100644 --- a/sourcespec/ssp_build_spectra.py +++ b/sourcespec/ssp_build_spectra.py @@ -515,16 +515,71 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): 'applied': False } + freq_range_cfg = getattr(config, 'spectral_snr_mask_freq_range', None) + freq_range = None + range_idx = None + if freq_range_cfg: + try: + if isinstance(freq_range_cfg, (list, tuple)): + values = list(freq_range_cfg) + else: + cleaned = str(freq_range_cfg).replace('[', '').replace(']', '') + values = [ + float(val.strip()) for val in cleaned.split(',') if val.strip() + ] + if not values: + raise ValueError + if len(values) == 1: + f_range_min = float(values[0]) + f_range_max = 1.0 + else: + f_range_min = float(values[0]) + f_range_max = float(values[1]) + except (TypeError, ValueError): + logger.warning( + 'Invalid spectral SNR mask frequency range %r: ignoring', + freq_range_cfg) + else: + if f_range_max < f_range_min: + f_range_min, f_range_max = f_range_max, f_range_min + f_range_min = max(f_range_min, freq[0]) + f_range_max = min(max(f_range_min, f_range_max), 1.0) + if f_range_min < f_range_max: + # reduce to available frequency support + idx = np.where((freq >= f_range_min) & (freq <= f_range_max))[0] + if idx.size: + range_idx = idx + # store effective range for logging/plots + freq_range = (float(freq[idx[0]]), float(freq[idx[-1]])) + mask_info.update({ + 'freq_range_min': freq_range[0], + 'freq_range_max': freq_range[1] + }) + else: + logger.info( + '%s: spectral SNR mask frequency range [%g, %g] Hz ' + 'does not intersect spectrum: ignoring', + spec.get_id(), f_range_min, f_range_max) + else: + logger.warning( + 'Invalid spectral SNR mask frequency range %r: ignoring', + freq_range_cfg) + if mode == 'left_of_first': - # Index of first crossing where SNR >= th; if none, do nothing. - idx = np.argmax(sn >= th) - if sn[idx] < th: - # never crosses threshold - spec.stats.snr_mask_info = mask_info - spec_id = spec.get_id() - logger.info(f'{spec_id}: pre-smoothing SNR mask enabled ' - f'(threshold={th:g}), but no crossing found: skipped') - return + range_forced = range_idx is not None and np.any(sn[range_idx] < th) + if range_forced: + idx = range_idx[-1] + else: + # Index of first crossing where SNR >= th; if none, do nothing. + idx = np.argmax(sn >= th) + if sn[idx] < th: + # never crosses threshold + spec.stats.snr_mask_info = mask_info + spec_id = spec.get_id() + logger.info( + f'{spec_id}: pre-smoothing SNR mask enabled ' + f'(threshold={th:g}), but no crossing found: skipped') + return A0 = s[idx] f0 = freq[idx] # Clamp everything left of f0 to at most A0 @@ -534,6 +589,8 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): # Optional soft ramp (half-cosine) across a band of width 'ramp_dec' decades if ramp_dec > 0.0: f1 = f0 / (10.0 ** ramp_dec) + if freq_range is not None: + f1 = max(f1, freq_range[0]) band = np.where((freq >= f1) & (freq < f0))[0] if band.size: # w goes 0->1 across [f1, f0]; blend clamped (A0) with original @@ -545,6 +602,8 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): 'applied': True, 'f_cross': float(f0), }) + if range_forced: + mask_info['range_forced'] = True if ramp_dec > 0.0: mask_info['f_ramp_start'] = float(f1) spec.stats.snr_mask_info = mask_info @@ -552,6 +611,10 @@ def _pre_smoothing_snr_mask(config, spec, specnoise): logger.info(f'{spec_id}: pre-smoothing SNR mask applied ' f'(threshold={th:g}, mode={mode}, f_cross={f0:.4f} Hz, ' f'ramp_decades={ramp_dec:.3f})') + if range_forced and freq_range is not None: + logger.info( + '%s: pre-smoothing SNR mask forced by RSS range %.4f-%.4f Hz', + spec_id, freq_range[0], freq_range[1]) else: # 'binary' fallback: clamp any low-SNR bin. mask = sn < th