Commit | Line | Data |
---|---|---|
40f12a2f BA |
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 | { | |
282342ba BA |
15 | int L = INTEGER(getAttrib(M_, R_DimSymbol))[0], |
16 | D = INTEGER(getAttrib(M_, R_DimSymbol))[1], | |
40f12a2f BA |
17 | w = INTEGER_VALUE(w_), |
18 | half_w = (w-1)/2, | |
19 | i, | |
20 | nt; //number of terms in the current sum (max: w) | |
282342ba | 21 | double *M = REAL(M_), |
40f12a2f BA |
22 | cs, //current sum (max: w terms) |
23 | left; //leftward term in the current moving sum | |
24 | ||
282342ba BA |
25 | SEXP fM_; //the filtered result |
26 | PROTECT(fM_ = allocMatrix(REALSXP, L, D)); | |
27 | double* fM = REAL(fM_); //pointer to the encapsulated vector | |
40f12a2f BA |
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; | |
282342ba | 34 | left = M[0]; |
40f12a2f BA |
35 | cs = 0.; |
36 | for (i=half_w; i>=0; i--) | |
282342ba | 37 | cs += M[i]; |
40f12a2f BA |
38 | |
39 | // Left border | |
40 | for (i=0; i<half_w; i++) | |
41 | { | |
282342ba BA |
42 | fM[i] = cs / nt; //(partial) moving average at current index i |
43 | cs += M[i+half_w+1]; | |
40f12a2f BA |
44 | nt++; |
45 | } | |
46 | ||
47 | // Middle: now nt == w, i == half_w | |
48 | for (; i<L-half_w-1; i++) | |
49 | { | |
282342ba BA |
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 | |
40f12a2f BA |
53 | } |
54 | ||
55 | // (Last "fully averaged" index +) Right border | |
56 | for (; i<L; i++) | |
57 | { | |
282342ba BA |
58 | fM[i] = cs / nt; //(partial) moving average at current index i |
59 | cs -= M[i-half_w]; | |
40f12a2f BA |
60 | nt--; |
61 | } | |
62 | ||
63 | // Shift by L == increment column index by 1 (storage per column) | |
282342ba BA |
64 | M = M + L; |
65 | fM = fM + L; | |
40f12a2f BA |
66 | } |
67 | ||
68 | UNPROTECT(1); | |
282342ba | 69 | return fM_; |
40f12a2f | 70 | } |