#include <Rcpp.h>

using namespace Rcpp;

//' filter
//'
//' Filter time-series
//'
//' @param cwt Continuous wavelets transform
//'
//' @return The filtered CWT
// [[Rcpp::export]]
NumericMatrix epclustFilter(NumericMatrix cwt)
{
	int L = cwt.nrow(),
		D = cwt.ncol();
	NumericMatrix fcwt(L, D); //fill with 0... TODO: back to SEXP C-style?
	double *cwt_c = cwt.begin(),
		*fcwt_c = fcwt.begin();

	//TODO: coding style is terrible... no time for now.
	for (int col=0; col<D; col++)
	{
		double v1 = cwt_c[0];
		double ma = v1 + cwt[1] + cwt_c[2];
		for (int i=1; i<L-2; i++)
		{
			fcwt_c[i] = ma / 3.;
			ma = ma - v1 + cwt_c[i+2];
			v1 = cwt_c[i];
		}
		fcwt_c[0] = cwt_c[0];
		fcwt_c[L-2] = ma / 3.;
		fcwt_c[L-1] = cwt_c[L-1];
		cwt_c = cwt_c + L;
		fcwt_c = fcwt_c + L;
	}

	return fcwt;
}