-#include <stdlib.h>
#include <R.h>
-#include <Rinternals.h>
-#include <Rmath.h>
+#include <Rdefines.h>
-#include <stdio.h>
-
-SEXP filter(SEXP cwt_)
+// filterMA
+//
+// Filter [time-]series by replacing all values by the moving average of values
+// centered around current one. Border values are averaged with available data.
+//
+// @param M_ A real matrix of size LxD
+// @param w_ The (odd) number of values to average
+//
+// @return The filtered matrix, of same size as the input
+SEXP filterMA(SEXP M_, SEXP w_)
{
- int L = INTEGER(getAttrib(cwt_, R_DimSymbol))[0],
- D = INTEGER(getAttrib(cwt_, R_DimSymbol))[1];
- double *cwt = REAL(cwt_);
- SEXP R_fcwt;
- PROTECT(R_fcwt = allocMatrix(REALSXP, L, D));
- double* fcwt = REAL(R_fcwt);
-
- //TODO: coding style is terrible... no time for now.
- for (int col=0; col<D; col++)
+ int L = INTEGER(getAttrib(M_, R_DimSymbol))[0],
+ D = INTEGER(getAttrib(M_, R_DimSymbol))[1],
+ w = INTEGER_VALUE(w_),
+ half_w = (w-1)/2,
+ i,
+ nt; //number of terms in the current sum (max: w)
+ double *M = REAL(M_),
+ cs, //current sum (max: w terms)
+ left; //leftward term in the current moving sum
+
+ SEXP fM_; //the filtered result
+ PROTECT(fM_ = allocMatrix(REALSXP, L, D));
+ double* fM = REAL(fM_); //pointer to the encapsulated vector
+
+ // NOTE: unused loop parameter: shifting at the end of the loop is more efficient
+ for (int col=D-1; col>=0; col--)
{
- double v1 = cwt[0];
- double ma = v1 + cwt[1] + cwt[2];
- for (int i=1; i<L-2; i++)
+ // Initialization
+ nt = half_w + 1;
+ left = M[0];
+ cs = 0.;
+ for (i=half_w; i>=0; i--)
+ cs += M[i];
+
+ // Left border
+ for (i=0; i<half_w; i++)
{
- fcwt[i] = ma / 3.;
- ma = ma - v1 + cwt[i+2];
- v1 = cwt[i];
+ fM[i] = cs / nt; //(partial) moving average at current index i
+ cs += M[i+half_w+1];
+ nt++;
}
- fcwt[0] = cwt[0];
- fcwt[L-2] = ma / 3.;
- fcwt[L-1] = cwt[L-1];
- cwt = cwt + L;
- fcwt = fcwt + L;
+
+ // Middle: now nt == w, i == half_w
+ for (; i<L-half_w-1; i++)
+ {
+ fM[i] = cs / w; //moving average at current index i
+ cs = cs - left + M[i+half_w+1]; //remove oldest items, add next
+ left = M[i-half_w+1]; //update first value for next iteration
+ }
+
+ // (Last "fully averaged" index +) Right border
+ for (; i<L; i++)
+ {
+ fM[i] = cs / nt; //(partial) moving average at current index i
+ cs -= M[i-half_w];
+ nt--;
+ }
+
+ // Shift by L == increment column index by 1 (storage per column)
+ M = M + L;
+ fM = fM + L;
}
UNPROTECT(1);
- return R_fcwt;
+ return fM_;
}