/* SigLib - FIR Filter Design By The Window Method

   +------+        +-------+		+--------+		  +---------+
   | Freq |_______\| 128Pt |_______\| Array  |_______\| Hamming |
   | Resp |       /| IFFT  |	   /| Rotate |		 /| Window  |
   +------+        +-------+		+--------+		  +---------+
                                                           |
           +-------+       Signal   +--------+  Coeffs     |
           | White |_______________\| FIR    |/____________|
           | Noise |               /| Filter |\
           +-------+                +--------+
       __________________________________|
       |
      \|/
   +-------+		+--------+		  +--------+        +----------+
   | 512Pt |_______\| Sq Mag |_______\| dB and |_______\| Display  |
   |  FFT  |	   /|  Sum   |		 /| Scale  |       /| Spectrum |
   +-------+		+--------+		  +--------+        +----------+
*/

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

/* Define constants */
#define	RAISED_COS			0				/* Set to '1' for raised cosine filter or '0' for arbitrary LPF */
#define FILTER_NOISE		0				/* Set to '1' to filter noise or '0' impulse response */
#define	FILTER_LENGTH		((SLArrayIndex_t)128)		/* Use a power of 2 size for FFT reasons */
#define	LOG_FILTER_LENGTH	((SLArrayIndex_t)7)
#define	FFT_SIZE			((SLArrayIndex_t)512)
#define	LOG_FFT_SIZE		((SLArrayIndex_t)9)

/* Declare arrays and variables */
SLArrayIndex_t	FilterIndex;
SLData_t	SrcData;

SLData_t	*pFilterCoeffs;							/* Filter coefficients */
SLData_t	*pFilterState;							/* Filter state array */

SLData_t	*pFiltered, *pImagData, *pWindowCoeffs, *pResults, *pFFTCoeffs;

void	main(void);

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

	SLArrayIndex_t		i;

#if FILTER_NOISE
	time_t		ltime;
	time (&ltime);
	srand ((unsigned int) ltime);					/* Randomise the seed */
#endif

	pFilterCoeffs = SUF_VectorArrayAllocate (FILTER_LENGTH);
	pFilterState = SUF_VectorArrayAllocate (FILTER_LENGTH);

	pFiltered = SUF_VectorArrayAllocate (FFT_SIZE);
	pImagData = SUF_VectorArrayAllocate (FFT_SIZE);
	pResults = SUF_VectorArrayAllocate (FFT_SIZE);
	pFFTCoeffs = SUF_FftCoefficientAllocate (FFT_SIZE);

	h2DGraph =										/* Initialize graph */
		Create2DGraph ("FIR Filter Design By The Window Method",	/* Graph title */
					   "Time / Frequency",			/* 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);
	}

													/* Init. FFT for generating filter coeffs */
	SIF_Fft (pFFTCoeffs,							/* Pointer to FFT coefficients */
			 SIGLIB_NULL_ARRAY_INDEX_PTR,			/* Pointer to bit reverse address table - NOT USED */
			 FILTER_LENGTH);						/* FFT Size */
													/* Clear imaginary array for real inverse FFT */
	SDA_Clear (pImagData,							/* Pointer to destination array */
			   FILTER_LENGTH);						/* Array length */

													/* Generate filter frequency response */
#if RAISED_COS										/* Raised cosine low pass filter */
	for (i = 0; i < FILTER_LENGTH; i++)
		*(pFilterCoeffs+i) = SIGLIB_HALF* (SIGLIB_ONE + SDS_Cos (SIGLIB_TWO * SIGLIB_PI * ((SLData_t)(i)) / ((SLData_t)(FILTER_LENGTH))));

#else												/* Even-real low pass filter */
	SDA_Clear (pFilterCoeffs,						/* Pointer to destination array */
			   FILTER_LENGTH);						/* Array length */

	*pFilterCoeffs = SIGLIB_ONE;					/* Central frequency domain coefficient */

	for (i = 0; i < FILTER_LENGTH/4; i++)
		*(pFilterCoeffs+i+1) = SIGLIB_ONE;

	for (i = 0; i < FILTER_LENGTH/4; i++)
		*(pFilterCoeffs+FILTER_LENGTH-i-1) = SIGLIB_ONE;
#endif

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Desired filter response",		/* Title of the dataset */
				    pFilterCoeffs,					/* Array of Double dataset */
				    FILTER_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_MAGENTA,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	printf ("\nDesired filter response\nPlease hit <Carriage Return> to continue . . ."); getchar ();

													/* Perform IFFT */
	SDA_Cifft (pFilterCoeffs,						/* Pointer to real array */
			   pImagData,							/* Pointer to imaginary array */
			   pFFTCoeffs,							/* Pointer to FFT coefficients */
			   SIGLIB_NULL_ARRAY_INDEX_PTR,			/* Pointer to bit reverse address table - NOT USED */
			   FILTER_LENGTH,						/* FFT size */
			   LOG_FILTER_LENGTH); 					/* log2 FFT size */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Real IFFT'd coefficients",		/* Title of the dataset */
				    pFilterCoeffs,					/* Array of Double dataset */
				    FILTER_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_MAGENTA,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Imag IFFT'd coefficients",		/* Title of the dataset */
				    pFilterCoeffs,					/* Array of Double dataset */
				    FILTER_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_BLUE,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_ADD);					/* New graph */
	printf ("\nIFFT'd coefficients\nPlease hit <Carriage Return> to continue . . ."); getchar ();

	pWindowCoeffs = SUF_VectorArrayAllocate (FILTER_LENGTH);	/* Initialise windowing function */
	SIF_Window (pWindowCoeffs,						/* Pointer to window oefficient */
				SIGLIB_HANNING,						/* Window type */
				SIGLIB_ZERO,						/* Window coefficient */
				FILTER_LENGTH);						/* Window length */

													/* Shift the FFT coefficients */
	SDA_FftShift (pFilterCoeffs,					/* Pointer to source array */
				  pFilterCoeffs,					/* Pointer to destination array */
				  FILTER_LENGTH);					/* Array length */

										 			/* Window the filter coefficients */
	SDA_Window (pFilterCoeffs,						/* Pointer to source array */
				pFilterCoeffs,						/* Pointer to destination array */
				pWindowCoeffs,						/* Pointer to window oefficients */
				FILTER_LENGTH);						/* Window length */

													/* Scale coefficients to 1.0 */
	SDA_Scale (pFilterCoeffs,						/* Pointer to source array */
			   pFilterCoeffs,						/* Pointer to destination array */
			   SIGLIB_ONE,							/* Peak level */
			   FILTER_LENGTH);						/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Shifted, window'd and scaled coefficients",	/* Title of the dataset */
				    pFilterCoeffs,					/* Array of Double dataset */
				    FILTER_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_MAGENTA,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	printf ("\nShifted, window'd and scaled coefficients\nPlease hit <Carriage Return> to continue . . ."); getchar ();

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

	for (i = 0; i < FFT_SIZE; i++)					/* Filter data */
	{
#if FILTER_NOISE
				/* Generate a noisy signal to filter */
		SrcData = (((SLData_t)rand()) - ((SLData_t)16383.0)) / ((SLData_t)16384.0);
#else
		if (i == 0)
			SrcData = SIGLIB_ONE;
		else
			SrcData = SIGLIB_ZERO;
#endif

				/* Apply fir filter and store filtered data */
		*(pFiltered+i) =
			SDS_Fir (SrcData,						/* Input data sample to be filtered */
					 pFilterState,					/* Pointer to filter state array */
					 pFilterCoeffs,					/* Pointer to filter coefficients */
					 &FilterIndex,					/* Pointer to filter index register */
					 FILTER_LENGTH);				/* Filter length */
	}

													/* Init. FFT for calculating filter response */
	SIF_Fft (pFFTCoeffs,							/* Pointer to FFT coefficients */
			 SIGLIB_NULL_ARRAY_INDEX_PTR,			/* Pointer to bit reverse address table - NOT USED */
			 FFT_SIZE);								/* FFT Size */

													/* Perform real FFT */
	SDA_Rfft (pFiltered,							/* Pointer to real array */
			  pImagData,							/* Pointer to imaginary array */
			  pFFTCoeffs,							/* Pointer to FFT coefficients */
			  SIGLIB_NULL_ARRAY_INDEX_PTR,			/* Pointer to bit reverse address table - NOT USED */
			  FFT_SIZE,								/* FFT size */
			  LOG_FFT_SIZE);						/* log2 FFT size */

													/* Calculate real power from complex */
	SDA_LogMagnitude (pFiltered,					/* Pointer to real source array */
					  pImagData,					/* Pointer to imaginary source array */
					  pResults,						/* Pointer to log magnitude destination array */
					  FFT_SIZE);					/* Array length */

#if FILTER_NOISE
	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Spectrum of filtered noise signal",	/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    FFT_SIZE,						/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_MAGENTA,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	printf ("\nSpectrum of filtered noise signal\n);
#else
	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Spectrum of impulse response",	/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    FFT_SIZE,						/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_MAGENTA,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	printf ("\nSpectrum of impulse response\n");
#endif


	SUF_MemoryFree (pFiltered);						/* Free memory */
	SUF_MemoryFree (pImagData);
	SUF_MemoryFree (pResults);
	SUF_MemoryFree (pFilterCoeffs);
	SUF_MemoryFree (pFilterState);
	SUF_MemoryFree (pWindowCoeffs);
	SUF_MemoryFree (pFFTCoeffs);
}


