8ece9c3f439dc7eaf2e23f3962e816fec42e829a
[epclust.git] / epclust / src / filter.c
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 }