/* SigLib - Program to read .dat file and generate a waterfall diagram */


#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <siglib.h>			/* SigLib DSP library */
#include "GraphFunctions.h"
#include "readdat.h"

#define SAMPLE_LENGTH		1024L					/* Length of array read from input file */
#define FFT_LENGTH			1024L					/* Length of FFT performed */
#define LOG2_FFT_LENGTH		10L
#define RESULT_LENGTH		(FFT_LENGTH >> 1)		/* Only need to store the lower 1/2 of the FFT output */
#define OVERLAP_LENGTH		(SAMPLE_LENGTH >> 2)	/* 25 % overlap */

#define GRAPH_X_AXIS_LENGTH		((SLArrayIndex_t)950)

#define SAMPLE_RATE			10162					/* Sample rate of data in vibration data file */
#define NUMBER_OF_SAMPLES	101547					/* Number of samples in vibration data file */

#define DB_SCALE			(1e5)					/* Scaling for dB */

SLData_t	*pDataArray, *pFDPRealData, *pFDPImagData, *pFDPResults;
SLData_t	*pFDPFFTCoeffs, *pWindowCoeffs, *pOverlapArray;

char		FileName[80];

void main (int argc, char *argv[]);

void main (int argc, char *argv[])
{
	GraphObject *h3DGraph;							/* Declare graph object */

	SLArrayIndex_t	SampleCount;
	FILE			*pInputFile;
	SLArrayIndex_t	i;
	SLArrayIndex_t	FrameNumber = 0L;
	SLArrayIndex_t	OverlapSrcArrayIndex;

	SLData_t		Speed;
	SLData_t		SampleRate;
	SLArrayIndex_t	TotalNumberOfSamples;

	SLData_t		WindowInverseCoherentGain = SIGLIB_ONE;

	pDataArray = SUF_VectorArrayAllocate (SAMPLE_LENGTH);		/* Input data array */
	pOverlapArray = SUF_VectorArrayAllocate (OVERLAP_LENGTH);	/* Overlap data array */
	pWindowCoeffs = SUF_VectorArrayAllocate (FFT_LENGTH);		/* Window coeffs data array */
	pFDPRealData = SUF_VectorArrayAllocate (FFT_LENGTH);		/* Real data array */
	pFDPImagData = SUF_VectorArrayAllocate (FFT_LENGTH);		/* Imaginary data array */
	pFDPResults = SUF_VectorArrayAllocate (FFT_LENGTH);			/* Results data array */
	pFDPFFTCoeffs = SUF_FftCoefficientAllocate (FFT_LENGTH);	/* FFT coefficient data array */

	if ((pDataArray == NULL) || (pOverlapArray == NULL) || (pWindowCoeffs == NULL) ||
		(pFDPRealData == NULL) || (pFDPImagData == NULL) ||
		(pFDPResults == NULL) || (pFDPFFTCoeffs == NULL))
	{
		printf ("Memory allocation error\n");
		exit (0);
	}

	SIF_CopyWithOverlap (&OverlapSrcArrayIndex);	/* Pointer to source array index */

													/* Generate window table */
	SIF_Window (pWindowCoeffs,						/* Window coefficients pointer */
				SIGLIB_BMAN_HARRIS,					/* Window type */
				SIGLIB_ZERO,						/* Window coefficient */
				FFT_LENGTH);						/* Window length */

													/* Calculate window inverse coherent gain */
	WindowInverseCoherentGain =
		SDA_WindowInverseCoherentGain(pWindowCoeffs,	/* Source array pointer */
									  FFT_LENGTH);		/* Window size */

													/* Initialise FFT */
	SIF_Fft (pFDPFFTCoeffs,							/* Pointer to FFT coefficients */
			 SIGLIB_NULL_ARRAY_INDEX_PTR,			/* Pointer to bit reverse address table */
			 FFT_LENGTH);							/* FFT Size */

	if (argc != 2)
	{
		printf ("\nUsage error  :\nspecgram filename\n\n");
		exit (1);									/* Exit - usage error */
	}

	strcpy (FileName, argv[1]);

	printf ("Source file = %s\n", FileName);

	if ((pInputFile = fopen(FileName, "r")) == NULL)	/* Note this file is text */
	{
		printf ("Error opening input vibration data file\n");
		exit (1);
	}

	h3DGraph =										/* Initialize graph */
		Create3DGraph ("Machine Spectrogram",		/* Graph title */
					   "Time",						/* X-Axis label */
					   "Frequency",					/* Y-Axis label */
					   GRAPH_X_AXIS_LENGTH,			/* X-axis length */
					   0.0,							/* Minimum signal magnitude */
					   80.0,						/* Maximum signal magnitude */
					   0.0,							/* Y-Axis minimum value */
					   SampleRate / SIGLIB_TWO,		/* Y-Axis maximum value */
					   SV_3D_COLOUR,				/* Graph mode */
					   "localhost");				/* Graph server */
	if (h3DGraph == NULL)							/* Graph creation failed - e.g is server running ? */
	{
		printf ("\nGraph creation failure. Please check that the server is running\n");
		exit (1);
	}

	SampleRate = (SLData_t)SAMPLE_RATE;
	TotalNumberOfSamples = NUMBER_OF_SAMPLES;


	while ((SampleCount =
				read_vibration_data (pDataArray,
									 pInputFile,
									 &Speed,
									 SAMPLE_LENGTH))
					== SAMPLE_LENGTH)
	{
//		printf ("Speed = %lf, Data = %lf\n", Speed, *pDataArray);

													/* Apply the overlap to the data */
		while (SDA_CopyWithOverlap (pDataArray,				/* Pointer to source array */
									pFDPRealData,			/* Pointer to destination array */
									pOverlapArray,			/* Pointer to overlap array */
									&OverlapSrcArrayIndex,	/* Pointer to source array index */
									SAMPLE_LENGTH,			/* Source array length */
									OVERLAP_LENGTH,			/* Overlap length */
									SAMPLE_LENGTH) <		/* Destination array length */
										SAMPLE_LENGTH)
		{

//SUF_Debugfprintf ("\nFN = %ld\n", FrameNumber);
//SUF_DebugPrintArray (pDataArray, SAMPLE_LENGTH);

													/* Apply window to real data */
			SDA_Window (pFDPRealData,				/* Source array pointer */
						pFDPRealData,				/* Destination array pointer */
						pWindowCoeffs,				/* Window array pointer */
						FFT_LENGTH);  				/* Window size */
													/* Perform FFT */
			SDA_Rfft (pFDPRealData,					/* Real array pointer */
					  pFDPImagData,					/* Pointer to imaginary array */
					  pFDPFFTCoeffs,				/* Pointer to FFT coefficients */
					  SIGLIB_NULL_ARRAY_INDEX_PTR,	/* Pointer to bit reverse address table */
					  FFT_LENGTH,					/* FFT size */
					  LOG2_FFT_LENGTH);				/* log2 FFT size */

													/* Normalize FFT scaling value */
			SDA_ComplexScalarMultiply (pFDPRealData,	/* Pointer to real source array */
									   pFDPImagData,	/* Pointer to imaginary source array */
						  			   (DB_SCALE * SIGLIB_TWO * WindowInverseCoherentGain)/((SLData_t)FFT_LENGTH),	/* Multiplier */
									   pFDPRealData,	/* Pointer to real destination array */
									   pFDPImagData,	/* Pointer to imaginary destination array */
									   FFT_LENGTH);		/* Array lengths */

													/* Calc real power fm complex */
			SDA_LogMagnitude (pFDPRealData,			/* Pointer to real source array */
							  pFDPImagData,			/* Pointer to imaginary source array */
							  pFDPResults,			/* Pointer to log magnitude destination array */
							  RESULT_LENGTH);		/* Array length */
													/* Clip off noise */
			SDA_Offset (pFDPResults,				/* Pointer to source array */
						-50.0,						/* Offset */
						pFDPResults,				/* Pointer to destination array */
						RESULT_LENGTH);				/* Array length */

			SDA_Clip (pFDPResults,					/* Source array address */
					  pFDPResults,					/* Destination array address */
					  SIGLIB_ZERO,					/* Value to clip signal to */
					  SIGLIB_CLIP_BELOW,			/* Clip type */
					  RESULT_LENGTH);				/* Array length */

													/* Plot results */
			Display3DGraph (h3DGraph,				/* Graph handle */
							"Machine Spectrogram",	/* Title of the dataset */
							pFDPResults,			/* Array of Double dataset */
							RESULT_LENGTH);			/* Number of data points */


			FrameNumber++;
		}
	}

	fclose (pInputFile);							/* Close files */

	free (pDataArray);								/* Free memory */
	free (pOverlapArray);
	free (pWindowCoeffs);
	free (pFDPRealData);
	free (pFDPImagData);
	free (pFDPResults);
	free (pFDPFFTCoeffs);
}