ENH: allow user to configure a trace margin to remove filter artifacts

swiss-seismological-service · Apr 16, 2022 · b461b4b · b461b4b
1 parent f225780
commit b461b4b
Show file tree

Hide file tree

Showing 6 changed files with 85 additions and 49 deletions.
diff --git a/apps/scrtdd/descriptions/scrtdd.xml b/apps/scrtdd/descriptions/scrtdd.xml
@@ -189,6 +189,9 @@
               <parameter name="filterString" type="string" default="&quot;ITAPER(1)&gt;&gt;BW_HLP(2,1,20)&quot;">
                 <description>SeisComP string based filter definition. Set to "" to disable filtering.</description>
               </parameter>
+              <parameter name="margin" type="double" default="1" unit="sec">
+                <description>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</description>
+              </parameter> 
               <parameter name="resampling" type="double" default="400" unit="Hz">
                 <description>Resample all traces at this samplig interval (hz) Set it to 0 disable resampling.</description>
               </parameter>

diff --git 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
@@ -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<Waveform::Processor> 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<Waveform::Loader> &baseLdr)
   }
   else if (_cfg.snr.minSnr > 0)
   {
-    _wfAccess.snrFilter->setAuxProcessor(
-        shared_ptr<Waveform::Processor>(new Waveform::BasicProcessor(baseLdr)));
+    _wfAccess.snrFilter->setAuxProcessor(shared_ptr<Waveform::Processor>(
+        new Waveform::BasicProcessor(baseLdr, _cfg.wfFilter.extraTraceLen)));
   }
   else
   {
-    _wfAccess.memCache->setAuxProcessor(
-        shared_ptr<Waveform::Processor>(new Waveform::BasicProcessor(baseLdr)));
+    _wfAccess.memCache->setAuxProcessor(shared_ptr<Waveform::Processor>(
+        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<Waveform::Processor> proc(new Waveform::BasicProcessor(loader));
+  shared_ptr<Waveform::Processor> 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
@@ -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
@@ -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<const Trace> 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<const Trace> 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<const Trace> BasicProcessor::get(const TimeWindow &tw,
   if (trans == Transform::TRANSVERSAL || trans == Transform::RADIAL)
   {
     shared_ptr<const Trace> trV(
-        loadWaveform(comps.names[ThreeComponents::Vertical]));
+        loadAndProcess(tw, ph, comps.names[ThreeComponents::Vertical],
+                       filterStr, resampleFreq));
     shared_ptr<const Trace> trH1(
-        loadWaveform(comps.names[ThreeComponents::FirstHorizontal]));
+        loadAndProcess(tw, ph, comps.names[ThreeComponents::FirstHorizontal],
+                       filterStr, resampleFreq));
     shared_ptr<const Trace> 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<const Trace> trH1(
-        loadWaveform(comps.names[ThreeComponents::FirstHorizontal]));
+        loadAndProcess(tw, ph, comps.names[ThreeComponents::FirstHorizontal],
+                       filterStr, resampleFreq));
     shared_ptr<const Trace> 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<const Trace> BasicProcessor::get(const TimeWindow &tw,
   return trace;
 }
 
-std::shared_ptr<Trace> BasicProcessor::process(const Trace &trace,
-                                               const std::string &filterStr,
-                                               double resampleFreq) const
+std::shared_ptr<Trace>
+BasicProcessor::loadAndProcess(const TimeWindow &tw,
+                               Catalog::Phase ph, // copied
+                               const string &channelCode,
+                               const string &filterStr,
+                               double resampleFreq) const
 {
-  shared_ptr<Trace> copy(new Trace(trace));
+  ph.channelCode = channelCode;
+
+  const TimeWindow twToLoad =
+      tw.extend(secToDur(_extraTraceLen), secToDur(_extraTraceLen));
+
+  shared_ptr<const Trace> trace = _auxLdr->get(twToLoad, ph);
+  if (!trace) return nullptr;
+
+  shared_ptr<Trace> 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
@@ -200,7 +200,9 @@ class Processor
 class BasicProcessor : public Processor
 {
 public:
-  BasicProcessor(const std::shared_ptr<Loader> &auxLdr) : _auxLdr(auxLdr) {}
+  BasicProcessor(const std::shared_ptr<Loader> &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<Trace> process(const Trace &trace,
-                                 const std::string &filterStr,
-                                 double resampleFreq) const;
+  std::shared_ptr<Trace> loadAndProcess(const TimeWindow &tw,
+                                        Catalog::Phase ph,
+                                        const std::string &channelCode,
+                                        const std::string &filterStr,
+                                        double resampleFreq) const;
+
   std::shared_ptr<Loader> _auxLdr;
+  double _extraTraceLen;
 };
 
 class MemCachedProc : public Processor