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 | { | |
15 | int L = INTEGER(getAttrib(cwt_, R_DimSymbol))[0], | |
16 | D = INTEGER(getAttrib(cwt_, 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 *cwt = REAL(cwt_), | |
22 | cs, //current sum (max: w terms) | |
23 | left; //leftward term in the current moving sum | |
24 | ||
25 | SEXP fcwt_; //the filtered result | |
26 | PROTECT(fcwt_ = allocMatrix(REALSXP, L, D)); | |
27 | double* fcwt = REAL(fcwt_); //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 = cwt[0]; | |
35 | cs = 0.; | |
36 | for (i=half_w; i>=0; i--) | |
37 | cs += cwt[i]; | |
38 | ||
39 | // Left border | |
40 | for (i=0; i<half_w; i++) | |
41 | { | |
42 | fcwt[i] = cs / nt; //(partial) moving average at current index i | |
43 | cs += cwt[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 | fcwt[i] = cs / w; //moving average at current index i | |
51 | cs = ms - left + cwt[i+half_w+1]; //remove oldest items, add next | |
52 | left = cwt[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 | fcwt[i] = cs / nt; //(partial) moving average at current index i | |
59 | cs -= cwt[i-half_w]; | |
60 | nt--; | |
61 | } | |
62 | ||
63 | // Shift by L == increment column index by 1 (storage per column) | |
64 | cwt = cwt + L; | |
65 | fcwt = fcwt + L; | |
66 | } | |
67 | ||
68 | UNPROTECT(1); | |
69 | return fcwt_; | |
70 | } |