drop enercast submodule; drop Rcpp requirement; fix doc, complete code, fix fix fix
[epclust.git] / epclust / src / filter.c
index 97dbef2..8ece9c3 100644 (file)
 // @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],
+       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 *cwt = REAL(cwt_),
+       double *M = REAL(M_),
                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
+       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--)
        {
                // Initialization
                nt = half_w + 1;
-               left = cwt[0];
+               left = M[0];
                cs = 0.;
                for (i=half_w; i>=0; i--)
-                       cs += cwt[i];
+                       cs += M[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];
+                       fM[i] = cs / nt; //(partial) moving average at current index i
+                       cs += M[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
+                       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++)
                {
-                       fcwt[i] = cs / nt; //(partial) moving average at current index i
-                       cs -= cwt[i-half_w];
+                       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)
-               cwt = cwt + L;
-               fcwt = fcwt + L;
+               M = M + L;
+               fM = fM + L;
        }
 
        UNPROTECT(1);
-       return fcwt_;
+       return fM_;
 }