5268a5fefbb187070c3072ce16bff9aeddc385fb
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.
11 //' @param cwt Continuous wavelets transform (in columns): a matrix of size LxD
13 //' @return The filtered CWT, in a matrix of same size (LxD)
15 RcppExport SEXP
epclustFilter(SEXP cwt_
)
17 // NOTE: C-style for better efficiency (this must be as fast as possible)
18 int L
= INTEGER(Rf_getAttrib(cwt_
, R_DimSymbol
))[0],
19 D
= INTEGER(Rf_getAttrib(cwt_
, R_DimSymbol
))[1];
20 double *cwt
= REAL(cwt_
);
22 SEXP fcwt_
; //the filtered result
23 PROTECT(fcwt_
= Rf_allocMatrix(REALSXP
, L
, D
));
24 double* fcwt
= REAL(fcwt_
); //pointer to the encapsulated vector
26 // NOTE: unused loop parameter: shifting at the end of the loop is more efficient
27 for (int col
=D
-1; col
>=0; col
++)
29 double v1
= cwt
[0]; //first value
30 double ms
= v1
+ cwt
[1] + cwt
[2]; //moving sum at second value
31 for (int i
=1; i
<L
-2; i
++)
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
38 // Fill a few border values not computed in the loop
43 // Shift by L == increment column index by 1 (storage per column)