epclust/src/filter.cpp

   1 #include <Rcpp.h>
   2
   3 using namespace Rcpp;
   4
   5 //' filter
   6 //'
   7 //' Filter [time-]series by replacing all values by the moving average of previous, current and
   8 //' next value. Possible extensions: larger window, gaussian kernel... (but would be costly!).
   9 //' Note: border values are unchanged.
  10 //'
  11 //' @param cwt Continuous wavelets transform (in columns): a matrix of size LxD
  12 //'
  13 //' @return The filtered CWT, in a matrix of same size (LxD)
  14 // [[Rcpp::export]]
  15 RcppExport SEXP epclustFilter(SEXP cwt_)
  16 {
  17     // NOTE: C-style for better efficiency (this must be as fast as possible)
  18     int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0],
  19         D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1];
  20     double *cwt = REAL(cwt_);
  21
  22     SEXP fcwt_; //the filtered result
  23     PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
  24     double* fcwt = REAL(fcwt_); //pointer to the encapsulated vector
  25
  26     // NOTE: unused loop parameter: shifting at the end of the loop is more efficient
  27     for (int col=D-1; col>=0; col++)
  28     {
  29         double v1 = cwt[0]; //first value
  30         double ms = v1 + cwt[1] + cwt[2]; //moving sum at second value
  31         for (int i=1; i<L-2; i++)
  32         {
  33             fcwt[i] = ms / 3.; //ms / 3: moving average at current index i
  34             ms = ms - v1 + cwt[i+2]; //update ms: remove oldest items, add next
  35             v1 = cwt[i]; //update first value for next iteration
  36         }
  37
  38         // Fill a few border values not computed in the loop
  39         fcwt[0] = cwt[0];
  40         fcwt[L-2] = ms / 3.;
  41         fcwt[L-1] = cwt[L-1];
  42
  43         // Shift by L == increment column index by 1 (storage per column)
  44         cwt = cwt + L;
  45         fcwt = fcwt + L;
  46     }
  47
  48     UNPROTECT(1);
  49     return fcwt_;
  50 }