use Rcpp; ongoing debug for parallel synchrones computation

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

diff --git a/epclust/src/filter.cpp b/epclust/src/filter.cpp
new file mode 100644 (file)
index 0000000..ee24af6
--- /dev/null
+++ b/epclust/src/filter.cpp
@@ -0,0 +1,40 @@
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+//' filter
+//'
+//' Filter time-series
+//'
+//' @param cwt Continuous wavelets transform
+//'
+//' @return The filtered CWT
+// [[Rcpp::export]]
+NumericMatrix epclustFilter(NumericMatrix 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();
+
+       //TODO: coding style is terrible... no time for now.
+       for (int col=0; col<D; col++)
+       {
+               double v1 = cwt_c[0];
+               double ma = v1 + cwt[1] + cwt_c[2];
+               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_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;
+       }
+
+       return fcwt;
+}