| 1 | #include <Rcpp.h> |
| 2 | |
| 3 | #include <R.h> // Rprintf() |
| 4 | //#include <R_ext/Utils.h> // user interrupts |
| 5 | #include <Rdefines.h> |
| 6 | #include <Rinternals.h> |
| 7 | |
| 8 | #include <stdio.h> |
| 9 | |
| 10 | using namespace Rcpp; |
| 11 | |
| 12 | //' filter |
| 13 | //' |
| 14 | //' Filter time-series |
| 15 | //' |
| 16 | //' @param cwt Continuous wavelets transform |
| 17 | //' |
| 18 | //' @return The filtered CWT |
| 19 | // [[Rcpp::export]] |
| 20 | RcppExport SEXP epclustFilter(SEXP cwt_) |
| 21 | { |
| 22 | int L = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[0], |
| 23 | D = INTEGER(Rf_getAttrib(cwt_, R_DimSymbol))[1]; |
| 24 | double *cwt = REAL(cwt_); |
| 25 | SEXP fcwt_; |
| 26 | PROTECT(fcwt_ = Rf_allocMatrix(REALSXP, L, D)); |
| 27 | double* fcwt = REAL(fcwt_); //(double*)malloc(L*D*sizeof(double)); |
| 28 | |
| 29 | //TODO: coding style is terrible... no time for now. |
| 30 | for (int col=0; col<D; col++) |
| 31 | { |
| 32 | double v1 = cwt[0]; |
| 33 | double ma = v1 + cwt[1] + cwt[2]; |
| 34 | for (int i=1; i<L-2; i++) |
| 35 | { |
| 36 | fcwt[i] = ma / 3.; |
| 37 | ma = ma - v1 + cwt[i+2]; |
| 38 | v1 = cwt[i]; |
| 39 | } |
| 40 | fcwt[0] = cwt[0]; |
| 41 | fcwt[L-2] = ma / 3.; |
| 42 | fcwt[L-1] = cwt[L-1]; |
| 43 | cwt = cwt + L; |
| 44 | fcwt = fcwt + L; |
| 45 | } |
| 46 | |
| 47 | // REAL(fcwt_) = fcwt; |
| 48 | UNPROTECT(1); |
| 49 | |
| 50 | return fcwt_; |
| 51 | } |