libstfinv [WP]: implement taper function

src/libs/libstfinv/stfinvfourier.cc

@@ -34,8 +34,10 @@
  *  - 14/10/2015   V1.3   new end-user usage functions
  *  - 28/06/2016   V1.4   introduce taper parameter to taper correction filter
  *                        response in the time domain
- *  - 22/07/2016   V1.5   provide separate FFT processor addressing just the
- *                        source time function correction filter
+ *  - 22/07/2016   V1.5   thof:
+ *                        - provide separate FFT processor addressing just the
+ *                          source time function correction filter
+ *                        - implement taper function
  * 
  * ============================================================================
  */
@@ -43,6 +45,7 @@
   "STFINV_STFINVFOURIER_CC   V1.5"
 
 #include <sstream>
+#include <cmath>
 #include <stfinv/stfinvfourier.h>
 #include <stfinv/stfinvfourier_summary_usage.h>
 #include <stfinv/stfinvfourier_description_usage.h>
@@ -577,7 +580,56 @@ namespace stfinv {
       Tfftengine::TAseries thestfseries=this->stfseries();
       thestfseries *= Mfftenginestf.scale_series(this->dt());
       // apply taper
-      //
+      /*
+       * Concept of application to periodic time series allowing for negative
+       * values of taper times:
+       *
+       * All samples between t2 and t3 remain unaltered. It is reasonable to
+       * process the taper starting at t3 passing t4 to t1 and then ending at
+       * t2. Behind this is the concept that time series are implicitely
+       * periodic in discrete Fourier transform.
+       */
+      // duration of time series
+      double T=this->dt()*this->nsamples();
+      // make sure taper definition is not longer than time series.
+      STFINV_assert((this->Mtt4-this->Mtt1)<T,
+                    "Taper width is larger than length of time series");
+      // set taper time values under conpect of periodic series
+      double tt3=this->Mtt3;
+      double tt4=this->Mtt4;
+      double tt1=this->Mtt1+T;
+      double tt2=this->Mtt2+T;
+      // sample to start at
+      int l3=int(std::floor(tt3/T));
+      int l2=int(std::ceil(tt2/T));
+      for (int l=l3; l<=l2; ++l)
+      {
+        // time of sample
+        double t=l*this->dt();
+        // index to series
+        int rl = l%this->nsamples();
+        if (rl < 0) { rl += this->nsamples(); }
+        STFINV_assert(((rl >= thestfseries.f(0)) && 
+                       (rl <= thestfseries.l(0))),
+                      "Index out of range. Internal programming error! "
+                      "Report as bug!")
+        // taper factor
+        double factor=1.;
+        if ( (t>=tt3) && (t<=tt4) && (tt4>tt3) )
+        {
+          factor=0.5+0.5*cos(M_PI*(t-tt3)/(tt4-tt3));
+        }
+        else if ( (t>tt4) && (t<tt1) )
+        {
+          factor=0.;
+        }
+        else if ( (t>=tt1) && (t<=tt2) && (tt2>tt1) )
+        {
+          factor=0.5-0.5*cos(M_PI*(t-tt1)/(tt2-tt1));
+        }
+        // apply to series
+        thestfseries(rl) = thestfseries(rl)*factor;
+      }
       // transform correction filter to Fourier domain
       Mfftenginestf.r2c();
       this->stfspec() *= Mfftenginestf.scale_spectrum(this->dt());