add comments, fix some things. TODO: comment tests, finish computeWerDists, test it

[epclust.git] / epclust / src / filter.cpp

diff --git a/epclust/src/filter.cpp b/epclust/src/filter.cpp
index ee24af6..5268a5f 100644 (file)
--- a/epclust/src/filter.cpp
+++ b/epclust/src/filter.cpp
@@ -4,37 +4,47 @@ using namespace Rcpp;
 
 //' filter
 //'
-//' Filter time-series
+//' Filter [time-]series by replacing all values by the moving average of previous, current and
+//' next value. Possible extensions: larger window, gaussian kernel... (but would be costly!).
+//' Note: border values are unchanged.
 //'
-//' @param cwt Continuous wavelets transform
+//' @param cwt Continuous wavelets transform (in columns): a matrix of size LxD
 //'
-//' @return The filtered CWT
+//' @return The filtered CWT, in a matrix of same size (LxD)
 // [[Rcpp::export]]
-NumericMatrix epclustFilter(NumericMatrix cwt)
+RcppExport SEXP epclustFilter(SEXP cwt_)
 {
-       int L = cwt.nrow(),
-               D = cwt.ncol();
-       NumericMatrix fcwt(L, D); //fill with 0... TODO: back to SEXP C-style?
-       double *cwt_c = cwt.begin(),
-               *fcwt_c = fcwt.begin();
+       // NOTE: C-style for better efficiency (this must be as fast as possible)
+       int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0],
+               D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1];
+       double *cwt = REAL(cwt_);
 
-       //TODO: coding style is terrible... no time for now.
-       for (int col=0; col<D; col++)
+       SEXP fcwt_; //the filtered result
+       PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
+       double* fcwt = REAL(fcwt_); //pointer to the encapsulated vector
+
+       // NOTE: unused loop parameter: shifting at the end of the loop is more efficient
+       for (int col=D-1; col>=0; col++)
        {
-               double v1 = cwt_c[0];
-               double ma = v1 + cwt[1] + cwt_c[2];
+               double v1 = cwt[0]; //first value
+               double ms = v1 + cwt[1] + cwt[2]; //moving sum at second value
                for (int i=1; i<L-2; i++)
                {
-                       fcwt_c[i] = ma / 3.;
-                       ma = ma - v1 + cwt_c[i+2];
-                       v1 = cwt_c[i];
+                       fcwt[i] = ms / 3.; //ms / 3: moving average at current index i
+                       ms = ms - v1 + cwt[i+2]; //update ms: remove oldest items, add next
+                       v1 = cwt[i]; //update first value for next iteration
                }
-               fcwt_c[0] = cwt_c[0];
-               fcwt_c[L-2] = ma / 3.;
-               fcwt_c[L-1] = cwt_c[L-1];
-               cwt_c = cwt_c + L;
-               fcwt_c = fcwt_c + L;
+
+               // Fill a few border values not computed in the loop
+               fcwt[0] = cwt[0];
+               fcwt[L-2] = ms / 3.;
+               fcwt[L-1] = cwt[L-1];
+
+               // Shift by L == increment column index by 1 (storage per column)
+               cwt = cwt + L;
+               fcwt = fcwt + L;
        }
 
-       return fcwt;
+       UNPROTECT(1);
+       return fcwt_;
 }