From b461b4be8fcf477f38701e458cd0a1baa98060ae Mon Sep 17 00:00:00 2001 From: Luca Scarabello Date: Sat, 16 Apr 2022 16:23:37 +0200 Subject: [PATCH] ENH: allow user to configure a trace margin to remove filter artifacts --- apps/scrtdd/descriptions/scrtdd.xml | 3 ++ apps/scrtdd/rtdd.cpp | 11 +++- libs/hdd/dd.cpp | 22 ++++---- libs/hdd/dd.h | 4 ++ libs/hdd/waveform.cpp | 80 ++++++++++++++++------------- libs/hdd/waveform.h | 14 +++-- 6 files changed, 85 insertions(+), 49 deletions(-) diff --git a/apps/scrtdd/descriptions/scrtdd.xml b/apps/scrtdd/descriptions/scrtdd.xml index 2d97bde2..30b07769 100644 --- a/apps/scrtdd/descriptions/scrtdd.xml +++ b/apps/scrtdd/descriptions/scrtdd.xml @@ -189,6 +189,9 @@ SeisComP string based filter definition. Set to "" to disable filtering. + + Extra seconds to add to the waveform ends before applying the filter. Useful to initialize the filter and to discard potential filter artifacts at the beginning and end of trace + Resample all traces at this samplig interval (hz) Set it to 0 disable resampling. diff --git a/apps/scrtdd/rtdd.cpp b/apps/scrtdd/rtdd.cpp index ca7dc2fe..14484182 100644 --- a/apps/scrtdd/rtdd.cpp +++ b/apps/scrtdd/rtdd.cpp @@ -837,6 +837,14 @@ bool RTDD::validateParameters() prof->ddCfg.wfFilter.filterStr = "ITAPER(1)>>BW_HLP(2,1,20)"; } try + { + prof->ddCfg.wfFilter.extraTraceLen = configGetDouble(prefix + "margin"); + } + catch (...) + { + prof->ddCfg.wfFilter.extraTraceLen = 1.0; + } + try { prof->ddCfg.wfFilter.resampleFreq = configGetDouble(prefix + "resampling"); @@ -1031,7 +1039,8 @@ bool RTDD::validateParameters() } prof->ddCfg.recordStreamURL = recordStreamURL(); - prof->ddCfg.diskTraceMinLen = configGetDouble("performance.cachedWaveformLength"); + prof->ddCfg.diskTraceMinLen = + configGetDouble("performance.cachedWaveformLength"); // no reason to make those configurable prof->singleEventClustering.minWeight = 0; diff --git a/libs/hdd/dd.cpp b/libs/hdd/dd.cpp index f5286d23..e378ed93 100644 --- a/libs/hdd/dd.cpp +++ b/libs/hdd/dd.cpp @@ -141,13 +141,13 @@ void DD::createWaveformCache() { _wfAccess.diskCache.reset( new Waveform::DiskCachedLoader(currLdr, _cacheDir)); - _wfAccess.extraLen.reset( - new Waveform::ExtraLenLoader(_wfAccess.diskCache, _cfg.diskTraceMinLen)); + _wfAccess.extraLen.reset(new Waveform::ExtraLenLoader( + _wfAccess.diskCache, _cfg.diskTraceMinLen)); currLdr = _wfAccess.extraLen; } shared_ptr currProc( - new Waveform::BasicProcessor(currLdr)); + new Waveform::BasicProcessor(currLdr, _cfg.wfFilter.extraTraceLen)); if (_cfg.snr.minSnr > 0) { @@ -168,13 +168,13 @@ void DD::replaceWaveformCacheLoader(const shared_ptr &baseLdr) } else if (_cfg.snr.minSnr > 0) { - _wfAccess.snrFilter->setAuxProcessor( - shared_ptr(new Waveform::BasicProcessor(baseLdr))); + _wfAccess.snrFilter->setAuxProcessor(shared_ptr( + new Waveform::BasicProcessor(baseLdr, _cfg.wfFilter.extraTraceLen))); } else { - _wfAccess.memCache->setAuxProcessor( - shared_ptr(new Waveform::BasicProcessor(baseLdr))); + _wfAccess.memCache->setAuxProcessor(shared_ptr( + new Waveform::BasicProcessor(baseLdr, _cfg.wfFilter.extraTraceLen))); } } @@ -1735,7 +1735,8 @@ DD::preloadNonCatalogWaveforms(Catalog &catalog, if (_useCatalogWaveformDiskCache && _waveformCacheAll) { diskLoader.reset(new Waveform::DiskCachedLoader(batchLoader, _tmpCacheDir)); - loader.reset(new Waveform::ExtraLenLoader(diskLoader, _cfg.diskTraceMinLen)); + loader.reset( + new Waveform::ExtraLenLoader(diskLoader, _cfg.diskTraceMinLen)); } // Load a large enough waveform to be able to check the SNR after @@ -1747,13 +1748,16 @@ DD::preloadNonCatalogWaveforms(Catalog &catalog, double extraBefore = std::min({_cfg.snr.noiseStart, _cfg.snr.signalStart}) - maxDelay; extraBefore = extraBefore > 0 ? 0 : std::abs(extraBefore); + extraBefore += _cfg.wfFilter.extraTraceLen; double extraAfter = std::max({_cfg.snr.noiseEnd, _cfg.snr.signalEnd}) + maxDelay; extraAfter = extraAfter < 0 ? 0 : extraAfter; + extraAfter += _cfg.wfFilter.extraTraceLen; loader.reset(new Waveform::ExtraLenLoader(loader, extraBefore, extraAfter)); } - shared_ptr proc(new Waveform::BasicProcessor(loader)); + shared_ptr proc( + new Waveform::BasicProcessor(loader, _cfg.wfFilter.extraTraceLen)); // // loop through reference event phases diff --git a/libs/hdd/dd.h b/libs/hdd/dd.h index bd061ec2..69be7d25 100644 --- a/libs/hdd/dd.h +++ b/libs/hdd/dd.h @@ -63,6 +63,10 @@ struct Config { std::string filterStr = "ITAPER(1)>>BW_HLP(2,1,20)"; // "" -> no filtering double resampleFreq = 400; // 0 -> no resampling + // Extra waveform to load before and after the required window. This extra + // len is useful to initialize the filter and to discard potential filter + // artifacts at the beginning and end of trace + double extraTraceLen = 1; // seconds } wfFilter; struct diff --git a/libs/hdd/waveform.cpp b/libs/hdd/waveform.cpp index a363a7b4..e02bb605 100644 --- a/libs/hdd/waveform.cpp +++ b/libs/hdd/waveform.cpp @@ -156,20 +156,9 @@ void BatchLoader::load() TimeWindow ExtraLenLoader::traceTimeWindowToLoad(const TimeWindow &neededTW, const UTCTime &pickTime) const { - TimeWindow twToLoad = neededTW; - - // Make sure to load at least `_traceMinLenh` seconds of the waveform. - if (_beforePickLen > 0 || _afterPickLen > 0) - { - const UTCTime::duration additionalTimeBefore = secToDur(_beforePickLen); - const UTCTime::duration additionalTimeAfter = secToDur(_afterPickLen); - - if (twToLoad.startTime() > pickTime - additionalTimeBefore) - twToLoad.setStartTime(pickTime - additionalTimeBefore); - - if (twToLoad.endTime() < pickTime + additionalTimeAfter) - twToLoad.setEndTime(pickTime + additionalTimeAfter); - } + TimeWindow extraTW(pickTime - secToDur(_beforePickLen), + pickTime + secToDur(_afterPickLen)); + TimeWindow twToLoad = extraTW.merge(neededTW); // // round the start/end time to the nearest second @@ -299,18 +288,9 @@ shared_ptr BasicProcessor::get(const TimeWindow &tw, double resampleFreq, Transform trans) { - auto loadWaveform = [this, &tw, &ph, &filterStr, - resampleFreq](const string &channelCode) { - Catalog::Phase copy(ph); - copy.channelCode = channelCode; - shared_ptr trace = _auxLdr->get(tw, copy); - if (trace) trace = process(*trace, filterStr, resampleFreq); - return trace; - }; - if (trans == Transform::NONE) { - return loadWaveform(ph.channelCode); + return loadAndProcess(tw, ph, ph.channelCode, filterStr, resampleFreq); } // @@ -338,20 +318,25 @@ shared_ptr BasicProcessor::get(const TimeWindow &tw, if (trans == Transform::TRANSVERSAL || trans == Transform::RADIAL) { shared_ptr trV( - loadWaveform(comps.names[ThreeComponents::Vertical])); + loadAndProcess(tw, ph, comps.names[ThreeComponents::Vertical], + filterStr, resampleFreq)); shared_ptr trH1( - loadWaveform(comps.names[ThreeComponents::FirstHorizontal])); + loadAndProcess(tw, ph, comps.names[ThreeComponents::FirstHorizontal], + filterStr, resampleFreq)); shared_ptr trH2( - loadWaveform(comps.names[ThreeComponents::SecondHorizontal])); + loadAndProcess(tw, ph, comps.names[ThreeComponents::SecondHorizontal], + filterStr, resampleFreq)); if (trH1 && trH2 && trV) trace = transformRT(tw, ph, ev, sta, comps, *trV, *trH1, *trH2, trans); } else if (trans == Transform::L2) { shared_ptr trH1( - loadWaveform(comps.names[ThreeComponents::FirstHorizontal])); + loadAndProcess(tw, ph, comps.names[ThreeComponents::FirstHorizontal], + filterStr, resampleFreq)); shared_ptr trH2( - loadWaveform(comps.names[ThreeComponents::SecondHorizontal])); + loadAndProcess(tw, ph, comps.names[ThreeComponents::SecondHorizontal], + filterStr, resampleFreq)); if (trH1 && trH2) trace = transformL2(tw, ph, comps, *trH1, *trH2); } @@ -368,19 +353,44 @@ shared_ptr BasicProcessor::get(const TimeWindow &tw, return trace; } -std::shared_ptr BasicProcessor::process(const Trace &trace, - const std::string &filterStr, - double resampleFreq) const +std::shared_ptr +BasicProcessor::loadAndProcess(const TimeWindow &tw, + Catalog::Phase ph, // copied + const string &channelCode, + const string &filterStr, + double resampleFreq) const { - shared_ptr copy(new Trace(trace)); + ph.channelCode = channelCode; + + const TimeWindow twToLoad = + tw.extend(secToDur(_extraTraceLen), secToDur(_extraTraceLen)); + + shared_ptr trace = _auxLdr->get(twToLoad, ph); + if (!trace) return nullptr; + + shared_ptr copy(new Trace(*trace)); try { filter(*copy, true, filterStr, resampleFreq); } catch (exception &e) { - logWarning("Errow while filtering waveform: %s", e.what()); - copy.reset(); + logDebug("Errow while filtering waveform: %s", e.what()); + return nullptr; + } + + if (!copy->slice(tw)) + { + logDebug("Error while processing phase data '%s': cannot slice trace from " + "%s length %.2f sec. Trace data from %s length %.2f sec, samples " + "%zu sampfreq %f", + string(ph).c_str(), + HDD::UTCClock::toString(tw.startTime()).c_str(), + HDD::durToSec(tw.length()), + HDD::UTCClock::toString(copy->startTime()).c_str(), + HDD::durToSec(copy->timeWindow().length()), copy->sampleCount(), + copy->samplingFrequency()); + return nullptr; } return copy; } diff --git a/libs/hdd/waveform.h b/libs/hdd/waveform.h index 0be4779b..68a1d684 100644 --- a/libs/hdd/waveform.h +++ b/libs/hdd/waveform.h @@ -200,7 +200,9 @@ class Processor class BasicProcessor : public Processor { public: - BasicProcessor(const std::shared_ptr &auxLdr) : _auxLdr(auxLdr) {} + BasicProcessor(const std::shared_ptr &auxLdr, double extraTraceLen) + : _auxLdr(auxLdr), _extraTraceLen(extraTraceLen) + {} virtual ~BasicProcessor() = default; @@ -216,10 +218,14 @@ class BasicProcessor : public Processor Transform trans) override; private: - std::shared_ptr process(const Trace &trace, - const std::string &filterStr, - double resampleFreq) const; + std::shared_ptr loadAndProcess(const TimeWindow &tw, + Catalog::Phase ph, + const std::string &channelCode, + const std::string &filterStr, + double resampleFreq) const; + std::shared_ptr _auxLdr; + double _extraTraceLen; }; class MemCachedProc : public Processor