/* 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. This program uses the chirp z-transform and the complex FFT and compares the two results. */ /* Include files */ #include #include #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)64) /* Input array length */ #define OUTPUT_LENGTH ((SLArrayIndex_t)128) #define FFT_SIZE 256L #define LOG_FFT_SIZE ((SLArrayIndex_t)8) #define FT_SIZE OUTPUT_LENGTH /* Declare global variables */ int waveform; SLData_t Radius, Decay, StartFreq, EndFreq; SLComplexPolar_s ContourStart; SLData_t phinc, w1inc, w2inc; /* Variables used to calculate W */ SLComplexRect_s A_1, W1, W_1, W12, W_12; /* Complex contour coeffs */ SLData_t *pInput, *pResults; SLData_t *pRealData, *pImagData, *pFFTCoeffs; SLData_t *pAWNr, *pAWNi, *pvLr, *pvLi, *pWMr, *pWMi; SLData_t *pCZTRealWork, *pCZTImagWork; SLData_t *pRealDataFT, *pImagDataFT; SLData_t SinePhase; void main(void) { GraphObject *h2DGraph; /* Declare graph object */ pInput = SUF_VectorArrayAllocate (FFT_SIZE); pRealData = SUF_VectorArrayAllocate (FFT_SIZE); pImagData = SUF_VectorArrayAllocate (FFT_SIZE); pRealDataFT = SUF_VectorArrayAllocate (FFT_SIZE); pImagDataFT = SUF_VectorArrayAllocate (FFT_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) || (pRealDataFT == NULL) || (pImagDataFT == 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); } 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); } SinePhase = SIGLIB_ZERO; SDA_SignalGenerate (pInput, /* Pointer to destination array */ SIGLIB_SINE_WAVE, /* Signal type - Sine wave */ SIGLIB_ONE, /* Signal peak level */ SIGLIB_FILL, /* Fill (overwrite) or add to existing array contents */ ((SLData_t)0.01562), /* 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 */ FFT_SIZE); /* Output array length */ SDA_Clear (pRealData, /* Pointer to destination array */ FT_SIZE); /* Array length */ SDA_Copy (pInput, /* Pointer to source array */ pRealData, /* Pointer to destination array */ INPUT_LENGTH); /* Array length */ /* Perform DFT */ SDA_Rft (pRealData, /* Pointer to real source array */ pRealDataFT, /* Pointer to real destination array */ pImagDataFT, /* Pointer to imaginary destination array */ FT_SIZE); /* Array length */ /* 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 */ Radius = SIGLIB_ONE; Decay = SIGLIB_ZERO; StartFreq = SIGLIB_ZERO; EndFreq = SAMPLE_RATE; /* Calculate A0 values */ /* Create a polar number from magnitude and angle */ ContourStart = SCV_Polar (Radius, (SIGLIB_TWO_PI * StartFreq) / SAMPLE_RATE); A_1 = SCV_Inverse (SCV_PolarToRectangular (ContourStart)); /* Calculate W0 values - W^(N^2/2) */ phinc = (-SIGLIB_TWO_PI * (EndFreq - StartFreq)) / (((SLData_t)OUTPUT_LENGTH) * SAMPLE_RATE); w1inc = SDS_Exp ((SIGLIB_PI * Decay) / SAMPLE_RATE); w2inc = SDS_Sqrt (w1inc); W1.real = w1inc * SDS_Cos (phinc); /* W */ W1.imag = w1inc * SDS_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, w1inc = %lf, w2inc = %lf\n\n", phinc, w1inc, w2inc); printf ("A_1.real = %lf, A_1.imag = %lf\n\n", A_1.real, A_1.imag); printf ("W1.real = %lf, W1.imag = %lf\n", W1.real, W1.imag); printf ("W12.real = %lf, W12.imag = %lf\n\n", W12.real, W12.imag); printf ("W_1.real = %lf, W_1.imag = %lf\n", W_1.real, W_1.imag); printf ("W_12.real = %lf, W_12.imag = %lf\n\n", W_12.real, W_12.imag); getchar (); #endif 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 */ "vLr", /* Title of the dataset */ pvLr, /* Array of Double dataset */ FFT_SIZE, /* Number of data points */ SV_GRAPH_LINE, /* Graph type */ SV_BLUE, /* Colour */ SV_HIDE_MARKERS, /* Marker enable / disable */ SV_GRAPH_NEW); /* New graph */ printf ("\nvLr\nPlease hit to continue . . ."); getchar (); Display2DGraph (h2DGraph, /* Graph handle */ "vLi", /* Title of the dataset */ pvLi, /* Array of Double dataset */ FFT_SIZE, /* Number of data points */ SV_GRAPH_LINE, /* Graph type */ SV_BLUE, /* Colour */ SV_HIDE_MARKERS, /* Marker enable / disable */ SV_GRAPH_NEW); /* New graph */ printf ("\nvLi\nPlease hit to continue . . ."); getchar (); /* VL data FFT */ SDA_Cfft (pvLr, /* Pointer to real array */ pvLi, /* 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 */ /* Ensure zero padded samples */ SDA_Clear (pCZTRealWork+INPUT_LENGTH, /* Pointer to destination array */ FFT_SIZE-INPUT_LENGTH); /* Array length */ SDA_Clear (pCZTImagWork+INPUT_LENGTH, /* Pointer to destination array */ FFT_SIZE-INPUT_LENGTH); /* Array length */ /* Complex window = complex mpy with real input 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 */ 2.0 * ((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 */ Display2DGraph (h2DGraph, /* Graph handle */ "Chirp z-Transform Real Result", /* Title of the dataset */ pRealData, /* 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 Real Result\n"); Display2DGraph (h2DGraph, /* Graph handle */ "Fourier Transform Real Result", /* Title of the dataset */ pRealDataFT, /* Array of Double dataset */ OUTPUT_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 ("Fourier Transform Real Result\nPlease hit to continue . . ."); getchar (); Display2DGraph (h2DGraph, /* Graph handle */ "Chirp z-Transform Imaginary Result", /* Title of the dataset */ pImagData, /* 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 Imaginary Result\n"); Display2DGraph (h2DGraph, /* Graph handle */ "Fourier Transform Imaginary Result", /* Title of the dataset */ pImagDataFT, /* Array of Double dataset */ OUTPUT_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 ("Fourier Transform Imaginary Result\n"); SUF_MemoryFree (pInput); /* Free memory */ SUF_MemoryFree (pRealData); SUF_MemoryFree (pImagData); SUF_MemoryFree (pRealDataFT); SUF_MemoryFree (pImagDataFT); 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); }