epclust/src/filter.cpp

   1 #include <Rcpp.h>
   2
   3 #include <R.h>           //  Rprintf()
   4 //#include <R_ext/Utils.h> //  user interrupts
   5 #include <Rdefines.h>
   6 #include <Rinternals.h>
   7
   8 #include <stdio.h>
   9
  10 using namespace Rcpp;
  11
  12 //' filter
  13 //'
  14 //' Filter time-series
  15 //'
  16 //' @param cwt Continuous wavelets transform
  17 //'
  18 //' @return The filtered CWT
  19 // [[Rcpp::export]]
  20 RcppExport SEXP epclustFilter(SEXP cwt_)
  21 {
  22     int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0],
  23         D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1];
  24     double *cwt = REAL(cwt_);
  25     SEXP fcwt_;
  26     PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
  27         double* fcwt = REAL(fcwt_); //(double*)malloc(L*D*sizeof(double));
  28
  29     //TODO: coding style is terrible... no time for now.
  30     for (int col=0; col<D; col++)
  31     {
  32         double v1 = cwt[0];
  33         double ma = v1 + cwt[1] + cwt[2];
  34         for (int i=1; i<L-2; i++)
  35         {
  36             fcwt[i] = ma / 3.;
  37             ma = ma - v1 + cwt[i+2];
  38             v1 = cwt[i];
  39         }
  40         fcwt[0] = cwt[0];
  41         fcwt[L-2] = ma / 3.;
  42         fcwt[L-1] = cwt[L-1];
  43         cwt = cwt + L;
  44         fcwt = fcwt + L;
  45     }
  46
  47 //  REAL(fcwt_) = fcwt;
  48     UNPROTECT(1);
  49
  50     return fcwt_;
  51 }