| 1 | #include <R.h> |
| 2 | #include <Rdefines.h> |
| 3 | |
| 4 | // filterMA |
| 5 | // |
| 6 | // Filter [time-]series by replacing all values by the moving average of values |
| 7 | // centered around current one. Border values are averaged with available data. |
| 8 | // |
| 9 | // @param M_ A real matrix of size LxD |
| 10 | // @param w_ The (odd) number of values to average |
| 11 | // |
| 12 | // @return The filtered matrix, of same size as the input |
| 13 | SEXP filterMA(SEXP M_, SEXP w_) |
| 14 | { |
| 15 | int L = INTEGER(getAttrib(M_, R_DimSymbol))[0], |
| 16 | D = INTEGER(getAttrib(M_, R_DimSymbol))[1], |
| 17 | w = INTEGER_VALUE(w_), |
| 18 | half_w = (w-1)/2, |
| 19 | i, |
| 20 | nt; //number of terms in the current sum (max: w) |
| 21 | double *M = REAL(M_), |
| 22 | cs, //current sum (max: w terms) |
| 23 | left; //leftward term in the current moving sum |
| 24 | |
| 25 | SEXP fM_; //the filtered result |
| 26 | PROTECT(fM_ = allocMatrix(REALSXP, L, D)); |
| 27 | double* fM = REAL(fM_); //pointer to the encapsulated vector |
| 28 | |
| 29 | // NOTE: unused loop parameter: shifting at the end of the loop is more efficient |
| 30 | for (int col=D-1; col>=0; col--) |
| 31 | { |
| 32 | // Initialization |
| 33 | nt = half_w + 1; |
| 34 | left = M[0]; |
| 35 | cs = 0.; |
| 36 | for (i=half_w; i>=0; i--) |
| 37 | cs += M[i]; |
| 38 | |
| 39 | // Left border |
| 40 | for (i=0; i<half_w; i++) |
| 41 | { |
| 42 | fM[i] = cs / nt; //(partial) moving average at current index i |
| 43 | cs += M[i+half_w+1]; |
| 44 | nt++; |
| 45 | } |
| 46 | |
| 47 | // Middle: now nt == w, i == half_w |
| 48 | for (; i<L-half_w-1; i++) |
| 49 | { |
| 50 | fM[i] = cs / w; //moving average at current index i |
| 51 | cs = cs - left + M[i+half_w+1]; //remove oldest items, add next |
| 52 | left = M[i-half_w+1]; //update first value for next iteration |
| 53 | } |
| 54 | |
| 55 | // (Last "fully averaged" index +) Right border |
| 56 | for (; i<L; i++) |
| 57 | { |
| 58 | fM[i] = cs / nt; //(partial) moving average at current index i |
| 59 | cs -= M[i-half_w]; |
| 60 | nt--; |
| 61 | } |
| 62 | |
| 63 | // Shift by L == increment column index by 1 (storage per column) |
| 64 | M = M + L; |
| 65 | fM = fM + L; |
| 66 | } |
| 67 | |
| 68 | UNPROTECT(1); |
| 69 | return fM_; |
| 70 | } |