/* SigLib Frequency Domain Overlapped FIR Filter Example
This program filters a signal in the time domain and then in
the frequency domain using the overlap-add and overlap-save
techniques */

#include <stdio.h>
#include <siglib.h>
#include "GraphFunctions.h"

/* Define constants */
#define	FILTER_LENGTH			((SLArrayIndex_t)64)
#define	SAMPLE_LENGTH			((SLArrayIndex_t)512)
#define	INVERSE_SAMPLE_LENGTH	(SIGLIB_ONE / ((SLData_t)SAMPLE_LENGTH))

#define	FFT_LENGTH				((SLArrayIndex_t)128)
#define	LOG_FFT_LENGTH			7L
#define TIME_DOMAIN_DATA_LENGTH	((SLArrayIndex_t)64)		/* Each iteration processes this many data samples */

			/* Initialise filter coefficients */
SLData_t	TimeDomainCoeffs[FILTER_LENGTH] = {
	-3.783E-3,  2.803E-3,  2.648E-3,  2.891E-3,  3.397E-3,  4.094E-3, 
 	 4.942E-3,  5.917E-3,  7.005E-3,  8.196E-3,  9.479E-3,  1.084E-2, 
 	 1.229E-2,  1.379E-2,  1.535E-2,  1.695E-2,  1.857E-2,  2.021E-2, 
 	 2.183E-2,  2.344E-2,  2.500E-2,  2.651E-2,  2.795E-2,  2.929E-2, 
 	 3.053E-2,  3.165E-2,  3.264E-2,  3.348E-2,  3.417E-2,  3.469E-2, 
 	 3.504E-2,  3.522E-2,  3.522E-2,  3.504E-2,  3.469E-2,  3.417E-2, 
 	 3.348E-2,  3.264E-2,  3.165E-2,  3.053E-2,  2.929E-2,  2.795E-2, 
 	 2.651E-2,  2.500E-2,  2.344E-2,  2.183E-2,  2.021E-2,  1.857E-2, 
 	 1.695E-2,  1.535E-2,  1.379E-2,  1.229E-2,  1.084E-2,  9.479E-3, 
 	 8.196E-3,  7.005E-3,  5.917E-3,  4.942E-3,  4.094E-3,  3.397E-3, 
 	 2.891E-3,  2.648E-3,  2.803E-3, -3.783E-3
	};


SLData_t	pFilterState[FILTER_LENGTH];
SLArrayIndex_t	FilterIndex;
SLData_t	*pSrc1, *pSrc2, *pSrc3;
SLData_t	*RealFreqDomainCoeffs, *ImagFreqDomainCoeffs;
SLData_t	*OverlapArrayPtr;
SLData_t	*TempArrayPtr;
SLData_t	SinePhase, CosinePhase;
SLData_t	*pFFTCoeffs;
SLData_t	InverseFFTSize;

void main(void)
{
	GraphObject *h2DGraph;							/* Declare graph object */

	SLArrayIndex_t	i;
	SLData_t		MSE;

													/* Allocate memory */
	pSrc1 = SUF_VectorArrayAllocate (SAMPLE_LENGTH + FFT_LENGTH);	/* Allow for overlap - overlap and save */
	pSrc2 = SUF_VectorArrayAllocate (SAMPLE_LENGTH);
	pSrc3 = SUF_VectorArrayAllocate (SAMPLE_LENGTH + FFT_LENGTH);	/* Allow for overlap - overlap and add */

	RealFreqDomainCoeffs = SUF_VectorArrayAllocate (FFT_LENGTH);	/* Allocate working arrays */
	ImagFreqDomainCoeffs = SUF_VectorArrayAllocate (FFT_LENGTH);
	OverlapArrayPtr = SUF_VectorArrayAllocate (FFT_LENGTH);			/* Allocate overlap storage */
	TempArrayPtr = SUF_VectorArrayAllocate (FFT_LENGTH);
	pFFTCoeffs = SUF_FftCoefficientAllocate (FFT_LENGTH);

													/* Clear array to ensure zero padded - overlap and save */
	SDA_Clear (pSrc1,								/* Pointer to destination array */
			   SAMPLE_LENGTH + FFT_LENGTH);			/* Array length */

	h2DGraph =										/* Initialize graph */
		Create2DGraph ("Frequency Domain Filtering",	/* Graph title */
					   "Time",						/* X-Axis label */
					   "Magnitude",					/* Y-Axis label */
					   SV_AUTO_SCALE,				/* Scaling mode */
					   SV_SIGNED,					/* Sign mode */
					   SV_GRAPH_LINE,				/* Graph type */
					   "localhost");				/* Graph server */
	if (h2DGraph == NULL)							/* Graph creation failed - e.g is server running ? */
	{
		printf ("\nGraph creation failure. Please check that the server is running\n");
		exit (1);
	}


	SinePhase = SIGLIB_ZERO;
	CosinePhase = SIGLIB_ZERO;

										/* Generate a noisy sinewave */
	SDA_SignalGenerate (pSrc1,						/* Pointer to destination array */
						SIGLIB_SINE_WAVE,			/* Signal type - Sine wave */
						SIGLIB_HALF,				/* Signal peak level */
						SIGLIB_FILL,				/* Fill (overwrite) or add to existing array contents */
						((SLData_t)0.01),			/* Signal frequency */
						SIGLIB_ZERO,				/* D.C. Offset */
						SIGLIB_ZERO,				/* Unused */
						SIGLIB_ZERO,				/* Signal end value - Unused */
						&SinePhase,					/* Signal phase - maintained across array boundaries */
						SIGLIB_NULL_DATA_PTR,		/* Unused */
						SAMPLE_LENGTH);				/* Output array length */

	SDA_SignalGenerate (pSrc1,						/* Pointer to destination array */
						SIGLIB_COS_WAVE,			/* Signal type - Cosine wave */
						SIGLIB_HALF,				/* Signal peak level */
						SIGLIB_ADD,					/* Fill (overwrite) or add to existing array contents */
						((SLData_t)0.01),			/* Signal frequency */
						SIGLIB_ZERO,				/* D.C. Offset */
						SIGLIB_ZERO,				/* Unused */
						SIGLIB_ZERO,				/* Signal end value - Unused */
						&CosinePhase,				/* Signal phase - maintained across array boundaries */
						SIGLIB_NULL_DATA_PTR,		/* Unused */
						SAMPLE_LENGTH);				/* Output array length */

	SDA_SignalGenerate (pSrc1,						/* Pointer to destination array */
						SIGLIB_WHITE_NOISE,			/* Signal type - random white noise */
						((SLData_t)0.2),			/* Signal peak level */
						SIGLIB_ADD,					/* Fill (overwrite) or add to existing array contents */
						SIGLIB_ZERO,				/* Signal frequency - Unused */
						SIGLIB_ZERO,				/* D.C. Offset */
						SIGLIB_ZERO,				/* Unused */
						SIGLIB_ZERO,				/* Signal end value - Unused */
						SIGLIB_NULL_DATA_PTR,		/* Unused */
						SIGLIB_NULL_DATA_PTR,		/* Unused */
						SAMPLE_LENGTH);				/* Output array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Unfiltered Signal",			/* Title of the dataset */
				    pSrc1,							/* Array of Double dataset */
				    SAMPLE_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_BLUE,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */

													/* Initialise FIR filter */
	SIF_Fir (pFilterState,							/* Pointer to filter state array */
			 &FilterIndex,							/* Pointer to filter index register */
			 FILTER_LENGTH);						/* Filter length */

													/* Apply fir filter */
	SDA_Fir (pSrc1,									/* Input array to be filtered */
			 pSrc2,									/* Filtered output array */
			 pFilterState,							/* Pointer to filter state array */
			 TimeDomainCoeffs,						/* Pointer to filter coefficients */
			 &FilterIndex,							/* Pointer to filter index register */
			 FILTER_LENGTH,							/* Filter length */
			 SAMPLE_LENGTH);						/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Time Domain Filtered Signal",	/* Title of the dataset */
				    pSrc2,							/* Array of Double dataset */
				    SAMPLE_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_RED,							/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_ADD);					/* New graph */
	printf ("\nTime Domain Filtered Signal\nPlease hit <Carriage Return> to continue . . ."); getchar ();


	SIF_FirOverlapAdd (TimeDomainCoeffs,			/* Time Domain Coefficients pointer */
					   RealFreqDomainCoeffs,		/* Ptr to Real Freq Domain Coefficients pointer */
					   ImagFreqDomainCoeffs,		/* Ptr to Imaginary Freq Domain Coefficients pointer */
					   OverlapArrayPtr,				/* Ptr to overlap array pointer */
					   pFFTCoeffs,					/* Pointer to FFT coeffs */
					   SIGLIB_NULL_ARRAY_INDEX_PTR,	/* Pointer to FFT bit reverse address table */
					   &InverseFFTSize,				/* Pointer to inverse FFT Length */
					   FFT_LENGTH,					/* FFT Length */
					   LOG_FFT_LENGTH,				/* Log10 FFT Length */
					   FILTER_LENGTH);				/* Filter length */

								/* Apply frequency domain (overlap and add) FIR filter */
	for (i = 0; i < SAMPLE_LENGTH; i += TIME_DOMAIN_DATA_LENGTH)
	{
		SDA_FirOverlapAdd (pSrc1+i,						/* Source data pointer */
						   pSrc3+i,						/* Destination data pointer */
						   RealFreqDomainCoeffs,		/* Real Freq Domain Coefficients pointer */
						   ImagFreqDomainCoeffs,		/* Imaginary Freq Domain Coefficients pointer */
						   OverlapArrayPtr,				/* Overlap array pointer */
						   TempArrayPtr,				/* Temporary array pointer */
						   pFFTCoeffs,					/* Pointer to FFT coeffs */
						   SIGLIB_NULL_ARRAY_INDEX_PTR,	/* Pointer to FFT bit reverse address table */
						   InverseFFTSize,				/* Inverse FFT Length */
						   FFT_LENGTH,					/* FFT Length */
						   LOG_FFT_LENGTH,				/* Log 10 FFT Length */
						   FILTER_LENGTH,				/* Filter length */
						   TIME_DOMAIN_DATA_LENGTH);	/* Array length */
	}
														/* Calculate mean square error */
	MSE =
		SDA_MeanSquareError (pSrc2,						/* Pointer to source array 1 */
							 pSrc3,						/* Pointer to source array 2 */
							 INVERSE_SAMPLE_LENGTH,		/* Inverse of the array length */
							 SAMPLE_LENGTH);			/* Array length */

	printf ("\nFrequency Domain (Overlap And Add) Filtered Signal MSE = %e\n", MSE);
	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Frequency Domain (Overlap And Add) Filtered Signal",	/* Title of the dataset */
				    pSrc3,							/* Array of Double dataset */
				    SAMPLE_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_MAGENTA,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_ADD);					/* New graph */
	printf ("\nFrequency Domain (Overlap And Add) Filtered Signal\nPlease hit <Carriage Return> to continue . . ."); getchar ();


	SIF_FirOverlapSave (TimeDomainCoeffs,				/* Time Domain Coefficients pointer */
						RealFreqDomainCoeffs,			/* Ptr to Real Freq Domain Coefficients pointer */
						ImagFreqDomainCoeffs,			/* Ptr to Imaginary Freq Domain Coefficients pointer */
						OverlapArrayPtr,				/* Ptr to overlap array pointer */
						pFFTCoeffs,						/* Pointer to FFT coeffs */
						SIGLIB_NULL_ARRAY_INDEX_PTR,	/* Pointer to FFT bit reverse address table */
						&InverseFFTSize,				/* Pointer to inverse FFT Length */
						FFT_LENGTH,						/* FFT Length */
						LOG_FFT_LENGTH,					/* Log10 FFT Length */
						FILTER_LENGTH);					/* Filter length */

								/* Apply frequency domain (overlap and add) FIR filter */
	for (i = 0; i < SAMPLE_LENGTH; i += TIME_DOMAIN_DATA_LENGTH)
	{
		SDA_FirOverlapSave (pSrc1+i,						/* Source data pointer */
							pSrc3+i,						/* Destination data pointer */
							RealFreqDomainCoeffs,			/* Real Freq Domain Coefficients pointer */
							ImagFreqDomainCoeffs,			/* Imaginary Freq Domain Coefficients pointer */
							OverlapArrayPtr,				/* Overlap array pointer */
							TempArrayPtr,					/* Temporary array pointer */
							pFFTCoeffs,						/* Pointer to FFT coeffs */
							SIGLIB_NULL_ARRAY_INDEX_PTR,	/* Pointer to FFT bit reverse address table */
							InverseFFTSize,					/* Inverse FFT Length */
							FFT_LENGTH,						/* FFT Length */
							LOG_FFT_LENGTH,					/* Log 10 FFT Length */
							FILTER_LENGTH,					/* Filter length */
							TIME_DOMAIN_DATA_LENGTH);		/* Array length */
	}
															/* Calculate mean square error */
	MSE =
		SDA_MeanSquareError (pSrc2,							/* Pointer to source array 1 */
							 pSrc3,							/* Pointer to source array 2 */
							 INVERSE_SAMPLE_LENGTH,			/* Inverse of the array length */
							 SAMPLE_LENGTH);				/* Array length */

	printf ("Frequency Domain (Overlap And Save) Filtered Signal MSE = %e\n", MSE);
	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Frequency Domain (Overlap And Save) Filtered Signal",	/* Title of the dataset */
				    pSrc3,							/* Array of Double dataset */
				    SAMPLE_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_GREEN,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_ADD);					/* New graph */

	SUF_MemoryFree (pSrc1);							/* Free memory */
	SUF_MemoryFree (pSrc2);
	SUF_MemoryFree (pSrc3);

}





