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 }