/* Bilinear Transform IIR Filter Design Example.
Generates a low pass filter and transforms the frequency */

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

/* Define constants */
#define	PRINT_RESULTS			0				/* Set to 1 to print results as well as plot them */

							/* Select which transforms to perform */
#define	TRANSFORM_L_L			1				/* Low pass to low pass transformation */
#define	TRANSFORM_L_H			1				/* Low pass to high pass transformation */
#define	TRANSFORM_L_P			1				/* Low pass to band pass transformation */
#define	TRANSFORM_L_S			1				/* Low pass to band stop transformation */

#define	SAMPLE_RATE				SIGLIB_ONE			/* Normalized to 1.0 Hz for convenience */
#define PREWARP_MATCH_FREQUENCY	((SLData_t)0.2)		/* Low pass filter cut-off frequency */
#define	NUMBER_OF_ZEROS			5L
#define	NUMBER_OF_POLES			5L
#define	NUMBER_OF_FILTER_COEFFS	(SIGLIB_IIR_COEFFS_PER_BIQUAD*NUMBER_OF_ZEROS)

#define FILTER_STAGES			NUMBER_OF_ZEROS		/* Number of 2nd-order filter stages */

#define	LPF_TO_LPF_NEW_FC		((SLData_t)0.3)
#define	LPF_TO_BPF_NEW_FC1		((SLData_t)0.1)
#define	LPF_TO_BPF_NEW_FC2		((SLData_t)0.3)
#define	LPF_TO_BSF_NEW_FC1		((SLData_t)0.1)
#define	LPF_TO_BSF_NEW_FC2		((SLData_t)0.3)
#define	LPF_TO_HPF_NEW_FC		((SLData_t)0.3)

#define	IMPULSE_RESPONSE_LENGTH	((SLArrayIndex_t)1024)

#define	FFT_SIZE				IMPULSE_RESPONSE_LENGTH
#define	LOG_FFT_SIZE			((SLArrayIndex_t)10)

#define	PLOT_LENGTH				(IMPULSE_RESPONSE_LENGTH/2L)

/* Declare global variables */
SLComplexRect_s SPlaneZeros [NUMBER_OF_ZEROS];
SLComplexRect_s SPlanePoles [NUMBER_OF_POLES];
SLComplexRect_s ZPlaneZeros [NUMBER_OF_ZEROS];
SLComplexRect_s ZPlanePoles [NUMBER_OF_POLES];
SLComplexRect_s TransformedZPlaneZeros [2*NUMBER_OF_ZEROS];		/* Need twice number of poles and zeros for BPFs and BSFs */
SLComplexRect_s TransformedZPlanePoles [2*NUMBER_OF_POLES];

SLData_t	*pFilterState, *pIIRCoeffs;
SLData_t	*pSrc, *pRealData, *pImagData, *pResults, *pFFTCoeffs;

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

													/* Allocate memory */
												/* Need twice number of poles and zeros for BPFs and BSFs */
	pIIRCoeffs = SUF_IirCoefficientAllocate (2 * FILTER_STAGES);
	pFilterState = SUF_IirStateArrayAllocate (2 * FILTER_STAGES);
	pRealData = SUF_VectorArrayAllocate (FFT_SIZE);
	pImagData = SUF_VectorArrayAllocate (FFT_SIZE);
	pResults = SUF_VectorArrayAllocate (PLOT_LENGTH);
	pSrc = SUF_VectorArrayAllocate (IMPULSE_RESPONSE_LENGTH);
	pFFTCoeffs = SUF_FftCoefficientAllocate (FFT_SIZE);

	if ((pIIRCoeffs == NULL) || (pRealData == NULL) || (pImagData == NULL) || (pResults == NULL) ||
		(pSrc == NULL) || (pFFTCoeffs == NULL))
	{
		printf ("\n\nMemory allocation failed\n\n");
		exit (0);
	}
								/* Initialise s-plane zero and pole positions */
								/* Create rectangular numbers from real and imaginary components */
	SPlaneZeros[0] = SCV_Rectangular (((SLData_t)0.0), ((SLData_t)3.0));
	SPlaneZeros[1] = SCV_Rectangular (((SLData_t)0.0), ((SLData_t)3.0));
	SPlaneZeros[2] = SCV_Rectangular (((SLData_t)0.0), ((SLData_t)3.0));
	SPlaneZeros[3] = SCV_Rectangular (((SLData_t)0.0), ((SLData_t)3.0));
	SPlaneZeros[4] = SCV_Rectangular (((SLData_t)0.0), ((SLData_t)3.0));

	SPlanePoles[0] = SCV_Rectangular (((SLData_t)-0.1), ((SLData_t)0.9));
	SPlanePoles[1] = SCV_Rectangular (((SLData_t)-0.3), ((SLData_t)0.8));
	SPlanePoles[2] = SCV_Rectangular (((SLData_t)-0.6), ((SLData_t)0.7));
	SPlanePoles[3] = SCV_Rectangular (((SLData_t)-0.8), ((SLData_t)0.6));
	SPlanePoles[4] = SCV_Rectangular (((SLData_t)-0.9), ((SLData_t)0.3));

#if	PRINT_RESULTS
	printf ("\nComplex s-plane zeros\n");			/* Print s-plane poles and zeros */
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
		print_rectangular (SPlaneZeros[i]);

	printf ("\nComplex s-plane poles\n");
	for (i = 0; i < NUMBER_OF_POLES; i++)
		print_rectangular (SPlanePoles[i]);
#endif

													/* Convert z-plane poles and zeros to s-plane */
													/* using the bilinear transform */
	SDA_BilinearTransform (SPlaneZeros,				/* S-plane zeros */
						   SPlanePoles,				/* S-plane poles */
						   ZPlaneZeros,				/* Z-plane zeros */
						   ZPlanePoles,				/* Z-plane poles */
						   SAMPLE_RATE,				/* Sample rate */
						   PREWARP_MATCH_FREQUENCY,	/* Pre-warp frequency */
						   SIGLIB_ON,				/* Pre-warp switch */
						   NUMBER_OF_ZEROS,			/* Number of zeros */
						   NUMBER_OF_POLES);		/* Number of poles */

#if	PRINT_RESULTS
	printf ("\nComplex z-plane zeros\n");			/* Print z-plane poles and zeros */
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
		print_rectangular (ZPlaneZeros[i]);

	printf ("\nComplex z-plane poles\n");
	for (i = 0; i < NUMBER_OF_POLES; i++)
		print_rectangular (ZPlanePoles[i]);
	printf ("Hit any key to see pole-zero plot . . ."); getchar();
#endif

	hPZGraph =										/* Initialize graph */
		CreatePoleZeroGraph("Pole-Zero Plot",		/* Graph title */
							"Real",					/* X-Axis label */
							"Imaginary",			/* Y-Axis label */
							SV_COMPLEX_ZERO,		/* PZ Data type */
							"localhost");			/* Graph server */
	if (hPZGraph == NULL)							/* Graph creation failed - e.g is server running ? */
	{
		printf ("\nGraph creation failure. Please check that the server is running\n");
		exit (1);
	}

	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Poles",							/* Title of dataset (for use on legends) */
				  (Complex_s *)ZPlanePoles,			/* Array of complex values */
				  NUMBER_OF_POLES,					/* Items of data */
				  SV_BLUE,							/* Colour */
				  SV_GRAPH_NEW,						/* Add / new mode */
				  SV_CONJUGATE_POLE);				/* Pole-zero mode */
	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Zeros",							/* Title of dataset (for use on legends) */
				  (Complex_s *)ZPlaneZeros,			/* Array of complex values */
				  NUMBER_OF_ZEROS,					/* Items of data */
				  SV_RED,							/* Colour */
				  SV_GRAPH_ADD,						/* Add / new mode */
				  SV_CONJUGATE_ZERO);				/* Pole-zero mode */
	printf ("Pole-zero Plot\nHit any key to continue . . ."); getchar ();


													/* Convert poles and zeros to coefficients */
	SDA_IirZplaneToCoeffs (ZPlaneZeros,				/* Z-plane zeros */
			ZPlanePoles,							/* Z-plane zeros */
			pIIRCoeffs,								/* IIR filter coefficients */
			NUMBER_OF_ZEROS,						/* Number of zeros */
			NUMBER_OF_POLES);						/* Number of poles */

#if	PRINT_RESULTS
	printf ("\nIIR filter coefficients\n");			/* Print filter coefficients */
	printf ("\nb(0), b(1), b(2), a(1), a(2)\n");
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
	{
		printf ("%le, %le, %le, %le, %le\n", *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+1), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+2),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+3), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+4));
	}
	printf ("Hit any key to see impulse and frequency response . . ."); getchar();
#endif

													/* Generate impulse response */
	h2DGraph =										/* Initialize graph */
		Create2DGraph ("Bilinear Transform IIR Filter Design",		/* 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);
	}

	SIF_Iir (pFilterState,							/* Pointer to filter state array */
			 FILTER_STAGES);						/* Number of second order stages */
													/* 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 */

				/* Generate test impulse */
	SDA_SignalGenerate (pSrc,						/* Pointer to destination array */
						SIGLIB_IMPULSE,				/* Signal type - Impulse function */
						SIGLIB_ONE,					/* Signal peak level */
						SIGLIB_FILL,				/* 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 */
						IMPULSE_RESPONSE_LENGTH);	/* Output array length */

													/* Apply iir filter and store filtered data */
	SDA_Iir (pSrc,									/* Input array to be filtered */
			 pRealData,								/* Filtered output array */
			 pFilterState,							/* Pointer to filter state array */
			 pIIRCoeffs,							/* Pointer to filter coefficients array */
			 FILTER_STAGES,							/* Number of stages */
			 IMPULSE_RESPONSE_LENGTH);				/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Impulse Response",				/* Title of the dataset */
				    pRealData,						/* Array of Double dataset */
				    PLOT_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 ("\nImpulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();

													/* Generate frequency response */
	SDA_Rfft (pRealData,							/* 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 (pRealData,					/* Pointer to real source array */
					  pImagData,					/* Pointer to imaginary source array */
					  pResults,						/* Pointer to log magnitude destination array */
					  PLOT_LENGTH);					/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Frequency Response",			/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    PLOT_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 ("\nFrequency Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();


#if	TRANSFORM_L_L				/* Low pass to low pass transformation */
	
				/* Transform the filter response and re-display */
	SDA_IirZplaneLpfToLpf (ZPlaneZeros,				/* Source Z-plane zeros */
						   ZPlanePoles,				/* Source Z-plane poles */
						   TransformedZPlaneZeros,	/* Destination Z-plane zeros */
						   TransformedZPlanePoles,	/* Destination Z-plane poles */
						   PREWARP_MATCH_FREQUENCY,	/* Source cut-off frequency */
						   LPF_TO_LPF_NEW_FC,		/* Destination cut-off frequency */
						   SAMPLE_RATE,				/* System sample rate */
						   NUMBER_OF_ZEROS,			/* Number of zeros in input and output */
						   NUMBER_OF_POLES);		/* Number of poles in input and output */

#if	PRINT_RESULTS
	printf ("\nTransformed complex z-plane zeros\n");	/* Print z-plane poles and zeros */
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
		print_rectangular (TransformedZPlaneZeros[i]);

	printf ("\nTransformed complex z-plane poles\n");
	for (i = 0; i < NUMBER_OF_POLES; i++)
		print_rectangular (TransformedZPlanePoles[i]);
	printf ("Hit any key to see pole-zero plot . . ."); getchar();
#endif

	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Poles",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlanePoles,	/* Array of complex values */
				  NUMBER_OF_POLES,					/* Items of data */
				  SV_BLUE,							/* Colour */
				  SV_GRAPH_NEW,						/* Add / new mode */
				  SV_CONJUGATE_POLE);				/* Pole-zero mode */
	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Zeros",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlaneZeros,	/* Array of complex values */
				  NUMBER_OF_ZEROS,					/* Items of data */
				  SV_RED,							/* Colour */
				  SV_GRAPH_ADD,						/* Add / new mode */
				  SV_CONJUGATE_ZERO);				/* Pole-zero mode */
	printf ("\nPole-zero Plot\nHit any key to continue . . ."); getchar ();


													/* Convert poles and zeros to coefficients */
	SDA_IirZplaneToCoeffs (TransformedZPlaneZeros,	/* Z-plane zeros */
						   TransformedZPlanePoles,	/* Z-plane zeros */
						   pIIRCoeffs,				/* IIR filter coefficients */
						   NUMBER_OF_ZEROS,			/* Number of zeros */
						   NUMBER_OF_POLES);		/* Number of poles */

#if	PRINT_RESULTS
	printf ("\nIIR filter coefficients\n");			/* Print filter coefficients */
	printf ("\nb(0), b(1), b(2), a(1), a(2)\n");
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
	{
		printf ("%le, %le, %le, %le, %le\n", *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+1), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+2),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+3), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+4));
	}
	printf ("\nHit any key to see transformed impulse and frequency response . . ."); getchar();
#endif

													/* Re-initialise filter state array */
	SIF_Iir (pFilterState,							/* Pointer to filter state array */
			 FILTER_STAGES);						/* Number of second order stages */

				/* Generate test impulse */
	SDA_SignalGenerate (pSrc,						/* Pointer to destination array */
						SIGLIB_IMPULSE,				/* Signal type - Impulse function */
						SIGLIB_ONE,					/* Signal peak level */
						SIGLIB_FILL,				/* 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 */
						IMPULSE_RESPONSE_LENGTH);	/* Output array length */

													/* Apply iir filter and store filtered data */
	SDA_Iir (pSrc,									/* Input array to be filtered */
			 pRealData,								/* Filtered output array */
			 pFilterState,							/* Pointer to filter state array */
			 pIIRCoeffs,							/* Pointer to filter coefficients array */
			 FILTER_STAGES,							/* Number of stages */
			 IMPULSE_RESPONSE_LENGTH);				/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Transformed Impulse Response",	/* Title of the dataset */
				    pRealData,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();

													/* Save the old frequency response */
	SDA_Copy (pResults,								/* Pointer to source array */
			  pSrc,									/* Pointer to destination array */
			  FFT_SIZE);							/* Array length */

													/* Generate frequency response */
													/* Perform real FFT */
	SDA_Rfft (pRealData,							/* 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 (pRealData,					/* Pointer to real source array */
					  pImagData,					/* Pointer to imaginary source array */
					  pResults,						/* Pointer to log magnitude destination array */
					  PLOT_LENGTH);					/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Original Frequency Response",	/* Title of the dataset */
				    pSrc,							/* Array of Double dataset */
				    PLOT_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 */
				    "Transformed Frequency Response",	/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed impulse response\nPlease hit <Carriage Return> to continue . . ."); getchar ();
#endif

#if	TRANSFORM_L_H				/* Low pass to high pass transformation */
	
				/* Transform the filter response and re-display */
	SDA_IirZplaneLpfToHpf (ZPlaneZeros,				/* Source Z-plane zeros */
						   ZPlanePoles,				/* Source Z-plane poles */
						   TransformedZPlaneZeros,	/* Destination Z-plane zeros */
						   TransformedZPlanePoles,	/* Destination Z-plane poles */
						   PREWARP_MATCH_FREQUENCY,	/* Source cut-off frequency */
						   LPF_TO_LPF_NEW_FC,		/* Destination cut-off frequency */
						   SAMPLE_RATE,				/* System sample rate */
						   NUMBER_OF_ZEROS,			/* Number of zeros in input and output */
						   NUMBER_OF_POLES);		/* Number of poles in input and output */

#if	PRINT_RESULTS
	printf ("\nTransformed complex z-plane zeros\n");	/* Print z-plane poles and zeros */
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
		print_rectangular (TransformedZPlaneZeros[i]);

	printf ("\nTransformed complex z-plane poles\n");
	for (i = 0; i < NUMBER_OF_POLES; i++)
		print_rectangular (TransformedZPlanePoles[i]);
	printf ("Hit any key to see pole-zero plot . . ."); getchar();
#endif

	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Poles",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlanePoles,	/* Array of complex values */
				  NUMBER_OF_POLES,					/* Items of data */
				  SV_BLUE,							/* Colour */
				  SV_GRAPH_NEW,						/* Add / new mode */
				  SV_CONJUGATE_POLE);				/* Pole-zero mode */
	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Zeros",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlaneZeros,	/* Array of complex values */
				  NUMBER_OF_ZEROS,					/* Items of data */
				  SV_RED,							/* Colour */
				  SV_GRAPH_ADD,						/* Add / new mode */
				  SV_CONJUGATE_ZERO);				/* Pole-zero mode */
	printf ("Pole-zero Plot\nHit any key to continue . . ."); getchar ();


													/* Convert poles and zeros to coefficients */
	SDA_IirZplaneToCoeffs (TransformedZPlaneZeros,	/* Z-plane zeros */
						   TransformedZPlanePoles,	/* Z-plane zeros */
						   pIIRCoeffs,				/* IIR filter coefficients */
						   NUMBER_OF_ZEROS,			/* Number of zeros */
						   NUMBER_OF_POLES);		/* Number of poles */

#if	PRINT_RESULTS
	printf ("\nIIR filter coefficients\n");			/* Print filter coefficients */
	printf ("\nb(0), b(1), b(2), a(1), a(2)\n");
	for (i = 0; i < NUMBER_OF_ZEROS; i++)
	{
		printf ("%le, %le, %le, %le, %le\n", *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+1), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+2),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+3), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+4));
	}
	printf ("Hit any key to see transformed impulse and frequency response . . ."); getchar();
#endif

													/* Re-initialise filter state array */
	SIF_Iir (pFilterState,							/* Pointer to filter state array */
			 FILTER_STAGES);						/* Number of second order stages */

				/* Generate test impulse */
	SDA_SignalGenerate (pSrc,						/* Pointer to destination array */
						SIGLIB_IMPULSE,				/* Signal type - Impulse function */
						SIGLIB_ONE,					/* Signal peak level */
						SIGLIB_FILL,				/* 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 */
						IMPULSE_RESPONSE_LENGTH);	/* Output array length */

													/* Apply iir filter and store filtered data */
	SDA_Iir (pSrc,									/* Input array to be filtered */
			 pRealData,								/* Filtered output array */
			 pFilterState,							/* Pointer to filter state array */
			 pIIRCoeffs,							/* Pointer to filter coefficients array */
			 FILTER_STAGES,							/* Number of stages */
			 IMPULSE_RESPONSE_LENGTH);				/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Transformed Impulse Response",	/* Title of the dataset */
				    pRealData,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();


													/* Generate frequency response */
													/* Perform real FFT */
	SDA_Rfft (pRealData,							/* 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 (pRealData,					/* Pointer to real source array */
					  pImagData,					/* Pointer to imaginary source array */
					  pResults,						/* Pointer to log magnitude destination array */
					  PLOT_LENGTH);					/* Array length */
	
	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Original Frequency Response",	/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();
#endif

#if	TRANSFORM_L_P				/* Low pass to band pass transformation */
	
				/* Transform the filter response and re-display */
	SDA_IirZplaneLpfToBpf (ZPlaneZeros,				/* Source Z-plane zeros */
						   ZPlanePoles,				/* Source Z-plane poles */
						   TransformedZPlaneZeros,	/* Destination Z-plane zeros */
						   TransformedZPlanePoles,	/* Destination Z-plane poles */
						   PREWARP_MATCH_FREQUENCY,	/* Source cut-off frequency */
						   LPF_TO_BPF_NEW_FC1,		/* Destination lower cut-off frequency */
						   LPF_TO_BPF_NEW_FC2,		/* Destination upper cut-off frequency */
						   SAMPLE_RATE,				/* System sample rate */
						   NUMBER_OF_ZEROS,			/* Number of zeros in input and output */
						   NUMBER_OF_POLES);		/* Number of poles in input and output */

#if	PRINT_RESULTS
	printf ("\nTransformed complex z-plane zeros\n");		/* Print z-plane poles and zeros */
	for (i = 0; i < 2*NUMBER_OF_ZEROS; i++)
		print_rectangular (TransformedZPlaneZeros[i]);

	printf ("\nTransformed complex z-plane poles\n");
	for (i = 0; i < 2*NUMBER_OF_POLES; i++)
		print_rectangular (TransformedZPlanePoles[i]);
	printf ("Hit any key to see pole-zero plot . . ."); getchar();
#endif

	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Poles",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlanePoles,	/* Array of complex values */
				  NUMBER_OF_POLES,					/* Items of data */
				  SV_BLUE,							/* Colour */
				  SV_GRAPH_NEW,						/* Add / new mode */
				  SV_CONJUGATE_POLE);				/* Pole-zero mode */
	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Zeros",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlaneZeros,	/* Array of complex values */
				  NUMBER_OF_ZEROS,					/* Items of data */
				  SV_RED,							/* Colour */
				  SV_GRAPH_ADD,						/* Add / new mode */
				  SV_CONJUGATE_ZERO);				/* Pole-zero mode */
	printf ("Pole-zero Plot\nHit any key to continue . . ."); getchar ();


													/* Convert poles and zeros to coefficients */
	SDA_IirZplaneToCoeffs (TransformedZPlaneZeros,	/* Z-plane zeros */
						   TransformedZPlanePoles,	/* Z-plane zeros */
						   pIIRCoeffs,				/* IIR filter coefficients */
						   2*NUMBER_OF_ZEROS,		/* Number of zeros */
						   2*NUMBER_OF_POLES);		/* Number of poles */

#if	PRINT_RESULTS
	printf ("\nIIR filter coefficients\n");			/* Print filter coefficients */
	printf ("\nb(0), b(1), b(2), a(1), a(2)\n");
	for (i = 0; i < 2*NUMBER_OF_ZEROS; i++)
	{
		printf ("%le, %le, %le, %le, %le\n", *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+1), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+2),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+3), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+4));
	}
	printf ("\nHit any key to see transformed impulse and frequency response . . ."); getchar();
#endif

													/* Re-initialise filter state array */
	SIF_Iir (pFilterState,							/* Pointer to filter state array */
			 2L*FILTER_STAGES);						/* Number of second order stages */

													/* Generate test impulse */
	SDA_SignalGenerate (pSrc,						/* Pointer to destination array */
						SIGLIB_IMPULSE,				/* Signal type - Impulse function */
						SIGLIB_ONE,					/* Signal peak level */
						SIGLIB_FILL,				/* 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 */
						IMPULSE_RESPONSE_LENGTH);	/* Output array length */

													/* Apply iir filter and store filtered data */
	SDA_Iir (pSrc,									/* Input array to be filtered */
			 pRealData,								/* Filtered output array */
			 pFilterState,							/* Pointer to filter state array */
			 pIIRCoeffs,							/* Pointer to filter coefficients array */
			 2L*FILTER_STAGES,						/* Number of stages */
			 IMPULSE_RESPONSE_LENGTH);				/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Transformed Impulse Response",	/* Title of the dataset */
				    pRealData,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();


													/* Generate frequency response */
													/* Perform real FFT */
	SDA_Rfft (pRealData,							/* 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 (pRealData,					/* Pointer to real source array */
					  pImagData,					/* Pointer to imaginary source array */
					  pResults,						/* Pointer to log magnitude destination array */
					  PLOT_LENGTH);					/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Original Frequency Response",	/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();
#endif

#if	TRANSFORM_L_S				/* Low pass to band stop transformation */
	
				/* Transform the filter response and re-display */
	SDA_IirZplaneLpfToBsf (ZPlaneZeros,				/* Source Z-plane zeros */
						   ZPlanePoles,				/* Source Z-plane poles */
						   TransformedZPlaneZeros,	/* Destination Z-plane zeros */
						   TransformedZPlanePoles,	/* Destination Z-plane poles */
						   PREWARP_MATCH_FREQUENCY,	/* Source cut-off frequency */
						   LPF_TO_BPF_NEW_FC1,		/* Destination lower cut-off frequency */
						   LPF_TO_BPF_NEW_FC2,		/* Destination upper cut-off frequency */
						   SAMPLE_RATE,				/* System sample rate */
						   NUMBER_OF_ZEROS,			/* Number of zeros in input and output */
						   NUMBER_OF_POLES);		/* Number of poles in input and output */

#if	PRINT_RESULTS
	printf ("\nTransformed complex z-plane zeros\n");	/* Print z-plane poles and zeros */
	for (i = 0; i < 2*NUMBER_OF_ZEROS; i++)
		print_rectangular (TransformedZPlaneZeros[i]);

	printf ("\nTransformed complex z-plane poles\n");
	for (i = 0; i < 2*NUMBER_OF_POLES; i++)
		print_rectangular (TransformedZPlanePoles[i]);
	printf ("Hit any key to see pole-zero plot . . ."); getchar();
#endif

	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Poles",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlanePoles,	/* Array of complex values */
				  NUMBER_OF_POLES,					/* Items of data */
				  SV_BLUE,							/* Colour */
				  SV_GRAPH_NEW,						/* Add / new mode */
				  SV_CONJUGATE_POLE);				/* Pole-zero mode */
	DisplayPZPlot(hPZGraph,							/* Graph object */
				  "Zeros",							/* Title of dataset (for use on legends) */
				  (Complex_s *)TransformedZPlaneZeros,	/* Array of complex values */
				  NUMBER_OF_ZEROS,					/* Items of data */
				  SV_RED,							/* Colour */
				  SV_GRAPH_ADD,						/* Add / new mode */
				  SV_CONJUGATE_ZERO);				/* Pole-zero mode */
	printf ("Pole-zero Plot\nHit any key to continue . . ."); getchar ();


													/* Convert poles and zeros to coefficients */
	SDA_IirZplaneToCoeffs (TransformedZPlaneZeros,	/* Z-plane zeros */
						   TransformedZPlanePoles,	/* Z-plane zeros */
						   pIIRCoeffs,				/* IIR filter coefficients */
						   2*NUMBER_OF_ZEROS,		/* Number of zeros */
						   2*NUMBER_OF_POLES);		/* Number of poles */

#if	PRINT_RESULTS
	printf ("\nIIR filter coefficients\n");			/* Print filter coefficients */
	printf ("\nb(0), b(1), b(2), a(1), a(2)\n");
	for (i = 0; i < 2*NUMBER_OF_ZEROS; i++)
	{
		printf ("%le, %le, %le, %le, %le\n", *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+1), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+2),
				*(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+3), *(pIIRCoeffs+(i*IIR_COEFFS_PER_BIQUAD)+4));
	}
	printf ("\nHit any key to see transformed impulse and frequency response . . ."); getchar();
#endif

													/* Re-initialise filter state array */
	SIF_Iir (pFilterState,							/* Pointer to filter state array */
			 2L*FILTER_STAGES);						/* Number of second order stages */

				/* Generate test impulse */
	SDA_SignalGenerate (pSrc,						/* Pointer to destination array */
						SIGLIB_IMPULSE,				/* Signal type - Impulse function */
						SIGLIB_ONE,					/* Signal peak level */
						SIGLIB_FILL,				/* 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 */
						IMPULSE_RESPONSE_LENGTH);	/* Output array length */

													/* Apply iir filter and store filtered data */
	SDA_Iir (pSrc,									/* Input array to be filtered */
			 pRealData,								/* Filtered output array */
			 pFilterState,							/* Pointer to filter state array */
			 pIIRCoeffs,							/* Pointer to filter coefficients array */
			 2L*FILTER_STAGES,						/* Number of stages */
			 IMPULSE_RESPONSE_LENGTH);				/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Transformed Impulse Response",	/* Title of the dataset */
				    pRealData,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\nPlease hit <Carriage Return> to continue . . ."); getchar ();

													/* Generate frequency response */
													/* Perform real FFT */
	SDA_Rfft (pRealData,							/* 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 (pRealData,					/* Pointer to real source array */
					  pImagData,					/* Pointer to imaginary source array */
					  pResults,						/* Pointer to log magnitude destination array */
					  PLOT_LENGTH);					/* Array length */

	Display2DGraph (h2DGraph,						/* Graph handle */
				    "Transformed Frequency Response",	/* Title of the dataset */
				    pResults,						/* Array of Double dataset */
				    PLOT_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 ("\nTransformed Impulse Response\n");
#endif

	SUF_MemoryFree (pSrc);							/* Free memory */
	SUF_MemoryFree (pIIRCoeffs);
	SUF_MemoryFree (pFilterState);
	SUF_MemoryFree (pRealData);
	SUF_MemoryFree (pImagData);
	SUF_MemoryFree (pResults);
	SUF_MemoryFree (pFFTCoeffs);
}