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