| 1 | #include <stdlib.h> |
| 2 | #include <R.h> |
| 3 | #include <Rinternals.h> |
| 4 | #include <Rmath.h> |
| 5 | |
| 6 | #include <stdio.h> |
| 7 | |
| 8 | SEXP filter(SEXP cwt_) |
| 9 | { |
| 10 | int L = INTEGER(getAttrib(cwt_, R_DimSymbol))[0], |
| 11 | D = INTEGER(getAttrib(cwt_, R_DimSymbol))[1]; |
| 12 | double *cwt = REAL(cwt_); |
| 13 | SEXP R_fcwt; |
| 14 | PROTECT(R_fcwt = allocMatrix(REALSXP, L, D)); |
| 15 | double* fcwt = REAL(R_fcwt); |
| 16 | |
| 17 | //TODO: coding style is terrible... no time for now. |
| 18 | for (int col=0; col<D; col++) |
| 19 | { |
| 20 | double v1 = cwt[0]; |
| 21 | double ma = v1 + cwt[1] + cwt[2]; |
| 22 | for (int i=1; i<L-2; i++) |
| 23 | { |
| 24 | fcwt[i] = ma / 3.; |
| 25 | ma = ma - v1 + cwt[i+2]; |
| 26 | v1 = cwt[i]; |
| 27 | } |
| 28 | fcwt[0] = cwt[0]; |
| 29 | fcwt[L-2] = ma / 3.; |
| 30 | fcwt[L-1] = cwt[L-1]; |
| 31 | cwt = cwt + L; |
| 32 | fcwt = fcwt + L; |
| 33 | } |
| 34 | |
| 35 | UNPROTECT(1); |
| 36 | return R_fcwt; |
| 37 | } |