X-Git-Url: https://git.auder.net/?p=epclust.git;a=blobdiff_plain;f=epclust%2Fsrc%2Ffilter.cpp;h=5268a5fefbb187070c3072ce16bff9aeddc385fb;hp=cb5a8a0b5d58e553770401accff99438fe0af73f;hb=d9bb53c5e1392018bf67f92140edb10137f3423c;hpb=9f05a4a0b703deffd7bdb9cd99b0aaa2246a5c83

diff --git a/epclust/src/filter.cpp b/epclust/src/filter.cpp
index cb5a8a0..5268a5f 100644
--- a/epclust/src/filter.cpp
+++ b/epclust/src/filter.cpp
@@ -1,51 +1,50 @@
 #include <Rcpp.h>
 
-#include <R.h>           //  Rprintf()
-//#include <R_ext/Utils.h> //  user interrupts
-#include <Rdefines.h>
-#include <Rinternals.h>
-
-#include <stdio.h>
-
 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]]
 RcppExport SEXP epclustFilter(SEXP cwt_)
 {
+	// 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_);
-	SEXP fcwt_;
+
+	SEXP fcwt_; //the filtered result
 	PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D));
-		double* fcwt = REAL(fcwt_); //(double*)malloc(L*D*sizeof(double));
+	double* fcwt = REAL(fcwt_); //pointer to the encapsulated vector
 
-	//TODO: coding style is terrible... no time for now.
-	for (int col=0; col<D; col++)
+	// 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[0];
-		double ma = v1 + cwt[1] + cwt[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[i] = ma / 3.;
-			ma = ma - v1 + cwt[i+2];
-			v1 = cwt[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
 		}
+
+		// Fill a few border values not computed in the loop
 		fcwt[0] = cwt[0];
-		fcwt[L-2] = ma / 3.;
+		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;
 	}
 
-//	REAL(fcwt_) = fcwt;
 	UNPROTECT(1);
-
 	return fcwt_;
 }