Commit | Line | Data |
---|---|---|
363ae134 BA |
1 | #include <Rcpp.h> |
2 | ||
0fe757f7 BA |
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 | ||
363ae134 BA |
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]] | |
0fe757f7 | 20 | RcppExport SEXP epclustFilter(SEXP cwt_) |
363ae134 | 21 | { |
0fe757f7 BA |
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)); | |
363ae134 BA |
28 | |
29 | //TODO: coding style is terrible... no time for now. | |
30 | for (int col=0; col<D; col++) | |
31 | { | |
0fe757f7 BA |
32 | double v1 = cwt[0]; |
33 | double ma = v1 + cwt[1] + cwt[2]; | |
363ae134 BA |
34 | for (int i=1; i<L-2; i++) |
35 | { | |
0fe757f7 BA |
36 | fcwt[i] = ma / 3.; |
37 | ma = ma - v1 + cwt[i+2]; | |
38 | v1 = cwt[i]; | |
363ae134 | 39 | } |
0fe757f7 BA |
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; | |
363ae134 BA |
45 | } |
46 | ||
0fe757f7 BA |
47 | // REAL(fcwt_) = fcwt; |
48 | UNPROTECT(1); | |
49 | ||
50 | return fcwt_; | |
363ae134 | 51 | } |