#include <Origin.h>
#include <matrix.h>


void tfilter(int iMethod = 0)
{
	// Declare worksheet object from active page - quit if not valid worksheet
	Worksheet wks = Project.ActiveLayer();
	if(!wks)
	{
		printf("Active page is not a worksheet!\n");
		return;
	}
	
	// Declare x,y datasets from first two columns - quit if they are not valid
	Dataset dsX(wks, 0);
	Dataset dsY(wks, 1);	
	if(!dsX || !dsY)
	{
		printf("Could not define X, Y datasets!\n");
		return;
	}
	
	// Get size of data to be filtered
	int np = dsY.GetSize();
	
	// Declare vectors to be used with FFT call
	vector<complex> vecSignal(np), vecFFT(np);
	vecSignal=dsY;

	// Compute FFT - internally this calls the NAG library functions
	int nRet=FFT(vecSignal, vecFFT);
	if(nRet!=0)
	{
		printf("Error when computing FFT. Error code: %d\n", nRet);
		return;
	}
	
	// Get amplitude of FFT result
	vector vecFFTAmpl;
	vecFFT.GetAmplitude(vecFFTAmpl);
	
	// Create a temp wks and store FFT results in there
	Worksheet wksFFT;
	wksFFT.Create("Origin.OTW", CREATE_HIDDEN);
	Dataset dsFreq(wksFFT, 0);
	Dataset dsAmpl(wksFFT, 1);
	dsFreq.SetSize(np / 2);		
	// Compute freq column
	double dFreqStep = 1.0 / (dsX[1] - dsX[0]) / np;
	for(int ii = 0; ii < np / 2; ii++)
	{
		dsFreq[ii] = ii * dFreqStep;
	}
	// Fill amplitude
	dsAmpl = vecFFTAmpl;
	dsAmpl.SetSize(np / 2);

	// Create a graph and show user the FFT spectrum in log scale
	Curve crvFFT(dsFreq.GetName(), dsAmpl.GetName());
	GraphPage gpg;
	gpg.Create("Origin.OTP");
	GraphLayer gly = gpg.Layers();
	nRet = gly.AddPlot(crvFFT, IDM_PLOT_LINE);	
	Scale sc(gly.Y);                     // Get X scale object from graph layer
	sc.Type = LOG10_SPACE;              // Set axis scale as Log10
	gly.Rescale();

	// Ask user for threshold cutoff value
	string	strInput = InputBox("Enter cutoff value:");
	double dCutoff = atof(strInput);
	if(dCutoff <= 0)
	{
		printf("Invalid entry!\n");
		return;
	}

	// Remove all components above or below this value, based on method passed
	for(ii = 0; ii < np; ii++)
	{
		if(0 == iMethod)
		{
			// Drop all components below cutoff value
			if(vecFFTAmpl[ii] < dCutoff) vecFFT[ii] = 0.;
		}
		else
		{
			// Drop all components above cutoff value
			if(vecFFTAmpl[ii] > dCutoff) vecFFT[ii] = 0.;
		}
	}
	

	// Compute inverse FFT and extract real part
	nRet = IFFT(vecFFT, vecFFT);
	if(nRet!=0)
	{
		printf("Error when computing inverse FFT. Error code: %d\n", nRet);
		return;
	}

	// Add another column to worksheet next to data, and place result
	int iCol = wks.AddCol("Filt");
	Dataset dsResult(wks, iCol);
	vector vecReal(np);
	vecFFT.GetReal(vecReal);
	dsResult = vecReal;
	
	// Delete temp graph and wks
	wksFFT.Destroy();
	gpg.Destroy();
	
	// Create a new graph to show original data and filtered result
	Curve crvData(dsX.GetName(), dsY.GetName());
	Curve crvResult(dsX.GetName(), dsResult.GetName());
	GraphPage gpage;
	gpage.Create("Origin.OTP");
	GraphLayer glyr = gpage.Layers();
	// Add data to plot
	glyr.AddPlot(crvData, IDM_PLOT_LINE);
	// Add filter result and set color to red
	glyr.AddPlot(crvResult, IDM_PLOT_LINE);
	DataPlot dp2 = glyr.DataPlots(1);
	dp2.SetColor(1, TRUE);
	glyr.Rescale();	
}