epclust/src/filter.cpp

   1 #include <Rcpp.h>
   2
   3 using namespace Rcpp;
   4
   5 //' filter
   6 //'
   7 //' Filter time-series
   8 //'
   9 //' @param cwt Continuous wavelets transform
  10 //'
  11 //' @return The filtered CWT
  12 // [[Rcpp::export]]
  13 NumericMatrix epclustFilter(NumericMatrix cwt)
  14 {
  15     int L = cwt.nrow(),
  16         D = cwt.ncol();
  17     NumericMatrix fcwt(L, D); //fill with 0... TODO: back to SEXP C-style?
  18     double *cwt_c = cwt.begin(),
  19         *fcwt_c = fcwt.begin();
  20
  21     //TODO: coding style is terrible... no time for now.
  22     for (int col=0; col<D; col++)
  23     {
  24         double v1 = cwt_c[0];
  25         double ma = v1 + cwt[1] + cwt_c[2];
  26         for (int i=1; i<L-2; i++)
  27         {
  28             fcwt_c[i] = ma / 3.;
  29             ma = ma - v1 + cwt_c[i+2];
  30             v1 = cwt_c[i];
  31         }
  32         fcwt_c[0] = cwt_c[0];
  33         fcwt_c[L-2] = ma / 3.;
  34         fcwt_c[L-1] = cwt_c[L-1];
  35         cwt_c = cwt_c + L;
  36         fcwt_c = fcwt_c + L;
  37     }
  38
  39     return fcwt;
  40 }