/* SigLib Chirp z-Transform Example
The chirp z-transform algorithm is shown in the following diagram :
				_____
  vL		   |	 |	VL
  -------------| FFT |-------
			   |_____|		|
				_____		|      ______	
  xN	yN->yL |	 |	YL	|  GL |		 | gL->gM	   xM
  ----O--------| FFT |------O-----| IFFT |---------O------
	  |		   |_____|			  |______|		   |
	  |											   |
	  | AWNr, AWNi								   | wM

The input vector length is 150 samples, the output vector length is 100 samples
and the FFTs are 256 points, which must be >> (input + output lengths - 1).
The input data files have been sampled at 10 KHz.
In some applications a pre-window can assist the analysis e.g. Hanning.

This program individually initialises all the vectors for the sake of
demonstration. A neater way of coding this is to use the czt_init()
function.
*/

/* Include files */
#include <stdio.h>
#include <conio.h>
#include <siglib.h>
#include <nhl.h>
#include "GraphFunctions.h"

/* Define constants */
#define DEBUG			0						/* Set to 0 for no debug */
#define	SAMPLE_RATE		10000.0					/* 10 KHz */
#define	INPUT_LENGTH	((SLArrayIndex_t)150)	/* Input array length */
#define	OUTPUT_LENGTH	((SLArrayIndex_t)100)
#define	WINDOW_SIZE		INPUT_LENGTH
#define	FFT_SIZE		256L
#define	LOG_FFT_SIZE	((SLArrayIndex_t)8)

/* Declare global variables */
SLFixData_t		waveform;
SLData_t		Radius, Decay, StartFreq, EndFreq;
SLComplexPolar_s	ContourStart;
SLComplexRect_s	temp;

SLData_t		deltaomega, deltasigma;			/* Variables used to calculate W */
SLData_t		phinc, w1inc, w2inc;

SLComplexRect_s	A_1, W1, W_1, W12, W_12;		/* Complex contour coeffs */

SLData_t		*pInput, *pResults;
SLData_t		*pRealData, *pImagData, *pWindowCoeffs, *pFFTCoeffs;

SLData_t		*pAWNr, *pAWNi, *pvLr, *pvLi, *pWMr, *pWMi;
SLData_t		*pCZTRealWork, *pCZTImagWork;

void main(void);

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

	pInput = SUF_VectorArrayAllocate (INPUT_LENGTH);
	pRealData = SUF_VectorArrayAllocate (FFT_SIZE);
	pImagData = SUF_VectorArrayAllocate (FFT_SIZE);
	pWindowCoeffs = SUF_VectorArrayAllocate (WINDOW_SIZE);

	pAWNr = SUF_VectorArrayAllocate (INPUT_LENGTH);
	pAWNi = SUF_VectorArrayAllocate (INPUT_LENGTH);

	pvLr = SUF_VectorArrayAllocate (FFT_SIZE);
	pvLi = SUF_VectorArrayAllocate (FFT_SIZE);

	pCZTRealWork = SUF_VectorArrayAllocate (FFT_SIZE);
	pCZTImagWork = SUF_VectorArrayAllocate (FFT_SIZE);

	pWMr = SUF_VectorArrayAllocate (OUTPUT_LENGTH);
	pWMi = SUF_VectorArrayAllocate (OUTPUT_LENGTH);

	pResults = SUF_VectorArrayAllocate (OUTPUT_LENGTH);
	pFFTCoeffs = SUF_FftCoefficientAllocate (FFT_SIZE);


	if ((pInput == NULL) || (pRealData == NULL) || (pImagData == NULL) || (pWindowCoeffs == NULL) ||
		(pAWNr == NULL) || (pAWNi == NULL) || (pvLr == NULL) || (pvLi == NULL) || 
		(pWMr == NULL) || (pWMi == NULL) || (pResults == NULL) || (pFFTCoeffs == NULL))
	{
		printf ("\n\nMalloc failed\n\n");
		exit (0);
	}

	system ("cls");
	printf ("\n\n\n");

	printf ("Choose a waveform file\n\n");
	printf ("(Sample rate = 10 KHz)\n\n");
	printf ("AA.SIG ... (1)\n");
	printf ("EE.SIG ... (2)\n");
	printf ("ER.SIG ... (3)\n");
	printf ("I.SIG  ... (4)\n");
	printf ("OO.SIG ... (5)\n\n");
	printf (">");
	scanf ("%d", &waveform);

	switch (waveform)
	{
		case 1 :
			printf ("Opening AA.SIG\n\n");
			read_buffer (pInput, "aa.sig", INPUT_LENGTH);	/* Read data from disk */
			break;

		case 2 :
			printf ("Opening EE.SIG\n\n");
			read_buffer (pInput, "ee.sig", INPUT_LENGTH);	/* Read data from disk */
			break;

		case 3 :
			printf ("Opening ER.SIG\n\n");
			read_buffer (pInput, "er.sig", INPUT_LENGTH);	/* Read data from disk */
			break;

		case 4 :
			printf ("Opening I.SIG\n\n");
			read_buffer (pInput, "i.sig", INPUT_LENGTH);	/* Read data from disk */
			break;

		case 5 :
			printf ("Opening OO.SIG\n\n");
			read_buffer (pInput, "oo.sig", INPUT_LENGTH);	/* Read data from disk */
			break;

		default :
			printf ("\nInvalid result\n\n");
			exit (0);
			break;
	}

	printf ("\n\nEnter radius (r>0.9) >");
	scanf ("%lf", &Radius);

	printf ("\n\nEnter decay (Delta sigma / 2 * SIGLIB_PI) >");
	scanf ("%lf", &Decay);

	printf ("\n\nEnter start_freq >");
	scanf ("%lf", &StartFreq);

	printf ("\n\nEnter end_freq >");
	scanf ("%lf", &EndFreq);

	getchar ();						/* Clear the CR in the pipeline */

									/* Calculate A0 values */
	ContourStart = SCV_Polar (Radius, (StartFreq / SAMPLE_RATE) * SIGLIB_TWO_PI);
	temp = SCV_PolarToRectangular (ContourStart);
	A_1 = SCV_Inverse (temp);

									/* Calculate W0 values - W^(N^2/2) */
	deltaomega = (EndFreq - StartFreq) / ((SLData_t)OUTPUT_LENGTH - SIGLIB_ONE);
	deltasigma = Decay;

	phinc = SIGLIB_TWO * (-SIGLIB_PI) * deltaomega / SAMPLE_RATE;
	w1inc = SDS_Exp (SIGLIB_TWO * SIGLIB_PI * deltasigma / SAMPLE_RATE);
	w2inc = SDS_Sqrt (w1inc);


	W1.real = w1inc * SDS_Cos (phinc);									/* W */
	W1.imag = w1inc * sin (phinc);

	W12.real = w2inc * SDS_Cos (phinc / SIGLIB_TWO);					/* W^(1/2) */
	W12.imag = w2inc * SDS_Sin (phinc / SIGLIB_TWO);

	W_1.real = (SIGLIB_ONE / w1inc) * SDS_Cos (-phinc);					/* W^(-1) */
	W_1.imag = (SIGLIB_ONE / w1inc) * SDS_Sin (-phinc);

	W_12.real = (SIGLIB_ONE / w2inc) * SDS_Cos (-phinc / SIGLIB_TWO);	/* W^(-1/2) */
	W_12.imag = (SIGLIB_ONE / w2inc) * SDS_Sin (-phinc / SIGLIB_TWO);


#if DEBUG
	printf ("phinc = %lf\n", phinc);
	printf ("w1inc = %lf\n", w1inc);
	printf ("w2inc = %lf\n\n", w2inc);

	printf ("A_1.real = %lf\n", A_1.real);
	printf ("A_1.imag = %lf\n\n", A_1.imag);

	printf ("W1.real = %lf\n", W1.real);
	printf ("W1.imag = %lf\n", W1.imag);
	printf ("W12.real = %lf\n", W12.real);
	printf ("W12.imag = %lf\n\n", W12.imag);
	
	printf ("W_1.real = %lf\n", W_1.real);
	printf ("W_1.imag = %lf\n", W_1.imag);
	printf ("W_12.real = %lf\n", W_12.real);
	printf ("W_12.imag = %lf\n\n", W_12.imag);
	getchar ();
#endif

	h2DGraph =										/* Initialize graph */
		Create2DGraph ("Chirp z-Transform",			/* 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);
	}

													/* Initialise FFT */
	SIF_Fft (pFFTCoeffs,							/* Pointer to FFT coefficients */
			 SIGLIB_NULL_ARRAY_INDEX_PTR,			/* Pointer to bit reverse address table - NOT USED */
			 FFT_SIZE);								/* FFT Size */
													/* Gen. complex window coeffs */
	SIF_Awn (pAWNr,									/* Real coefficient pointer */
			 pAWNi,									/* Imaginary coefficient pointer */
			 A_1,									/* A ^ (-1) */
			 W1,									/* W */
			 W12,									/* W^(1/2) */
			 INPUT_LENGTH);							/* Array length */
													/* Gen. contour definition coeffs */
	SIF_Vl (pvLr,									/* Real coefficient pointer */
			pvLi,									/* Imaginary coefficient pointer */
			W_1,									/* W^(-1) */
			W_12,									/* W^(-1/2) */
			INPUT_LENGTH,							/* Input array length */
			OUTPUT_LENGTH,							/* Output array length */
			FFT_SIZE);								/* FFT array length */
													/* Gen. weighting coeffs */
	SIF_Wm (pWMr,									/* Real coefficient pointer */
			pWMi,									/* Imaginary coefficient pointer */
			W1,										/* W */
			W12,									/* W^(1/2) */
			OUTPUT_LENGTH);							/* Array length */


	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Source Signal",				/* Title of the dataset */
				    pInput,							/* Array of Double dataset */
				    INPUT_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_BLUE,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	printf ("\nSource Signal\nPlease hit <Carriage Return> to continue . . ."); getchar ();

													/* Ensure zero padded samples */
	SDA_Clear (pCZTRealWork,						/* Pointer to destination array */
			   FFT_SIZE);							/* Array length */
	SDA_Clear (pCZTImagWork,						/* Pointer to destination array */
			   FFT_SIZE);							/* Array length */

													/* Complex window = complex mpy with real data */
	SDA_ComplexWindow (pInput,						/* Pointer to real source array */
					   pInput,						/* Pointer to imaginary source array */
					   pCZTRealWork,				/* Pointer to real destination array */
					   pCZTImagWork,				/* Pointer to imaginary destination array */
					   pAWNr,						/* Real window array pointer */
					   pAWNi,						/* Imaginary window array pointer */
					   INPUT_LENGTH);				/* Window size */

													/* Frequency domain convolution */
	SDA_Cfft (pCZTRealWork,							/* Pointer to real array */
			  pCZTImagWork,							/* 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 */

	SDA_ComplexMultiply2 (pCZTRealWork,				/* Pointer to real source array 1 */
						  pCZTImagWork,				/* Pointer to imaginary source array 1 */
						  pvLr,						/* Pointer to real source array 2 */
						  pvLi,						/* Pointer to imaginary source array 2 */
						  pCZTRealWork,				/* Pointer to real destination array */
						  pCZTImagWork,				/* Pointer to imaginary destination array */
						  FFT_SIZE);				/* Array lengths */

													/* IFFT */
	SDA_Cifft (pCZTRealWork,						/* Pointer to real array */
			   pCZTImagWork,						/* 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 */

													/* Complex multiply */
	SDA_ComplexMultiply2 (pWMr,						/* Pointer to real source array 1 */
						  pWMi,						/* Pointer to imaginary source array 1 */
						  pCZTRealWork,				/* Pointer to real source array 2 */
						  pCZTImagWork,				/* Pointer to imaginary source array 2 */
						  pRealData,				/* Pointer to real destination array */
						  pImagData,				/* Pointer to imaginary destination array */
						  OUTPUT_LENGTH);			/* Array lengths */

													/* Scale chirp z-transform results */
	SDA_ComplexScalarDivide (pRealData,				/* Pointer to real source array */
							 pImagData,				/* Pointer to imaginary source array */
				  			 ((SLData_t)INPUT_LENGTH) * ((SLData_t)FFT_SIZE),	/* Divisor */
							 pRealData,				/* Pointer to real destination array */
							 pImagData,				/* Pointer to imaginary destination array */
							 OUTPUT_LENGTH);		/* Array lengths */


													/* Calculate real magnitude from complex */
	SDA_Magnitude (pRealData,						/* Pointer to real source array */
				   pImagData,						/* Pointer to imaginary source array */
				   pResults,						/* Pointer to magnitude destination array */
				   OUTPUT_LENGTH);					/* Array length */
													/* Scale results so peaks equal 100.0 */
	SDA_Scale (pResults,							/* Pointer to source array */
			   pResults,							/* Pointer to destination array */
			   ((SLData_t)100.0),					/* Peak level */
			   OUTPUT_LENGTH);						/* Array length */
													/* dB scaling */
	SDA_20Log10 (pResults,							/* Pointer to source array */
				 pResults,							/* Pointer to destination array */
				 OUTPUT_LENGTH);					/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "chirp z-Transform",			/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    OUTPUT_LENGTH,					/* Number of data points */
					SV_GRAPH_LINE,					/* Graph type */
				    SV_BLUE,						/* Colour */
					SV_HIDE_MARKERS,				/* Marker enable / disable */
					SV_GRAPH_NEW);					/* New graph */
	printf ("\nchirp z-Transform\n");

	SUF_MemoryFree (pInput);						/* Free memory */
	SUF_MemoryFree (pRealData);
	SUF_MemoryFree (pImagData);
	SUF_MemoryFree (pWindowCoeffs);
	SUF_MemoryFree (pAWNr);
	SUF_MemoryFree (pAWNi);
	SUF_MemoryFree (pvLr);
	SUF_MemoryFree (pvLi);
	SUF_MemoryFree (pCZTRealWork);
	SUF_MemoryFree (pCZTImagWork);
	SUF_MemoryFree (pWMr);
	SUF_MemoryFree (pWMi);
	SUF_MemoryFree (pResults);
	SUF_MemoryFree (pFFTCoeffs);
}