Commit | Line | Data |
---|---|---|
363ae134 BA |
1 | #include <Rcpp.h> |
2 | ||
3 | using namespace Rcpp; | |
4 | ||
5 | //' filter | |
6 | //' | |
d9bb53c5 BA |
7 | //' Filter [time-]series by replacing all values by the moving average of previous, current and |
8 | //' next value. Possible extensions: larger window, gaussian kernel... (but would be costly!). | |
9 | //' Note: border values are unchanged. | |
363ae134 | 10 | //' |
d9bb53c5 | 11 | //' @param cwt Continuous wavelets transform (in columns): a matrix of size LxD |
363ae134 | 12 | //' |
d9bb53c5 | 13 | //' @return The filtered CWT, in a matrix of same size (LxD) |
363ae134 | 14 | // [[Rcpp::export]] |
a52836b2 | 15 | RcppExport SEXP filterMA(SEXP cwt_) |
363ae134 | 16 | { |
d9bb53c5 | 17 | // NOTE: C-style for better efficiency (this must be as fast as possible) |
0fe757f7 BA |
18 | int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0], |
19 | D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1]; | |
20 | double *cwt = REAL(cwt_); | |
d9bb53c5 BA |
21 | |
22 | SEXP fcwt_; //the filtered result | |
0fe757f7 | 23 | PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D)); |
d9bb53c5 | 24 | double* fcwt = REAL(fcwt_); //pointer to the encapsulated vector |
363ae134 | 25 | |
d9bb53c5 | 26 | // NOTE: unused loop parameter: shifting at the end of the loop is more efficient |
a52836b2 | 27 | for (int col=D-1; col>=0; col--) |
363ae134 | 28 | { |
d9bb53c5 BA |
29 | double v1 = cwt[0]; //first value |
30 | double ms = v1 + cwt[1] + cwt[2]; //moving sum at second value | |
363ae134 BA |
31 | for (int i=1; i<L-2; i++) |
32 | { | |
d9bb53c5 BA |
33 | fcwt[i] = ms / 3.; //ms / 3: moving average at current index i |
34 | ms = ms - v1 + cwt[i+2]; //update ms: remove oldest items, add next | |
35 | v1 = cwt[i]; //update first value for next iteration | |
363ae134 | 36 | } |
d9bb53c5 BA |
37 | |
38 | // Fill a few border values not computed in the loop | |
0fe757f7 | 39 | fcwt[0] = cwt[0]; |
d9bb53c5 | 40 | fcwt[L-2] = ms / 3.; |
0fe757f7 | 41 | fcwt[L-1] = cwt[L-1]; |
d9bb53c5 BA |
42 | |
43 | // Shift by L == increment column index by 1 (storage per column) | |
0fe757f7 BA |
44 | cwt = cwt + L; |
45 | fcwt = fcwt + L; | |
363ae134 BA |
46 | } |
47 | ||
0fe757f7 | 48 | UNPROTECT(1); |
0fe757f7 | 49 | return fcwt_; |
363ae134 | 50 | } |