X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=epclust%2Fsrc%2Ffilter.c;fp=epclust%2Fsrc%2Ffilter.c;h=97dbef2983fe2a4ea19a8bf30fcd16707f6e05c0;hb=40f12a2f66d06fd77183ea02b996f5c66f90761c;hp=0000000000000000000000000000000000000000;hpb=a52836b23adb4bfa6722642ec6426fb7b5f39650;p=epclust.git

diff --git a/epclust/src/filter.c b/epclust/src/filter.c
new file mode 100644
index 0000000..97dbef2
--- /dev/null
+++ b/epclust/src/filter.c
@@ -0,0 +1,70 @@
+#include <R.h>
+#include <Rdefines.h>
+
+// 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],
+		w = INTEGER_VALUE(w_),
+		half_w = (w-1)/2,
+		i,
+		nt; //number of terms in the current sum (max: w)
+	double *cwt = REAL(cwt_),
+		cs, //current sum (max: w terms)
+		left; //leftward term in the current moving sum
+
+	SEXP fcwt_; //the filtered result
+	PROTECT(fcwt_ = allocMatrix(REALSXP, L, D));
+	double* fcwt = REAL(fcwt_); //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--)
+	{
+		// Initialization
+		nt = half_w + 1;
+		left = cwt[0];
+		cs = 0.;
+		for (i=half_w; i>=0; i--)
+			cs += cwt[i];
+
+		// Left border
+		for (i=0; i<half_w; i++)
+		{
+			fcwt[i] = cs / nt; //(partial) moving average at current index i
+			cs += cwt[i+half_w+1];
+			nt++;
+		}
+
+		// Middle: now nt == w, i == half_w
+		for (; i<L-half_w-1; i++)
+		{
+			fcwt[i] = cs / w; //moving average at current index i
+			cs = ms - left + cwt[i+half_w+1]; //remove oldest items, add next
+			left = cwt[i-half_w+1]; //update first value for next iteration
+		}
+
+		// (Last "fully averaged" index +) Right border
+		for (; i<L; i++)
+		{
+			fcwt[i] = cs / nt; //(partial) moving average at current index i
+			cs -= cwt[i-half_w];
+			nt--;
+		}
+
+		// Shift by L == increment column index by 1 (storage per column)
+		cwt = cwt + L;
+		fcwt = fcwt + L;
+	}
+
+	UNPROTECT(1);
+	return fcwt_;
+}