// @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_;
}