diff --git a/libs/seiscomp/math/filter/CMakeLists.txt b/libs/seiscomp/math/filter/CMakeLists.txt index 8b738148c..32b3c1024 100644 --- a/libs/seiscomp/math/filter/CMakeLists.txt +++ b/libs/seiscomp/math/filter/CMakeLists.txt @@ -19,6 +19,7 @@ SET(FILTER_SOURCES seismometers.cpp sr.cpp sum.cpp + filterpickercf.cpp ) SET(FILTER_HEADERS @@ -46,6 +47,7 @@ SET(FILTER_HEADERS seismometers.h sr.h sum.h + filterpickercf.h ) SC_SETUP_LIB_SUBDIR(FILTER) diff --git a/libs/seiscomp/math/filter/filterpickercf.cpp b/libs/seiscomp/math/filter/filterpickercf.cpp new file mode 100644 index 000000000..88afcc40f --- /dev/null +++ b/libs/seiscomp/math/filter/filterpickercf.cpp @@ -0,0 +1,271 @@ +/*************************************************************************** + * Copyright (C) Donavin Liebgott * + * All rights reserved. * + * Contact: Donavin Liebgott (donavinliebgott@gmail.com) * + * * + * GNU Affero General Public License Usage * + * This file may be used under the terms of the GNU Affero * + * Public License version 3.0 as published by the Free Software Foundation * + * and appearing in the file LICENSE included in the packaging of this * + * file. Please review the following information to ensure the GNU Affero * + * Public License version 3.0 requirements will be met: * + * https://www.gnu.org/licenses/agpl-3.0.html. * + * * + * Other Usage * + * Alternatively, this file may be used in accordance with the terms and * + * conditions contained in a signed written agreement between you and * + * gempa GmbH. * + ***************************************************************************/ + + +#include +#include +#include +#include + + +namespace Seiscomp { +namespace Math { +namespace Filtering { + + +// ----------------------------------------------------------------------------- +// FilterPickerCF template implementation +// ----------------------------------------------------------------------------- + +template +FilterPickerCF::FilterPickerCF() { + _filters.reserve(_numBands); +} + + +template +FilterPickerCF::~FilterPickerCF() { + for (auto* filter : _filters) { + delete filter; + } + _filters.clear(); +} + + +template +void FilterPickerCF::reset() { + // Reset all bandpass filters + for (auto* filter : _filters) { + filter->reset(); + } + + // Reset envelope buffer + if (!_envelopeBuffer.empty()) { + std::fill(_envelopeBuffer.begin(), _envelopeBuffer.end(), 0); + } + _envelopeBufferPos = 0; +} + + +template +void FilterPickerCF::setSamplingFrequency(double fsamp) { + if (fsamp <= 0) { + return; + } + + _fsamp = fsamp; + + // Initialize envelope buffer + _envelopeBuffer.resize(_envelopeWindowSize, 0); + + // Initialize filter bank + initFilterBank(); +} + + +template +int FilterPickerCF::setParameters(int n, const double *params) { + if (n < 3) { + return -1; // Not enough parameters + } + + int numBands = static_cast(params[0]); + double minFreq = params[1]; + double maxFreq = params[2]; + + if (numBands <= 0 || minFreq <= 0 || maxFreq <= minFreq) { + return -2; // Invalid parameter values + } + + setFilterParams(numBands, minFreq, maxFreq); + + return n; +} + + +template +void FilterPickerCF::setFilterParams(int numBands, double minFreq, double maxFreq) { + _numBands = numBands; + _minFreq = minFreq; + _maxFreq = maxFreq; + + // Reinitialize if sampling frequency is already set + if (_fsamp > 0) { + setSamplingFrequency(_fsamp); + } +} + + +template +void FilterPickerCF::initFilterBank() { + // Clean up old filters + for (auto* filter : _filters) { + delete filter; + } + _filters.clear(); + _filterBandsConfig.clear(); + + if (_numBands <= 0) { + _numBands = 5; + } + + // Compute Nyquist frequency and adjust max frequency + double nyquist = _fsamp / 2.0; + double effectiveMaxFreq = std::min(_maxFreq, nyquist * 0.95); + + if (_minFreq <= 0 || effectiveMaxFreq <= _minFreq) { + SEISCOMP_WARNING("FilterPickerCF: invalid frequency range (fsamp=%.2f Hz), using adjusted defaults", + _fsamp); + _minFreq = 0.5; + effectiveMaxFreq = std::min(20.0, nyquist * 0.95); + } + + // Create logarithmically spaced frequency bands + double logMinFreq = log10(_minFreq); + double logMaxFreq = log10(effectiveMaxFreq); + double logStep = (logMaxFreq - logMinFreq) / _numBands; + + for (int i = 0; i < _numBands; ++i) { + FilterBankConfig config; + + // Center frequency (logarithmic spacing) + config.centerFreq = pow(10.0, logMinFreq + (i + 0.5) * logStep); + + // Band edges (approximately 1/3 octave bandwidth) + double bandwidth = pow(10.0, logStep / 2.0); + config.lowFreq = config.centerFreq / bandwidth; + config.highFreq = config.centerFreq * bandwidth; + + // Clamp to min/max frequencies and Nyquist + config.lowFreq = std::max(config.lowFreq, _minFreq); + config.highFreq = std::min(config.highFreq, effectiveMaxFreq); + + // Ensure lowFreq < highFreq after clamping + if (config.lowFreq >= config.highFreq) { + config.lowFreq = config.highFreq * 0.8; + } + + config.poles = 2; + _filterBandsConfig.push_back(config); + + // Create bandpass filter for this band + auto* filter = new Math::Filtering::IIR::ButterworthBandpass( + config.poles, config.lowFreq, config.highFreq, _fsamp); + _filters.push_back(filter); + } + + SEISCOMP_DEBUG("FilterPickerCF initialized with %d frequency bands (fsamp=%.2f Hz)", + _numBands, _fsamp); +} + + +template +void FilterPickerCF::apply(int n, TYPE *inout) { + if (n <= 0 || inout == nullptr) { + return; + } + + // Check if filter bank is initialized + if (_filters.empty() && _fsamp > 0) { + initFilterBank(); + } + + if (_filters.empty()) { + // Not initialized, pass through + return; + } + + // Validate envelope buffer + if (_envelopeWindowSize <= 0 || _envelopeBuffer.empty()) { + // Not initialized, pass through + return; + } + + // Process each sample incrementally + for (int i = 0; i < n; ++i) { + TYPE sample = inout[i]; + + // Apply filter bank and get maximum envelope across all bands + TYPE maxEnvelope = 0; + + for (size_t band = 0; band < _filters.size(); ++band) { + // Null check before dereferencing + if (_filters[band] == nullptr) { + continue; + } + + // Apply bandpass filter (maintains internal state) + TYPE filtered = sample; + _filters[band]->apply(1, &filtered); + + // Compute envelope (absolute value) + TYPE envelope = std::fabs(filtered); + + // Track maximum across bands (before smoothing) + if (envelope > maxEnvelope) { + maxEnvelope = envelope; + } + } + + // Apply short moving average smoothing to the maximum envelope + // Update circular buffer + _envelopeBuffer[_envelopeBufferPos] = maxEnvelope; + _envelopeBufferPos = (_envelopeBufferPos + 1) % _envelopeWindowSize; + + // Compute moving average + TYPE smoothedEnvelope = 0; + for (int j = 0; j < _envelopeWindowSize; ++j) { + smoothedEnvelope += _envelopeBuffer[j]; + } + + // Prevent division by zero + if (_envelopeWindowSize > 0) { + smoothedEnvelope /= _envelopeWindowSize; + } + + // Output the smoothed maximum envelope (characteristic function) + // The picker (AIC, STA/LTA, etc.) will perform detection on this + inout[i] = smoothedEnvelope; + } +} + + +template +InPlaceFilter* FilterPickerCF::clone() const { + auto* clone = new FilterPickerCF(); + clone->setFilterParams(_numBands, _minFreq, _maxFreq); + if (_fsamp > 0) { + clone->setSamplingFrequency(_fsamp); + } + return clone; +} + + +// Explicit template instantiations +template class SC_SYSTEM_CORE_API FilterPickerCF; +template class SC_SYSTEM_CORE_API FilterPickerCF; + + +// Register the FilterPickerCF InplaceFilter plugin +REGISTER_INPLACE_FILTER(FilterPickerCF, "FP"); + + +} // namespace Seiscomp::Math::Filter +} // namespace Seiscomp::Math +} // namespace Seiscomp diff --git a/libs/seiscomp/math/filter/filterpickercf.h b/libs/seiscomp/math/filter/filterpickercf.h new file mode 100644 index 000000000..787149671 --- /dev/null +++ b/libs/seiscomp/math/filter/filterpickercf.h @@ -0,0 +1,166 @@ +/*************************************************************************** + * Copyright (C) Donavin Liebgott * + * All rights reserved. * + * Contact: Donavin Liebgott (donavinliebgott@gmail.com) * + * * + * GNU Affero General Public License Usage * + * This file may be used under the terms of the GNU Affero * + * Public License version 3.0 as published by the Free Software Foundation * + * and appearing in the file LICENSE included in the packaging of this * + * file. Please review the following information to ensure the GNU Affero * + * Public License version 3.0 requirements will be met: * + * https://www.gnu.org/licenses/agpl-3.0.html. * + * * + * Other Usage * + * Alternatively, this file may be used in accordance with the terms and * + * conditions contained in a signed written agreement between you and * + * gempa GmbH. * + ***************************************************************************/ + + +#ifndef SEISCOMP_MATH_FILTER_FILTERPICKERCF_H +#define SEISCOMP_MATH_FILTER_FILTERPICKERCF_H + + +#include +#include +#include + + +namespace Seiscomp { +namespace Math { +namespace Filtering { + + +/** + * @brief FilterPicker Characteristic Function InplaceFilter + * + * Based on the FilterPicker algorithm by Lomax et al. (2012) + * Reference: Lomax, A., Satriano, C., & Vassallo, M. (2012). + * Automatic picker developments and optimization: FilterPicker - a robust, + * broadband picker for real-time seismic monitoring and earthquake early-warning. + * Seismological Research Letters, 83(3), 531-540. + * + * This filter computes the characteristic function of the FilterPicker algorithm. + * It operates on multiple frequency bands and outputs the maximum envelope across + * all bands. The output can be used by pickers (e.g., AIC, STA/LTA) for phase detection. + * + * The filter works incrementally, maintaining state between apply() calls. + * This is crucial for proper integration with SeisComP's filter chains and + * continuous data processing. + * + * INCREMENTAL PROCESSING: + * The filter produces identical results whether processing samples in bulk or + * one sample at a time. This is achieved by: + * - IIR filters (Butterworth) maintain internal state (v1, v2 memory) + * - Envelope smoothing buffer maintains state across apply() calls + * - All operations are sample-by-sample with no look-ahead + * + * Usage: + * - As a pre-filter for scautopick: picker.filter = "FP(numBands,minFreq,maxFreq)" + * - In filter chains: "BP(0.5,10)+FP(5,0.5,20)" + * - With AIC picker: picker = AIC, picker.filter = "FP(5,0.5,20)" + * - With STA/LTA picker: picker = STALTA, picker.filter = "FP(5,0.5,20)" + * + * Parameters: + * - numBands: Number of frequency bands (default: 5) + * - minFreq: Minimum frequency in Hz (default: 0.5) + * - maxFreq: Maximum frequency in Hz (default: 20.0) + * + * The filter outputs the maximum envelope across all frequency bands. + * The picker (AIC, STA/LTA, etc.) performs the actual detection on this output. + */ +template +class FilterPickerCF : public InPlaceFilter { + // ---------------------------------------------------------------------- + // Public types + // ---------------------------------------------------------------------- + public: + /** + * @brief Structure to hold filter bank configuration + */ + struct FilterBankConfig { + double centerFreq; //!< Center frequency in Hz + double lowFreq; //!< Low cutoff frequency in Hz + double highFreq; //!< High cutoff frequency in Hz + int poles; //!< Number of filter poles + }; + + // ---------------------------------------------------------------------- + // X'truction + // ---------------------------------------------------------------------- + public: + //! C'tor + FilterPickerCF(); + + //! D'tor + ~FilterPickerCF() override; + + + // ---------------------------------------------------------------------- + // InplaceFilter interface + // ---------------------------------------------------------------------- + public: + void setSamplingFrequency(double fsamp) override; + int setParameters(int n, const double *params) override; + void apply(int n, TYPE *inout) override; + InPlaceFilter* clone() const override; + + + // ---------------------------------------------------------------------- + // Public methods + // ---------------------------------------------------------------------- + public: + //! Reset the filter state + void reset(); + + //! Set filter parameters + void setFilterParams(int numBands, double minFreq, double maxFreq); + + //! Get the number of frequency bands + int numBands() const { return _numBands; } + + //! Get filter bank configuration + const std::vector& filterBanks() const { return _filterBandsConfig; } + + + // ---------------------------------------------------------------------- + // Private methods + // ---------------------------------------------------------------------- + private: + /** + * @brief Initialize the filter bank + */ + void initFilterBank(); + + + // ---------------------------------------------------------------------- + // Private members + // ---------------------------------------------------------------------- + private: + std::vector _filterBandsConfig; //!< Filter bank configuration + + // Algorithm parameters + int _numBands{5}; //!< Number of frequency bands + double _minFreq{0.5}; //!< Minimum frequency in Hz + double _maxFreq{20.0}; //!< Maximum frequency in Hz + + // Sampling rate + double _fsamp{0.0}; //!< Sampling frequency in Hz + + // Filter state for each band + std::vector*> _filters; + + // Envelope smoothing buffer + std::vector _envelopeBuffer; //!< Recent envelope values + size_t _envelopeBufferPos{0}; //!< Position in envelope buffer + int _envelopeWindowSize{10}; //!< Envelope smoothing window +}; + + +} // namespace Seiscomp::Math::Filter +} // namespace Seiscomp::Math +} // namespace Seiscomp + + +#endif // SEISCOMP_MATH_FILTER_FILTERPICKERCF_H diff --git a/libs/seiscomp/processing/CMakeLists.txt b/libs/seiscomp/processing/CMakeLists.txt index d4d771e7b..057f5129e 100644 --- a/libs/seiscomp/processing/CMakeLists.txt +++ b/libs/seiscomp/processing/CMakeLists.txt @@ -31,12 +31,14 @@ SET(PROC_HEADERS regions.h detector.h + picker.h secondarypicker.h amplitudeprocessor.h magnitudeprocessor.h ) +SC_ADD_SUBDIR_SOURCES(PROC filterpicker) SC_ADD_SUBDIR_SOURCES(PROC fx) SC_ADD_SUBDIR_SOURCES(PROC picker) SC_ADD_SUBDIR_SOURCES(PROC secondarypicker) diff --git a/libs/seiscomp/processing/filterpicker/CMakeLists.txt b/libs/seiscomp/processing/filterpicker/CMakeLists.txt new file mode 100644 index 000000000..4868bfa50 --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/CMakeLists.txt @@ -0,0 +1,9 @@ +# CMakeLists.txt for FilterPicker plugin +# Based on the Lomax et al. (2012) FilterPicker algorithm + +SET(FILTERPICKER_SOURCES filterpicker.cpp) +SET(FILTERPICKER_HEADERS filterpicker.h) + +FILE(GLOB descs "${CMAKE_CURRENT_SOURCE_DIR}/descriptions/*.xml") INSTALL(FILES ${descs} DESTINATION ${SC3_PACKAGE_APP_DESC_DIR}) + +SC_SETUP_LIB_SUBDIR(FILTERPICKER) diff --git a/libs/seiscomp/processing/filterpicker/README.md b/libs/seiscomp/processing/filterpicker/README.md new file mode 100644 index 000000000..9eee01b30 --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/README.md @@ -0,0 +1,156 @@ +# FilterPicker for SeisComP + +A C++ implementation of the FilterPicker algorithm for SeisComP, based on the work of Lomax et al. (2012). Adapted for SeisComP by Donavin Liebgott. + +## Overview + +FilterPicker is a general-purpose, broadband phase detector and picker algorithm applicable to real-time seismic monitoring and earthquake early-warning systems. It operates on multiple frequency bands and uses characteristic functions to detect and pick seismic phases. + +## Features + +- **Broadband detection**: Operates on multiple frequency bands simultaneously +- **Robust picking**: Uses characteristic function integration for reliable onset detection +- **Adaptive thresholding**: Optional noise-adaptive threshold scaling +- **Uncertainty estimation**: Provides realistic timing uncertainty estimates +- **Polarity detection**: Determines onset polarity (positive/negative) + +## Algorithm + +The FilterPicker algorithm works as follows: + +1. **Filter Bank**: Input data is filtered through a bank of bandpass filters with logarithmically spaced center frequencies +2. **Characteristic Function**: For each frequency band, a characteristic function (CF) is computed based on the envelope STA/LTA ratio +3. **Integration**: The maximum CF across all bands is integrated over a time window +4. **Detection**: A pick is declared when the integral exceeds a threshold +5. **Refinement**: The exact onset time and uncertainty are estimated from the CF shape + +## Installation + +### Prerequisites + +- SeisComP source code (scsrc for reference) +- CMake >= 3.10 +- C++ compiler with C++11 support + +### Build Instructions + +1. The FilterPicker module is located in: + ``` + scsrc/src/base/common/libs/seiscomp/processing/filterpicker/ + ``` + +2. The module is automatically included in the build when you compile SeisComP: + ```bash + cd scsrc/build + cmake .. + make -j4 + make install + ``` + +3. Verify the plugin is available: + ```bash + ls -x ~/seiscomp/share/plugins/filterpicker.so + ``` + +## Configuration + +### Basic Configuration + +To enable FilterPicker in scautopick, add to your `scautopick.cfg` or `global.cfg`: + +```ini +# Enable FilterPicker for P-phase picking +picker = FilterPicker +``` + +### Parameters + +#### FilterPicker (P-phase) + +| Parameter | Default | Description | +|-----------|---------|-------------| +| `picker.FilterPicker.numBands` | 5 | Number of frequency bands | +| `picker.FilterPicker.minFreq` | 0.5 | Minimum frequency (Hz) | +| `picker.FilterPicker.maxFreq` | 20.0 | Maximum frequency (Hz) | +| `picker.FilterPicker.windowSize` | 2 | Integration window (seconds) | +| `picker.FilterPicker.threshold` | 3.0 | Detection threshold | +| `picker.FilterPicker.thresholdOff` | 1.5 | Reset threshold | +| `picker.FilterPicker.adaptiveThreshold` | false | Enable adaptive thresholding | +| `picker.FilterPicker.noiseBegin` | -10.0 | Noise window start (seconds) | +| `picker.FilterPicker.signalBegin` | -5.0 | Signal window start (seconds) | +| `picker.FilterPicker.signalEnd` | 20.0 | Signal window end (seconds) | +| `picker.FilterPicker.minSNR` | 3.0 | Minimum SNR for valid pick | + +### Example Configuration + +See `filterpicker.cfg` for a complete example configuration file. + +## Usage + +### Real-time Processing + +```bash +seiscomp exec scautopick --config-filter="picker=FilterPicker" +``` + +### Offline Processing + +```bash +seiscomp exec scautopick --offline --playback -I data.mseed --ep > picks.xml +``` + +## Tuning Guidelines + +### For Local/Regional Events + +- Increase `maxFreq` to 30-50 Hz +- Use more frequency bands (6-8) +- Lower the threshold (2.0-2.5) + +### For Teleseismic Events + +- Decrease `minFreq` to 0.5-0.7 Hz +- Use fewer frequency bands (3-4) +- Increase windowSize to 3-5 seconds + +### For Noisy Environments + +- Enable adaptive thresholding +- Increase the threshold +- Increase minSNR + +## Performance + +The FilterPicker is designed for real-time operation. Performance depends on: + +- Number of frequency bands (more bands = more CPU) +- Sampling rate (higher rate = more CPU) +- Window size (larger window = more latency) + +Typical performance on modern hardware: +- 100 Hz data, 5 bands: < 10 ms per trace per second +- Memory usage: ~1 MB per active trace + +## References + +1. Lomax, A., Satriano, C., & Vassallo, M. (2012). Automatic picker developments and optimization: FilterPicker - a robust, broadband picker for real-time seismic monitoring and earthquake early-warning. *Seismological Research Letters*, 83(3), 531-540. https://doi.org/10.1785/gssrl.83.3.531 + +2. Vassallo, M., Satriano, C., & Lomax, A. (2012). Automatic picker developments and optimization: A strategy for improving the performances of automatic phase pickers. *Seismological Research Letters*, 83(3), 541-554. + +## Original Implementation + +The original FilterPicker implementation by A. Lomax is available at: +- http://alomax.free.fr/FilterPicker/ +- http://alomax.net/projects/java/net/alomax/timedom/FilterPicker5.java + +## License + +This implementation is part of SeisComP and is licensed under the GNU Affero General Public License version 3.0 (AGPL-3.0). + +## Author + +Implementation based on the FilterPicker algorithm by Anthony Lomax, adapted for SeisComP by Donavin Liebgott. + +## Support + +For issues and questions, please refer to the SeisComP documentation or community forums. diff --git a/libs/seiscomp/processing/filterpicker/descriptions/filterpicker.rst b/libs/seiscomp/processing/filterpicker/descriptions/filterpicker.rst new file mode 100644 index 000000000..4bbdda78a --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/descriptions/filterpicker.rst @@ -0,0 +1,452 @@ +.. _filterpicker: + +=========== +FilterPicker +=========== + +.. contents:: + :local: + :depth: 3 + +Overview +======== + +**FilterPicker** is a robust, broadband phase detector and picker algorithm for real-time seismic monitoring and earthquake early-warning systems. The implementation is based on the work of Lomax et al. (2012) and operates on multiple frequency bands using characteristic functions to detect and pick seismic phases. + +The module provides a primary picker for P-phase detection and picking, designed for real-time operation in the SeisComP :program:`scautopick` module. + + +Algorithm Description +===================== + +The FilterPicker algorithm works through the following stages: + +1. **Filter Bank Decomposition** + + Input data is decomposed through a bank of bandpass filters with logarithmically spaced center frequencies. This broadband approach ensures sensitivity to seismic phases across a wide frequency range. + +2. **Characteristic Function Computation** + + For each frequency band, a characteristic function (CF) is computed based on the envelope STA/LTA (Short-Term Average / Long-Term Average) ratio. The CF emphasizes the onset of seismic phases. + +3. **Integration** + + The maximum CF across all frequency bands is integrated over a configurable time window to enhance the detection signal. + +4. **Detection and Picking** + + A pick is declared when the characteristic function exceeds an adaptive threshold. The exact onset time is determined from the CF shape. + +5. **Uncertainty Estimation** + + Realistic timing uncertainty estimates are provided based on the rise time of the characteristic function. + +6. **Polarity Determination** + + The onset polarity (positive, negative, or undecidable) is determined from the initial motion of the seismic signal. + + +Installation +============ + +The FilterPicker module is part of the SeisComP base processing library and is automatically included when building SeisComP from source. + +Build Requirements +------------------ + +- SeisComP source code (scsrc) +- CMake >= 3.10 +- C++ compiler with C++11 support + +Build Instructions +------------------ + +1. The FilterPicker module is located in:: + + scsrc/src/base/common/libs/seiscomp/processing/filterpicker/ + +2. Build SeisComP as usual:: + + cd scsrc/build + cmake .. + make -j4 + make install + +3. Verify the plugin is available:: + + seiscomp exec scautopick --help + + +Configuration +============= + +Basic Setup +----------- + +To enable FilterPicker in :program:`scautopick`, add the following to your ``scautopick.cfg`` or profile configuration file: + +.. code-block:: ini + + # Enable FilterPicker for P-phase picking + picker = FilterPicker + + +Module Parameters +================= + +.. table:: FilterPicker configuration parameters + + ========================= ========= ================================================ + Parameter Default Description + ========================= ========= ================================================ + ``picker.FilterPicker.numBands`` 5 Number of frequency bands in the filter bank + ``picker.FilterPicker.minFreq`` 0.5 Minimum frequency in Hz + ``picker.FilterPicker.maxFreq`` 20.0 Maximum frequency in Hz + ``picker.FilterPicker.windowSize`` 2 Integration window size in seconds + ``picker.FilterPicker.threshold`` 2.0 Detection threshold (CF ratio) + ``picker.FilterPicker.thresholdOff`` 1.5 Reset threshold for detector + ``picker.FilterPicker.adaptiveThreshold`` ``true`` Enable adaptive thresholding based on noise + ``picker.FilterPicker.noiseBegin`` -10.0 Start of noise window relative to trigger (s) + ``picker.FilterPicker.signalBegin`` -5.0 Start of signal window relative to trigger (s) + ``picker.FilterPicker.signalEnd`` 20.0 End of signal window relative to trigger (s) + ``picker.FilterPicker.minSNR`` 2.0 Minimum SNR required for valid pick + ========================= ========= ================================================ + + +Parameter Descriptions +====================== + +Frequency Bands (``numBands``) +------------------------------ + +The number of frequency bands determines how many parallel filters are applied to the input data. More bands provide better frequency coverage but increase computational cost. + +- **Recommended range**: 4–8 bands +- **Default**: 5 bands +- **Trade-off**: More bands = better detection across frequencies, but higher CPU usage + + +Frequency Range (``minFreq``, ``maxFreq``) +------------------------------------------ + +Defines the frequency range covered by the filter bank. The bands are logarithmically spaced between these limits. + +- **Typical**: 0.5–20 Hz (broadband) +- **Local/regional events**: Consider 1–30 Hz +- **Teleseismic events**: Consider 0.1–5 Hz + + +Detection Threshold (``threshold``) +----------------------------------- + +The characteristic function ratio required to trigger a pick. Lower values increase sensitivity but may produce more false picks. + +- **Lower values** (1.5–2.0): More sensitive, suitable for quiet stations +- **Higher values** (2.5–3.5): More conservative, suitable for noisy environments +- **Default**: 2.0 + + +Adaptive Thresholding (``adaptiveThreshold``) +--------------------------------------------- + +When enabled, the detection threshold is automatically scaled based on the background noise level. This improves performance in varying noise conditions. + +- **Enabled (true)**: Recommended for most installations +- **Disabled (false)**: Use fixed threshold regardless of noise + + +Integration Window (``windowSize``) +----------------------------------- + +Time window for characteristic function integration. Longer windows provide more stable detection but increase latency. + +- **Typical**: 2 seconds (default) +- **Fast response**: Reduce to 1 second +- **Noisy environments**: Increase to 3–5 seconds + + +SNR Threshold (``minSNR``) +-------------------------- + +Minimum signal-to-noise ratio required for a valid pick. Picks with SNR below this threshold are rejected. + +- **Lower values** (1.5–2.0): Accept more picks, including weaker signals +- **Higher values** (3.0–5.0): Only accept high-quality picks +- **Default**: 2.0 + + +Example Configurations +====================== + +Standard Configuration +---------------------- + +A balanced configuration suitable for most regional seismic networks: + +.. code-block:: ini + + picker = FilterPicker + + picker.FilterPicker.numBands = 5 + picker.FilterPicker.minFreq = 0.5 + picker.FilterPicker.maxFreq = 20.0 + picker.FilterPicker.windowSize = 2 + picker.FilterPicker.threshold = 2.0 + picker.FilterPicker.adaptiveThreshold = true + picker.FilterPicker.minSNR = 2.0 + + +Local/Regional Events (Sensitive) +--------------------------------- + +Optimized for detecting local and regional events with higher frequency content: + +.. code-block:: ini + + picker.FilterPicker.numBands = 6 + picker.FilterPicker.minFreq = 1.0 + picker.FilterPicker.maxFreq = 30.0 + picker.FilterPicker.threshold = 1.5 + picker.FilterPicker.minSNR = 1.5 + picker.FilterPicker.adaptiveThreshold = true + + +Teleseismic Events (Conservative) +--------------------------------- + +Optimized for teleseismic events with lower frequency content: + +.. code-block:: ini + + picker.FilterPicker.numBands = 4 + picker.FilterPicker.minFreq = 0.1 + picker.FilterPicker.maxFreq = 5.0 + picker.FilterPicker.windowSize = 3 + picker.FilterPicker.threshold = 2.5 + picker.FilterPicker.minSNR = 3.0 + + +Noisy Environment +----------------- + +Configuration for stations with high background noise: + +.. code-block:: ini + + picker.FilterPicker.numBands = 5 + picker.FilterPicker.threshold = 3.0 + picker.FilterPicker.minSNR = 3.5 + picker.FilterPicker.adaptiveThreshold = true + picker.FilterPicker.windowSize = 3 + + +Usage +===== + +Real-time Processing +-------------------- + +Enable FilterPicker in your SeisComP configuration and start :program:`scautopick`: + +.. code-block:: bash + + seiscomp start scautopick + +Or run with command-line configuration override: + +.. code-block:: bash + + seiscomp exec scautopick --config-filter="picker=FilterPicker" + + +Offline Processing +------------------ + +Process archived data with FilterPicker: + +.. code-block:: bash + + seiscomp exec scautopick --offline --playback -I data.mseed --ep + + +Testing and Debugging +--------------------- + +Test FilterPicker with specific parameters: + +.. code-block:: bash + + seiscomp exec scautopick \ + --config-filter="picker.FilterPicker.numBands=6" \ + --config-filter="picker.FilterPicker.threshold=2.5" \ + --debug + + +Performance Considerations +========================== + +The FilterPicker is designed for real-time operation. Performance depends on several factors: + +Computational Cost +------------------ + +- **Number of frequency bands**: Each additional band requires one bandpass filter operation +- **Sampling rate**: Higher sampling rates increase the number of samples to process +- **Integration window size**: Larger windows require more memory but have minimal CPU impact + +Typical Performance +------------------- + +On modern hardware (as of 2024): + +- 100 Hz data, 5 bands: < 10 ms per trace per second +- Memory usage: ~1 MB per active trace + +Optimization Tips +----------------- + +1. Use the minimum number of frequency bands necessary for your application +2. For high sampling rate data (>100 Hz), consider decimating before picking +3. Enable adaptive thresholding to reduce false picks in varying noise conditions +4. Adjust the frequency range to match your expected signal characteristics + + +Tuning Guidelines +================= + +For Local/Regional Events +------------------------- + +- Increase ``maxFreq`` to 30–50 Hz +- Use more frequency bands (6–8) +- Lower the threshold (1.5–2.0) +- Reduce ``minSNR`` to 1.5–2.0 + +For Teleseismic Events +---------------------- + +- Decrease ``minFreq`` to 0.1–0.2 Hz +- Use fewer frequency bands (3–4) +- Increase ``windowSize`` to 3–5 seconds +- Increase threshold to 2.5–3.5 + +For Noisy Environments +---------------------- + +- Enable adaptive thresholding (``adaptiveThreshold = true``) +- Increase the threshold to 3.0–4.0 +- Increase ``minSNR`` to 3.0–5.0 +- Consider increasing ``windowSize`` for more stable detection + +For Earthquake Early Warning (EEW) +---------------------------------- + +- Minimize ``windowSize`` (1–2 seconds) for fastest response +- Use moderate threshold (2.0–2.5) to balance speed and reliability +- Enable adaptive thresholding for varying noise conditions +- Consider reducing ``numBands`` to 4 for faster processing + + +Output +====== + +Pick Properties +--------------- + +FilterPicker generates picks with the following properties: + +- **Time**: Onset time with sub-sample precision +- **Uncertainty**: Timing uncertainty estimated from CF rise time (typically 0.1–1.0 s) +- **SNR**: Signal-to-noise ratio of the detected phase +- **Polarity**: Initial motion polarity (positive, negative, or undecidable) +- **Method ID**: ``FilterPicker`` +- **Filter ID**: Description of the filter bank configuration used + + +Quality Indicators +------------------ + +The following indicators can be used to assess pick quality: + +- **SNR**: Higher values (>5) indicate confident picks +- **Uncertainty**: Lower values indicate sharper onsets +- **Polarity**: Confident polarity determination suggests good signal quality + + +Troubleshooting +=============== + +No Picks Generated +------------------ + +1. **Check threshold**: Lower the ``threshold`` parameter +2. **Check frequency range**: Ensure ``minFreq`` and ``maxFreq`` match your signal +3. **Check SNR threshold**: Lower ``minSNR`` if picks are being rejected +4. **Verify data quality**: Check for gaps, spikes, or other data issues + +Too Many False Picks +-------------------- + +1. **Increase threshold**: Raise the ``threshold`` parameter +2. **Increase SNR requirement**: Raise ``minSNR`` +3. **Enable adaptive thresholding**: Set ``adaptiveThreshold = true`` +4. **Adjust frequency range**: Narrow the band to exclude noise frequencies + +Picks Have Large Uncertainty +---------------------------- + +1. **Check signal quality**: Noisy signals produce uncertain picks +2. **Reduce integration window**: Smaller ``windowSize`` may help +3. **Adjust frequency bands**: Ensure bands cover the signal frequency content + + +References +========== + +Primary Reference +----------------- + +Lomax, A., Satriano, C., & Vassallo, M. (2012). +Automatic picker developments and optimization: FilterPicker - a robust, +broadband picker for real-time seismic monitoring and earthquake early-warning. +*Seismological Research Letters*, 83(3), 531-540. +https://doi.org/10.1785/gssrl.83.3.531 + +Additional References +--------------------- + +Vassallo, M., Satriano, C., & Lomax, A. (2012). +Automatic picker developments and optimization: A strategy for improving +the performances of automatic phase pickers. +*Seismological Research Letters*, 83(3), 541-554. + +Original Implementation +----------------------- + +The original FilterPicker implementation by A. Lomax is available at: + +- http://alomax.free.fr/FilterPicker/ +- http://alomax.net/projects/java/net/alomax/timedom/FilterPicker5.java + + +See Also +======== + +- :ref:`scautopick` - Automatic picking module +- :ref:`picker-gfz` - GFZ picker implementation +- :ref:`picker-aic` - AIC picker implementation +- :ref:`processing-picks` - Pick processing in SeisComP + + +Authors +======= + +FilterPicker algorithm developed by Anthony Lomax, Claudio Satriano, and Maurizio Vassallo. + +SeisComP implementation by Donavin Liebgott. + + +License +======= + +This implementation is part of SeisComP and is licensed under the GNU Affero General Public License version 3.0 (AGPL-3.0). diff --git a/libs/seiscomp/processing/filterpicker/descriptions/global_filterpicker.xml b/libs/seiscomp/processing/filterpicker/descriptions/global_filterpicker.xml new file mode 100644 index 000000000..8a809a9f7 --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/descriptions/global_filterpicker.xml @@ -0,0 +1,105 @@ + + + + global + + The FilterPicker plugin implements a robust, broadband picker for real-time + seismic monitoring based on the FilterPicker algorithm by Lomax et al. (2012). + It operates on multiple frequency bands using characteristic functions to detect + and pick seismic phases. The algorithm is particularly well suited for + earthquake early warning applications and real-time seismic monitoring. + Reference: Lomax, A., Satriano, C., & Vassallo, M. (2012). Automatic picker + developments and optimization: FilterPicker - a robust, broadband picker for + real-time seismic monitoring and earthquake early-warning. Seismological + Research Letters, 83(3), 531-540. + + + + + + + + FilterPicker P-phase picking parameters. This picker implements + a multi-band approach with characteristic function computation + for robust broadband phase detection. + + + + Number of frequency bands for the filter bank. More bands + provide better frequency coverage but increase computational + cost. Recommended range: 4-8 bands. + + + + + Minimum frequency of the filter bank in Hz. The bands are + logarithmically spaced between minFreq and maxFreq. + Typical: 0.5 Hz for P-phase (broadband). + + + + + Maximum frequency of the filter bank in Hz. The bands are + logarithmically spaced between minFreq and maxFreq. + Automatically limited to 95% of Nyquist. Typical: 20 Hz. + + + + + Integration window size in seconds. The characteristic + function is integrated over this window to enhance detection. + Fast response: 1s, Default: 2s, Noisy: 3-5s. + + + + + Detection threshold for the characteristic function ratio. + Lower values increase sensitivity but may produce more + false picks. Sensitive: 1.5-2.0, Default: 2.0, Conservative: 2.5-3.5. + + + + + Threshold for resetting the detector after a pick. The + detector will reset when the characteristic function falls + below this value. + + + + + Enable adaptive thresholding based on noise level. When + enabled, the detection threshold is automatically scaled based + on the background noise level. Recommended for most installations. + + + + + Start of noise window relative to trigger time in seconds. + This window is used to compute noise statistics for SNR + calculation and pick uncertainty estimation. + + + + + Start of signal window relative to trigger time in seconds. + The pick search begins at this time relative to trigger. + + + + + End of signal window relative to trigger time in seconds. + The pick search ends at this time. + + + + + Minimum signal-to-noise ratio required for a valid pick. + Picks with SNR below this threshold are rejected. + Accept more: 1.5-2.0, Default: 2.0, High quality: 3.0-5.0. + + + + + + + diff --git a/libs/seiscomp/processing/filterpicker/filterpicker.cfg b/libs/seiscomp/processing/filterpicker/filterpicker.cfg new file mode 100644 index 000000000..9d144947a --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/filterpicker.cfg @@ -0,0 +1,111 @@ +# FilterPicker Configuration for SeisComP scautopick +# Based on Lomax et al. (2012): FilterPicker - a robust, broadband picker +# for real-time seismic monitoring and earthquake early-warning. +# +# Reference: +# Lomax, A., Satriano, C., & Vassallo, M. (2012). +# Automatic picker developments and optimization: FilterPicker - a robust, +# broadband picker for real-time seismic monitoring and earthquake early-warning. +# Seismological Research Letters, 83(3), 531-540. +# https://doi.org/10.1785/gssrl.83.3.531 + +# ============================================================================= +# Basic Configuration +# ============================================================================= + +# Enable FilterPicker for P-phase picking +# Set this to 'FilterPicker' to use the FilterPicker algorithm +picker = FilterPicker + +# Detector filter chain (used for initial trigger detection) +# The FilterPicker will refine picks after the detector triggers +filter = "RMHP(10)>>ITAPER(30)>>BW(4,0.7,2)>>STALTA(2,80)" + +# ============================================================================= +# FilterPicker P-phase Parameters +# ============================================================================= + +# Number of frequency bands for the filter bank +# More bands provide better frequency coverage but increase computation time +# Recommended: 4-8 bands for most applications +picker.FilterPicker.numBands = 5 + +# Minimum frequency of the filter bank (Hz) +# Lower frequencies help detect distant events +picker.FilterPicker.minFreq = 0.5 + +# Maximum frequency of the filter bank (Hz) +# Higher frequencies improve local event detection +picker.FilterPicker.maxFreq = 20.0 + +# Integration window size (seconds) +# Time window for characteristic function integration +picker.FilterPicker.windowSize = 2 + +# Detection threshold +# Higher values reduce false picks but may miss weaker events +picker.FilterPicker.threshold = 3.0 + +# Threshold for resetting detector after a pick +# Should be lower than the detection threshold +picker.FilterPicker.thresholdOff = 1.5 + +# Enable adaptive thresholding based on noise level +# When true, threshold is scaled by background noise +picker.FilterPicker.adaptiveThreshold = false + +# ============================================================================= +# Time Window Parameters +# ============================================================================= + +# Start of noise window relative to trigger (seconds) +# Used for SNR calculation and adaptive thresholding +picker.FilterPicker.noiseBegin = -10.0 + +# Start of signal window relative to trigger (seconds) +picker.FilterPicker.signalBegin = -5.0 + +# End of signal window relative to trigger (seconds) +picker.FilterPicker.signalEnd = 20.0 + +# Minimum SNR required for a valid pick +picker.FilterPicker.minSNR = 3.0 + +# ============================================================================= +# scautopick Global Settings +# ============================================================================= + +# Time correction applied to picks (seconds) +timeCorrection = -0.8 + +# Record ringbuffer size (seconds) +ringBufferSize = 300 + +# Lead time to start picking before current time (seconds) +leadTime = 60 + +# Initialization time (blind period after startup) +initTime = 60 + +# Gap interpolation +gapInterpolation = false + +# Thresholds +thresholds.triggerOn = 3 +thresholds.triggerOff = 1.5 +thresholds.maxGapLength = 4.5 +thresholds.amplMaxTimeWindow = 10 +thresholds.deadTime = 30 +thresholds.minAmplOffset = 3 + +# Amplitudes to compute +amplitudes = MLv, mb, mB + +# Use all streams for picking +useAllStreams = false + +# Kill pending S-pickers when new detection found +killPendingSPickers = true + +# Send detections as picks (debug mode) +sendDetections = false diff --git a/libs/seiscomp/processing/filterpicker/filterpicker.cpp b/libs/seiscomp/processing/filterpicker/filterpicker.cpp new file mode 100644 index 000000000..ad45df642 --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/filterpicker.cpp @@ -0,0 +1,869 @@ +/*************************************************************************** + * Copyright (C) Donavin Liebgott * + * All rights reserved. * + * Contact: Donavin Liebgott (donavinliebgott@gmail.com) * + * * + * GNU Affero General Public License Usage * + * This file may be used under the terms of the GNU Affero * + * Public License version 3.0 as published by the Free Software Foundation * + * and appearing in the file LICENSE included in the packaging of this * + * file. Please review the following information to ensure the GNU Affero * + * Public License version 3.0 requirements will be met: * + * https://www.gnu.org/licenses/agpl-3.0.html. * + * * + * Other Usage * + * Alternatively, this file may be used in accordance with the terms and * + * conditions contained in a signed written agreement between you and * + * gempa GmbH. * + ***************************************************************************/ + + +#define SEISCOMP_COMPONENT FilterPicker + +#include +#include +#include "filterpicker.h" +#include +#include +#include + + +using namespace std; + + +namespace Seiscomp { +namespace Processing { + + +// Register the FilterPicker plugin +REGISTER_POSTPICKPROCESSOR(FilterPicker, "FilterPicker"); + + +namespace { + + +/** + * @brief Compute the envelope of a signal using the Hilbert transform approximation + * + * @param data Input signal + * @return Envelope of the signal + */ +vector computeEnvelope(const vector& data) { + int n = static_cast(data.size()); + vector envelope(n); + + // Simple envelope computation using absolute value + // For a more accurate envelope, use the Hilbert transform + for (int i = 0; i < n; ++i) { + envelope[i] = fabs(data[i]); + } + + // Apply a short moving average smoother (fixed window, not % of trace) + // Use 5-10 samples to preserve onset characteristics + int windowSize = min(10, max(3, n / 100)); // 1% of trace, min 3, max 10 + if (windowSize % 2 == 0) windowSize++; + + int halfWindow = windowSize / 2; + for (int i = 0; i < n; ++i) { + double sum = 0; + int count = 0; + for (int j = max(0, i - halfWindow); j <= min(n - 1, i + halfWindow); ++j) { + sum += envelope[j]; + count++; + } + // Prevent division by zero + if (count > 0) { + envelope[i] = sum / count; + } + } + + return envelope; +} + + +/** + * @brief Compute the characteristic function for a filtered trace + * + * The characteristic function emphasizes the onset of seismic phases + * by computing a ratio of short-term to long-term averages of the envelope. + * + * @param envelope Envelope of the filtered trace + * @param staWindow Short-term average window in samples + * @param ltaWindow Long-term average window in samples + * @return Characteristic function values + */ +vector computeCF(const vector& envelope, int staWindow, int ltaWindow) { + int n = static_cast(envelope.size()); + vector cf(n, 0.0); + + // Compute running sums for STA and LTA + vector sta(n, 0.0); + vector lta(n, 0.0); + + for (int i = 0; i < n; ++i) { + // Short-term average + double staSum = 0; + int staCount = 0; + for (int j = max(0, i - staWindow + 1); j <= i; ++j) { + staSum += envelope[j]; + staCount++; + } + sta[i] = staSum / staCount; + + // Long-term average + double ltaSum = 0; + int ltaCount = 0; + for (int j = max(0, i - ltaWindow + 1); j <= i; ++j) { + ltaSum += envelope[j]; + ltaCount++; + } + lta[i] = ltaSum / ltaCount; + } + + // Compute CF as STA/LTA ratio with regularization + double cfSum = 0.0; + double cfMax = 0.0; + for (int i = 0; i < n; ++i) { + if (lta[i] > 1e-10) { + cf[i] = sta[i] / lta[i]; + } else { + cf[i] = 1.0; // Default to 1 when LTA is zero + } + cfSum += cf[i]; + cfMax = max(cfMax, cf[i]); + } + + return cf; +} + + +} // anonymous namespace + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +FilterPicker::FilterPicker() +: _numBands(5), _minFreq(0.5), _maxFreq(20.0), + _windowSize(2), _threshold(2.0), _thresholdOff(1.5), + _useAdaptiveThreshold(false), + _lastFsamp(0.0) { + // Filter bank will be initialized in calculatePick() when we have valid fsamp +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +FilterPicker::FilterPicker(const Core::Time& trigger) +: Picker(trigger), + _numBands(5), _minFreq(0.5), _maxFreq(20.0), + _windowSize(2), _threshold(2.0), _thresholdOff(1.5), + _useAdaptiveThreshold(false), + _lastFsamp(0.0) { + // Filter bank will be initialized in calculatePick() when we have valid fsamp +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +FilterPicker::~FilterPicker() {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void FilterPicker::reset() { + Picker::reset(); + _filterBandsConfig.clear(); + _lastFsamp = 0.0; + // Filter bank will be re-initialized in calculatePick() when we have valid fsamp +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool FilterPicker::setup(const Settings &settings) { + if (!Picker::setup(settings)) { + return false; + } + + // Read FilterPicker-specific configuration parameters + + // Number of frequency bands + try { + _numBands = settings.getInt("picker.FilterPicker.numBands"); + SEISCOMP_DEBUG("FilterPicker: numBands = %d", _numBands); + } catch (...) { + _numBands = 5; // Default: 5 bands + SEISCOMP_DEBUG("FilterPicker: numBands default = %d", _numBands); + } + + // Minimum frequency + try { + _minFreq = settings.getDouble("picker.FilterPicker.minFreq"); + SEISCOMP_DEBUG("FilterPicker: minFreq = %.2f Hz", _minFreq); + } catch (...) { + _minFreq = 0.5; // Default: 0.5 Hz + SEISCOMP_DEBUG("FilterPicker: minFreq default = %.2f Hz", _minFreq); + } + + // Maximum frequency + try { + _maxFreq = settings.getDouble("picker.FilterPicker.maxFreq"); + SEISCOMP_DEBUG("FilterPicker: maxFreq = %.2f Hz", _maxFreq); + } catch (...) { + _maxFreq = 20.0; // Default: 20 Hz + SEISCOMP_DEBUG("FilterPicker: maxFreq default = %.2f Hz", _maxFreq); + } + + // Integration window size (in seconds) + try { + _windowSize = settings.getInt("picker.FilterPicker.windowSize"); + SEISCOMP_DEBUG("FilterPicker: windowSize = %d s", _windowSize); + } catch (...) { + _windowSize = 2; // Default: 2 seconds + SEISCOMP_DEBUG("FilterPicker: windowSize default = %d s", _windowSize); + } + + // Detection threshold + try { + _threshold = settings.getDouble("picker.FilterPicker.threshold"); + SEISCOMP_DEBUG("FilterPicker: threshold = %.2f", _threshold); + } catch (...) { + _threshold = 2.0; // Default: 2.0 (CF ratio) + SEISCOMP_DEBUG("FilterPicker: threshold default = %.2f", _threshold); + } + + // Threshold off (for reset) + try { + _thresholdOff = settings.getDouble("picker.FilterPicker.thresholdOff"); + SEISCOMP_DEBUG("FilterPicker: thresholdOff = %.2f", _thresholdOff); + } catch (...) { + _thresholdOff = 1.5; // Default: 1.5 + SEISCOMP_DEBUG("FilterPicker: thresholdOff default = %.2f", _thresholdOff); + } + + // Adaptive threshold + try { + _useAdaptiveThreshold = settings.getBool("picker.FilterPicker.adaptiveThreshold"); + SEISCOMP_DEBUG("FilterPicker: adaptiveThreshold = %s", + _useAdaptiveThreshold ? "true" : "false"); + } catch (...) { + _useAdaptiveThreshold = true; // Default: enabled + SEISCOMP_DEBUG("FilterPicker: adaptiveThreshold default = true"); + } + + // Configuration from parent class (noise/signal windows, SNR) + try { + _config.noiseBegin = settings.getDouble("picker.FilterPicker.noiseBegin"); + } catch (...) { + _config.noiseBegin = -10.0; + } + + try { + _config.signalBegin = settings.getDouble("picker.FilterPicker.signalBegin"); + } catch (...) { + _config.signalBegin = -5.0; + } + + try { + _config.signalEnd = settings.getDouble("picker.FilterPicker.signalEnd"); + } catch (...) { + _config.signalEnd = 20.0; + } + + try { + _config.snrMin = settings.getDouble("picker.FilterPicker.minSNR"); + } catch (...) { + _config.snrMin = 2.0; + } + + // Reinitialize filter bank with new parameters + // NOTE: fsamp may not be set yet (set when first record is processed) + // The filter bank will be initialized in calculatePick() when fsamp is available + if (_stream.fsamp > 0) { + initFilterBank(); + } + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +const string &FilterPicker::methodID() const { + static string method = "FilterPicker"; + return method; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +const std::string &FilterPicker::filterID() const { + return _filterString; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void FilterPicker::initFilterBank() { + _filterBandsConfig.clear(); + _filterString.clear(); + + if (_numBands <= 0) { + SEISCOMP_WARNING("FilterPicker: numBands must be positive, using default"); + _numBands = 5; + } + + // Compute Nyquist frequency and adjust max frequency if needed + double nyquist = _stream.fsamp / 2.0; + double effectiveMaxFreq = min(_maxFreq, nyquist * 0.95); // Use 95% of Nyquist for safety margin + + if (_minFreq <= 0 || effectiveMaxFreq <= _minFreq) { + SEISCOMP_WARNING("FilterPicker: invalid frequency range (fsamp=%.2f Hz), using adjusted defaults", + _stream.fsamp); + _minFreq = 0.5; + effectiveMaxFreq = min(20.0, nyquist * 0.95); + } + + // Create logarithmically spaced frequency bands + // This follows the FilterPicker design for broadband coverage + double logMinFreq = log10(_minFreq); + double logMaxFreq = log10(effectiveMaxFreq); + double logStep = (logMaxFreq - logMinFreq) / _numBands; + + ostringstream filterStr; + filterStr << "FilterPicker[" << _numBands << " bands: " << _minFreq << "-" << effectiveMaxFreq << " Hz]"; + + for (int i = 0; i < _numBands; ++i) { + FilterBankConfig config; + + // Center frequency (logarithmic spacing) + config.centerFreq = pow(10.0, logMinFreq + (i + 0.5) * logStep); + + // Band edges (approximately 1/3 octave bandwidth) + double bandwidth = pow(10.0, logStep / 2.0); + config.lowFreq = config.centerFreq / bandwidth; + config.highFreq = config.centerFreq * bandwidth; + + // Clamp to min/max frequencies and Nyquist + config.lowFreq = max(config.lowFreq, _minFreq); + config.highFreq = min(config.highFreq, effectiveMaxFreq); + + // Ensure lowFreq < highFreq after clamping + if (config.lowFreq >= config.highFreq) { + config.lowFreq = config.highFreq * 0.8; + } + + // Number of poles (typically 2-4 for Butterworth) + config.poles = 2; + + _filterBandsConfig.push_back(config); + + SEISCOMP_DEBUG("FilterPicker band %d: %.2f-%.2f Hz (center: %.2f Hz)", + i, config.lowFreq, config.highFreq, config.centerFreq); + } + + _filterString = filterStr.str(); + SEISCOMP_INFO("FilterPicker initialized with %d frequency bands (fsamp=%.2f Hz, Nyquist=%.2f Hz)", + _numBands, _stream.fsamp, nyquist); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +vector> FilterPicker::applyFilterBank(const double *data, int n) { + vector> filteredTraces(_numBands); + + for (int i = 0; i < _numBands; ++i) { + const FilterBankConfig& config = _filterBandsConfig[i]; + + // Copy input data + vector trace(data, data + n); + + // Apply Butterworth bandpass filter + try { + Math::Filtering::IIR::ButterworthBandpass filter( + config.poles, config.lowFreq, config.highFreq, _stream.fsamp); + filter.apply(n, trace.data()); + } catch (exception& e) { + SEISCOMP_WARNING("FilterPicker: Filter band %d failed: %s", i, e.what()); + // Fill with zeros on failure + trace.assign(n, 0.0); + } + + filteredTraces[i] = trace; + } + + return filteredTraces; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +FilterPicker::CharacteristicFunction FilterPicker::computeCharacteristicFunction( + const vector& filtered) { + + CharacteristicFunction cf; + int n = static_cast(filtered.size()); + + // Compute envelope + vector envelope = computeEnvelope(filtered); + cf.filtered = filtered; + + // Compute characteristic function (STA/LTA of envelope) + // Use standard STA/LTA window sizes for seismic phase detection + // STA: 0.5-1.0 s (short-term, reacts to onset) + // LTA: 10-20 s (long-term, represents background) + int staWindow = max(3, static_cast(_stream.fsamp * 0.5)); // 0.5 s STA + int ltaWindow = max(50, static_cast(_stream.fsamp * 10.0)); // 10 s LTA + + // Ensure LTA is at least 5x STA for proper ratio + ltaWindow = max(ltaWindow, staWindow * 5); + + cf.values = computeCF(envelope, staWindow, ltaWindow); + + // Compute integral and find maximum + cf.integral = 0.0; + cf.maxVal = 0.0; + cf.maxIdx = -1; + + for (int i = 0; i < n; ++i) { + cf.integral += cf.values[i]; + if (cf.values[i] > cf.maxVal) { + cf.maxVal = cf.values[i]; + cf.maxIdx = i; + } + } + + return cf; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +double FilterPicker::computeCFIntegral( + const vector& cfs, int startIdx, int windowSize) { + + if (cfs.empty() || startIdx < 0) return 0.0; + + int n = static_cast(cfs[0].values.size()); + if (n == 0 || startIdx >= n) return 0.0; + + // Ensure windowSize is positive and within bounds + int safeWindowSize = max(1, min(windowSize, n - startIdx)); + int endIdx = min(startIdx + safeWindowSize, n); + double maxIntegral = 0.0; + + // Compute integral of maximum CF across all bands + for (int i = startIdx; i < endIdx; ++i) { + double maxCF = 0.0; + for (const auto& cf : cfs) { + // Bounds check for each CF vector (they may have different sizes) + if (!cf.values.empty() && i < static_cast(cf.values.size())) { + maxCF = max(maxCF, cf.values[i]); + } + } + maxIntegral += maxCF; + } + + return maxIntegral; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool FilterPicker::findOnset(const vector& cfs, + double threshold, int windowSize, + int &onsetIdx, double &snr) { + if (cfs.empty()) return false; + + int n = static_cast(cfs[0].values.size()); + if (n < 10) return false; + + // windowSize is already in samples (passed from calculatePick) + // Ensure windowSamples is within valid bounds to prevent out-of-bounds access + int windowSamples = max(1, min(windowSize, n / 2)); + + // Always use adaptive thresholding based on noise level + // Use first third of trace as noise reference + int noiseStart = 0; + int noiseEnd = max(10, min(n / 3, static_cast(n / 10.0))); // First 10% or at least 10 samples + + double noiseSum = 0.0; + int noiseCount = 0; + double noiseMax = 0.0; + + for (int i = noiseStart; i < noiseEnd; ++i) { + double maxCF = 0.0; + for (const auto& cf : cfs) { + if (i < static_cast(cf.values.size())) { + maxCF = max(maxCF, cf.values[i]); + } + } + noiseSum += maxCF; + noiseMax = max(noiseMax, maxCF); + noiseCount++; + } + + double noiseMean = noiseCount > 0 ? noiseSum / noiseCount : 1.0; + double noiseStd = 0.0; + + // Compute noise standard deviation + for (int i = noiseStart; i < noiseEnd; ++i) { + double maxCF = 0.0; + for (const auto& cf : cfs) { + if (i < static_cast(cf.values.size())) { + maxCF = max(maxCF, cf.values[i]); + } + } + double diff = maxCF - noiseMean; + noiseStd += diff * diff; + } + noiseStd = noiseCount > 1 ? sqrt(noiseStd / (noiseCount - 1)) : noiseMean; + + // Use adaptive threshold: based on noise statistics + // The threshold is compared against the MAX CF value in the window, not the integral + // So we use point-wise threshold, not integrated + double adaptiveThreshold = max(threshold * noiseMean, noiseMean + threshold * noiseStd); + + // Ensure minimum threshold (CF ratio typically 2-10 for good onsets) + adaptiveThreshold = max(adaptiveThreshold, 1.5); + + SEISCOMP_DEBUG("FilterPicker: noise stats - mean=%.3f, std=%.3f, max=%.3f, threshold=%.3f", + noiseMean, noiseStd, noiseMax, adaptiveThreshold); + SEISCOMP_DEBUG("FilterPicker: n=%d, windowSamples=%d, loop end=%d, cfs.size=%zu", + n, windowSamples, n - windowSamples, cfs.size()); + + // Scan for onset: find first point where MAX CF (not integral) exceeds threshold + // Ensure loop bound is valid (n > windowSamples) + int maxCFIdx = -1; + double maxCFVal = 0.0; + int scanCount = 0; + int loopEnd = n - windowSamples; + + // Prevent negative or zero loop end + if (loopEnd <= 0) { + SEISCOMP_DEBUG("FilterPicker: invalid loop end (%d), n=%d, windowSamples=%d", + loopEnd, n, windowSamples); + return false; + } + + for (int i = 0; i < loopEnd; ++i) { + // Get max CF value across all bands at this sample + double maxCF = 0.0; + for (const auto& cf : cfs) { + if (i < static_cast(cf.values.size())) { + maxCF = max(maxCF, cf.values[i]); + } + } + scanCount++; + + // Debug first few iterations + if (i < 3) { + SEISCOMP_DEBUG("findOnset: i=%d, maxCF=%.3f", i, maxCF); + } + + if (maxCF >= adaptiveThreshold) { + onsetIdx = i; + + // Compute SNR as ratio of max signal to noise level + double signalMax = 0.0; + + // Signal window (after onset) + int signalEnd = min(i + windowSamples * 2, n); + for (int j = i; j < signalEnd; ++j) { + for (const auto& cf : cfs) { + if (j < static_cast(cf.values.size())) { + signalMax = max(signalMax, cf.values[j]); + } + } + } + + snr = noiseMax > 1e-10 ? signalMax / noiseMax : 0.0; + + SEISCOMP_DEBUG("FilterPicker: onset found at index %d, SNR = %.2f", onsetIdx, snr); + return true; + } + + // Track maximum CF value for fallback + if (maxCF > maxCFVal) { + maxCFVal = maxCF; + maxCFIdx = i; + } + } + + // Fallback: if no threshold crossing found, use the maximum CF point + // if it's significantly above noise level + if (maxCFIdx >= 0 && maxCFVal > adaptiveThreshold * 0.5) { + onsetIdx = maxCFIdx; + + // Compute SNR at this point + double signalMax = 0.0; + int signalEnd = min(maxCFIdx + windowSamples * 2, n); + for (int j = maxCFIdx; j < signalEnd; ++j) { + for (const auto& cf : cfs) { + if (j < static_cast(cf.values.size())) { + signalMax = max(signalMax, cf.values[j]); + } + } + } + snr = noiseMax > 1e-10 ? signalMax / noiseMax : 0.0; + + SEISCOMP_DEBUG("FilterPicker: onset found (fallback) at index %d, SNR = %.2f", onsetIdx, snr); + return true; + } + + SEISCOMP_DEBUG("FilterPicker: no onset found (maxCF=%.3f, threshold=%.3f)", + maxCFVal, adaptiveThreshold); + return false; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void FilterPicker::estimateUncertainty(const vector& cfs, + int onsetIdx, int &lowerUnc, int &upperUnc) { + if (cfs.empty() || onsetIdx < 0) { + lowerUnc = -1; + upperUnc = -1; + return; + } + + int n = static_cast(cfs[0].values.size()); + if (n == 0 || onsetIdx >= n) { + lowerUnc = -1; + upperUnc = -1; + return; + } + + // Find the rise time of the characteristic function + // Uncertainty is related to how sharp the onset is + + // Get the maximum CF value at onset + double maxCF = 0.0; + for (const auto& cf : cfs) { + if (onsetIdx < static_cast(cf.values.size())) { + maxCF = max(maxCF, cf.values[onsetIdx]); + } + } + + // Handle zero maxCF case + if (maxCF <= 0) { + lowerUnc = 1; + upperUnc = 1; + return; + } + + // Find points where CF reaches 10% and 90% of maximum + int idx10 = onsetIdx; + int idx90 = onsetIdx; + + double threshold10 = 0.1 * maxCF; + double threshold90 = 0.9 * maxCF; + + // Search backward for 10% point + for (int i = onsetIdx; i >= 0; --i) { + double cf = 0.0; + for (const auto& c : cfs) { + if (i < static_cast(c.values.size())) { + cf = max(cf, c.values[i]); + } + } + if (cf <= threshold10) { + idx10 = i; + break; + } + } + + // Search forward for 90% point + for (int i = onsetIdx; i < n; ++i) { + double cf = 0.0; + for (const auto& c : cfs) { + if (i < static_cast(c.values.size())) { + cf = max(cf, c.values[i]); + } + } + if (cf >= threshold90) { + idx90 = i; + break; + } + } + + // Uncertainty is proportional to rise time + int riseTime = idx90 - idx10; + lowerUnc = max(1, riseTime / 2); + upperUnc = max(1, riseTime / 2); + + SEISCOMP_DEBUG("FilterPicker: uncertainty estimated: lower=%d, upper=%d samples", + lowerUnc, upperUnc); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +Picker::Polarity FilterPicker::determinePolarity(const double *data, int onsetIdx, double fsamp) { + // Validate input parameters + if (data == nullptr || onsetIdx < 0 || fsamp <= 0) { + return UNDECIDABLE; + } + + // Calculate safe data size limit + int dataSizeLimit = static_cast(_stream.fsamp * _config.signalEnd); + if (dataSizeLimit <= 0) { + return UNDECIDABLE; + } + + if (onsetIdx >= dataSizeLimit) { + return UNDECIDABLE; + } + + // Look at the first few samples after onset to determine polarity + int windowSize = max(1, static_cast(fsamp * 0.1)); // 100 ms window + double sum = 0.0; + int count = 0; + + int endIndex = min(onsetIdx + windowSize, dataSizeLimit); + for (int i = onsetIdx; i < endIndex; ++i) { + sum += data[i]; + count++; + } + + // Prevent division by zero + if (count == 0) { + return UNDECIDABLE; + } + + double mean = sum / count; + + if (mean > 1e-10) { + return POSITIVE; + } else if (mean < -1e-10) { + return NEGATIVE; + } + + return UNDECIDABLE; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool FilterPicker::calculatePick(int n, const double *data, + int signalStartIdx, int signalEndIdx, + int &triggerIdx, int &lowerUncertainty, + int &upperUncertainty, double &snr, + OPT(Polarity) &polarity) { + + SEISCOMP_DEBUG("FilterPicker::calculatePick: n=%d, signalStart=%d, signalEnd=%d, fsamp=%.2f", + n, signalStartIdx, signalEndIdx, _stream.fsamp); + + // Validate input parameters + if (n <= 10) { + SEISCOMP_INFO("FilterPicker: not enough data (n=%d)", n); + return false; + } + + // Validate signal indices + if (signalStartIdx < 0 || signalEndIdx <= signalStartIdx || signalEndIdx > n) { + SEISCOMP_WARNING("FilterPicker: invalid signal indices (start=%d, end=%d, n=%d)", + signalStartIdx, signalEndIdx, n); + return false; + } + + // Check if we have a valid sampling rate + if (_stream.fsamp <= 0) { + SEISCOMP_WARNING("FilterPicker: invalid sampling rate (%.2f Hz), cannot pick", _stream.fsamp); + return false; + } + + // Initialize filter bank if not yet done or if sampling rate changed + if (_filterBandsConfig.empty() || _stream.fsamp != _lastFsamp) { + _lastFsamp = _stream.fsamp; + initFilterBank(); + } + + // Demean the data using the noise window + int nnoise = max(1, static_cast(n / 3)); + // Ensure nnoise doesn't exceed available data + nnoise = min(nnoise, signalStartIdx); + + double offset = 0.0; + for (int i = 0; i < nnoise; ++i) { + offset += data[i]; + } + offset /= nnoise; + + // Create demeaned copy with bounds-checked size + int demeanedSize = signalEndIdx - signalStartIdx; + if (demeanedSize <= 0) { + SEISCOMP_WARNING("FilterPicker: invalid demeaned size (%d)", demeanedSize); + return false; + } + + vector demeaned(demeanedSize); + for (int i = 0; i < demeanedSize; ++i) { + // Bounds check before access + if (signalStartIdx + i < n) { + demeaned[i] = data[signalStartIdx + i] - offset; + } else { + SEISCOMP_WARNING("FilterPicker: data access out of bounds (idx=%d, n=%d)", + signalStartIdx + i, n); + demeaned[i] = 0.0; + } + } + + // Apply filter bank + vector> filteredTraces = applyFilterBank(demeaned.data(), + demeanedSize); + + // Compute characteristic functions for all bands + vector cfs; + cfs.reserve(_numBands); + + for (const auto& filtered : filteredTraces) { + cfs.push_back(computeCharacteristicFunction(filtered)); + } + + // Find onset + // Ensure windowSamples is within valid bounds + int windowSamples = static_cast(_windowSize * _stream.fsamp); + windowSamples = max(1, min(windowSamples, demeanedSize / 2)); + + int onsetIdx = -1; + snr = -1.0; + + if (!findOnset(cfs, _threshold, windowSamples, onsetIdx, snr)) { + SEISCOMP_INFO("FilterPicker: no onset found"); + return false; + } + + // Check SNR threshold + if (snr < _config.snrMin) { + SEISCOMP_INFO("FilterPicker: SNR %.2f below threshold %.2f", snr, _config.snrMin); + return false; + } + + // Validate onset index before using it + if (onsetIdx < 0 || onsetIdx >= demeanedSize) { + SEISCOMP_WARNING("FilterPicker: invalid onset index (%d, size=%d)", + onsetIdx, demeanedSize); + return false; + } + + // Convert to absolute index + triggerIdx = onsetIdx + signalStartIdx; + + // Estimate uncertainty + estimateUncertainty(cfs, onsetIdx, lowerUncertainty, upperUncertainty); + + // Determine polarity + polarity = determinePolarity(demeaned.data(), onsetIdx, _stream.fsamp); + + SEISCOMP_INFO("FilterPicker: pick found at index %d (relative: %d), SNR=%.2f, " + "uncertainty=[%d,%d] samples", + triggerIdx, onsetIdx, snr, lowerUncertainty, upperUncertainty); + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +} // namespace Processing +} // namespace Seiscomp diff --git a/libs/seiscomp/processing/filterpicker/filterpicker.h b/libs/seiscomp/processing/filterpicker/filterpicker.h new file mode 100644 index 000000000..ddba8859a --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/filterpicker.h @@ -0,0 +1,220 @@ +/*************************************************************************** + * Copyright (C) Donavin Liebgott * + * All rights reserved. * + * Contact: Donavin Liebgott (donavinliebgott@gmail.com) * + * * + * GNU Affero General Public License Usage * + * This file may be used under the terms of the GNU Affero * + * Public License version 3.0 as published by the Free Software Foundation * + * and appearing in the file LICENSE included in the packaging of this * + * file. Please review the following information to ensure the GNU Affero * + * Public License version 3.0 requirements will be met: * + * https://www.gnu.org/licenses/agpl-3.0.html. * + * * + * Other Usage * + * Alternatively, this file may be used in accordance with the terms and * + * conditions contained in a signed written agreement between you and * + * gempa GmbH. * + ***************************************************************************/ + + +#ifndef SEISCOMP_PROCESSING_FILTERPICKER_H +#define SEISCOMP_PROCESSING_FILTERPICKER_H + + +#include +#include +#include +#include +#include +#include + + +namespace Seiscomp { +namespace Processing { + + +/** + * @brief FilterPicker - A robust, broadband picker for real-time seismic monitoring + * + * Based on the FilterPicker algorithm by Lomax et al. (2012) + * Reference: Lomax, A., Satriano, C., & Vassallo, M. (2012). + * Automatic picker developments and optimization: FilterPicker - a robust, + * broadband picker for real-time seismic monitoring and earthquake early-warning. + * Seismological Research Letters, 83(3), 531-540. + * + * The algorithm operates on multiple frequency bands and uses characteristic + * functions to detect and pick seismic phases. + */ +class SC_SYSTEM_CLIENT_API FilterPicker : public Picker { + // ---------------------------------------------------------------------- + // Public types + // ---------------------------------------------------------------------- + public: + /** + * @brief Structure to hold filter bank configuration + */ + struct FilterBankConfig { + double centerFreq; // Center frequency in Hz + double lowFreq; // Low cutoff frequency in Hz + double highFreq; // High cutoff frequency in Hz + int poles; // Number of filter poles + }; + + /** + * @brief Structure to hold characteristic function data for a frequency band + */ + struct CharacteristicFunction { + std::vector values; // CF values + std::vector filtered; // Filtered trace for this band + double integral; // Running integral + double maxVal; // Maximum value + int maxIdx; // Index of maximum + }; + + // ---------------------------------------------------------------------- + // X'truction + // ---------------------------------------------------------------------- + public: + //! C'tor + FilterPicker(); + FilterPicker(const Core::Time& trigger); + + //! D'tor + ~FilterPicker(); + + + // ---------------------------------------------------------------------- + // Public Interface + // ---------------------------------------------------------------------- + public: + bool setup(const Settings &settings) override; + + const std::string &methodID() const override; + const std::string &filterID() const override; + + //! Get the number of frequency bands + int numBands() const { return _numBands; } + + //! Get filter bank configuration + const std::vector& filterBanks() const { return _filterBandsConfig; } + + + // ---------------------------------------------------------------------- + // Protected Interface + // ---------------------------------------------------------------------- + protected: + //! Calculate pick using FilterPicker algorithm + bool calculatePick(int n, const double *data, + int signalStartIdx, int signalEndIdx, + int &triggerIdx, int &lowerUncertainty, + int &upperUncertainty, double &snr, + OPT(Polarity) &polarity) override; + + //! Reset the picker state + void reset() override; + + + // ---------------------------------------------------------------------- + // Private Interface + // ---------------------------------------------------------------------- + private: + /** + * @brief Initialize the filter bank + * + * Creates bandpass filters for each frequency band + */ + void initFilterBank(); + + /** + * @brief Apply filter bank to input data + * + * @param data Input data + * @param n Number of samples + * @return Vector of filtered traces for each band + */ + std::vector> applyFilterBank(const double *data, int n); + + /** + * @brief Compute characteristic function for a filtered trace + * + * The characteristic function is based on the envelope of the signal + * and emphasizes the onset of seismic phases. + * + * @param filtered Filtered trace + * @return Characteristic function values + */ + CharacteristicFunction computeCharacteristicFunction(const std::vector& filtered); + + /** + * @brief Compute the integral of the maximum characteristic function + * + * @param cfs Vector of characteristic functions for all bands + * @param startIdx Start index for integration + * @param windowSize Size of the integration window in samples + * @return Integral value + */ + double computeCFIntegral(const std::vector& cfs, + int startIdx, int windowSize); + + /** + * @brief Find the pick onset using the FilterPicker algorithm + * + * @param cfs Characteristic functions for all bands + * @param threshold Detection threshold + * @param windowSize Integration window size + * @param onsetIdx Output: index of the onset + * @param snr Output: signal-to-noise ratio + * @return true if onset found, false otherwise + */ + bool findOnset(const std::vector& cfs, + double threshold, int windowSize, + int &onsetIdx, double &snr); + + /** + * @brief Estimate pick uncertainty + * + * @param cfs Characteristic functions + * @param onsetIdx Onset index + * @param lowerUnc Output: lower uncertainty in samples + * @param upperUnc Output: upper uncertainty in samples + */ + void estimateUncertainty(const std::vector& cfs, + int onsetIdx, int &lowerUnc, int &upperUnc); + + /** + * @brief Determine polarity of the onset + * + * @param data Original data + * @param onsetIdx Onset index + * @param fsamp Sampling frequency + * @return Polarity (POSITIVE, NEGATIVE, or UNDECIDABLE) + */ + Polarity determinePolarity(const double *data, int onsetIdx, double fsamp); + + // ---------------------------------------------------------------------- + // Private Members + // ---------------------------------------------------------------------- + private: + std::vector _filterBandsConfig; // Filter bank configuration + std::string _filterString; // Filter string for filterID() + + // Algorithm parameters + int _numBands; // Number of frequency bands (default: 5) + double _minFreq; // Minimum frequency in Hz (default: 0.5) + double _maxFreq; // Maximum frequency in Hz (default: 20) + int _windowSize; // Integration window size in seconds (default: 2) + double _threshold; // Detection threshold (default: 3.0) + double _thresholdOff; // Threshold for reset (default: 1.5) + bool _useAdaptiveThreshold; // Use adaptive thresholding + + // Sampling rate tracking + double _lastFsamp{0.0}; // Last sampling rate used for filter bank initialization +}; + + +} // namespace Processing +} // namespace Seiscomp + + +#endif // SEISCOMP_PROCESSING_FILTERPICKER_H diff --git a/libs/seiscomp/processing/filterpicker/filterpicker_profile_example.cfg b/libs/seiscomp/processing/filterpicker/filterpicker_profile_example.cfg new file mode 100644 index 000000000..5f3282c0b --- /dev/null +++ b/libs/seiscomp/processing/filterpicker/filterpicker_profile_example.cfg @@ -0,0 +1,64 @@ +# FilterPicker Configuration Example +# Add these parameters to your scautopick profile file +# (e.g., ~/.seiscomp/etc/key/scautopick/profile_vquiet or station-specific files) + +# ============================================================================= +# FilterPicker Parameters +# ============================================================================= + +# Number of frequency bands for the filter bank (default: 5) +# More bands = better frequency coverage but more computation +picker.FilterPicker.numBands = 5 + +# Frequency range for the filter bank (default: 0.5-20 Hz) +# Adjust based on your station characteristics and target events +picker.FilterPicker.minFreq = 0.5 +picker.FilterPicker.maxFreq = 20.0 + +# Integration window size in seconds (default: 2) +# Larger windows = more stable but slower response +picker.FilterPicker.windowSize = 2 + +# Detection threshold - CF ratio required to trigger (default: 2.0) +# Lower = more sensitive, more false picks +# Higher = fewer picks, might miss weak events +picker.FilterPicker.threshold = 2.0 + +# Threshold for turning off detection (default: 1.5) +picker.FilterPicker.thresholdOff = 1.5 + +# Enable adaptive thresholding (default: true) +# Adapts threshold based on noise level +picker.FilterPicker.adaptiveThreshold = true + +# Minimum SNR required to accept a pick (default: 2.0) +# Lower = more picks accepted, including noisy ones +# Higher = only high-quality picks accepted +picker.FilterPicker.minSNR = 2.0 + +# Noise window start relative to trigger (default: -10 s) +picker.FilterPicker.noiseBegin = -10.0 + +# Signal window start relative to trigger (default: -5 s) +picker.FilterPicker.signalBegin = -5.0 + +# Signal window end relative to trigger (default: 20 s) +picker.FilterPicker.signalEnd = 20.0 + +# ============================================================================= +# Example: Sensitive configuration for local/regional events +# ============================================================================= +# picker.FilterPicker.numBands = 5 +# picker.FilterPicker.minFreq = 1.0 +# picker.FilterPicker.maxFreq = 15.0 +# picker.FilterPicker.threshold = 1.5 +# picker.FilterPicker.minSNR = 1.5 + +# ============================================================================= +# Example: Conservative configuration for teleseismic events +# ============================================================================= +# picker.FilterPicker.numBands = 4 +# picker.FilterPicker.minFreq = 0.5 +# picker.FilterPicker.maxFreq = 5.0 +# picker.FilterPicker.threshold = 2.5 +# picker.FilterPicker.minSNR = 3.0