projects
/
epclust.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
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
cb5a8a0
..
5268a5f
100644
(file)
--- a/
epclust/src/filter.cpp
+++ b/
epclust/src/filter.cpp
@@
-1,51
+1,50
@@
#include <Rcpp.h>
#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
//'
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_)
{
// [[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_);
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));
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 m
a = v1 + cwt[1] + cwt[2];
+ double v1 = cwt[0];
//first value
+ double m
s = v1 + cwt[1] + cwt[2]; //moving sum at second value
for (int i=1; i<L-2; i++)
{
for (int i=1; i<L-2; i++)
{
- fcwt[i] = m
a / 3.;
- m
a = ma - v1 + cwt[i+2];
- v1 = cwt[i];
+ fcwt[i] = m
s / 3.; //ms / 3: moving average at current index i
+ m
s = 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[0] = cwt[0];
- fcwt[L-2] = m
a
/ 3.;
+ fcwt[L-2] = m
s
/ 3.;
fcwt[L-1] = cwt[L-1];
fcwt[L-1] = cwt[L-1];
+
+ // Shift by L == increment column index by 1 (storage per column)
cwt = cwt + L;
fcwt = fcwt + L;
}
cwt = cwt + L;
fcwt = fcwt + L;
}
-// REAL(fcwt_) = fcwt;
UNPROTECT(1);
UNPROTECT(1);
-
return fcwt_;
}
return fcwt_;
}