Chirp-Z transform - source code needed

[loucy]

Chirp-Z transform

Anyone know of a source code for doing chirp-z transform? 2-D case?

Thanks!

 

[btbass]

Try to find siglib, loads of dsp source code, might have what you are looking for.
SigLib sample

/**************************************************************************
File Name : CHIRPZ.C | Author : JOHN EDWARDS
Siglib Library version : 6.00 | Modifications :
----------------------------------------+----------------------------------
Compiler : Independent | Start Date : 16/07/94
Options : | Latest Update : 02/06/01
---------------------------------------------------------------------------
Support for SigLib is available via Email : support@numerix-dsp.com

This file may be modified, in any way, providing that this header remains
within the file and the only additions are made to this header.

THIS FILE MAY ONLY BE USED IN ACCORDANCE WITH THE TERMS OF THE NUMERIX LTD.
END USER LICENSE AGREEMENT.
DELETION OF ANY INFORMATION IN THIS HEADER IS IN VIOLATION OF YOUR LICENSE.

Copyright (C) 1994 to 2001 Numerix Ltd.
---------------------------------------------------------------------------
Description : Chirp z-Transform routines, for SigLib DSP library.


****************************************************************************/

#define SIGLIB_SRC_FILE_CHIRPZ 1 /* Defines the source file that siglibms.h is being used in */

#include <siglib.h> /* Include SigLib header file */

#define DEBUG 0 /* Set to '1' to enable debug print outs */

/**/
/********************************************************
* Function : siglib_numerix_SIFCzt
*
* Parameters :
* SFLOAT * SIGLIB_PTR_DECL , - Pointer to AWNr coefficients
* SFLOAT * SIGLIB_PTR_DECL , - Pointer to AWNi coefficients
* SFLOAT * SIGLIB_PTR_DECL , - Pointer to WMr coefficients
* SFLOAT * SIGLIB_PTR_DECL , - Pointer to WMi coefficients
* SFLOAT * SIGLIB_PTR_DECL , - Pointer to vLr coefficients
* SFLOAT * SIGLIB_PTR_DECL , - Pointer to vLi coefficients
* SFLOAT * SIGLIB_PTR_DECL , - FFT coeff pointer
* SFIX * SIGLIB_PTR_DECL , - Bit reverse address table pointer
* const SFLOAT Radius, - Contour start radius
* const SFLOAT Decay, - Contour decay rate
* const SFLOAT StartFreq, - Contour start frequency
* const SFLOAT EndFreq, - Contour end frequency
* const SFLOAT SampleRate, - System sample rate
* const SFIX InLen, - Input buffer length
* const SFIX OutLen, - Output buffer length
* const SFIX FFTSize, - FFT size
* const SFIX Log2FFTSize - Log 2 FFT size
*
* Return value :
* Error code
*
* Description :
* Initialise the chirp z-transform functions.
*
********************************************************/

#if SIGLIB_INCLUDE_FN_siglib_numerix_SIFCzt

SFIX SIGLIB_FUNC_DECL siglib_numerix_SIFCzt (SFLOAT * SIGLIB_PTR_DECL pAWNr,
SFLOAT * SIGLIB_PTR_DECL pAWNi,
SFLOAT * SIGLIB_PTR_DECL pWMr,
SFLOAT * SIGLIB_PTR_DECL pWMi,
SFLOAT * SIGLIB_PTR_DECL pvLr,
SFLOAT * SIGLIB_PTR_DECL pvLi,
SFLOAT * SIGLIB_PTR_DECL pFFTCoeffs,
SFIX * SIGLIB_PTR_DECL pBitReverseAddressTable,
const SFLOAT Radius,
const SFLOAT Decay,
const SFLOAT StartFreq,
const SFLOAT EndFreq,
const SFLOAT SampleRate,
const SFIX InLen,
const SFIX OutLen,
const SFIX FFTSize,
const SFIX Log2FFTSize)

{
ComplexPolar_s contour_start;
SFLOAT deltaomega, deltasigma; /* Variables used to calculate W */
SFLOAT phinc, w1inc, w2inc;

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

SFIX SigLibErrorCode;


/* Calculate A0 values */
contour_start = siglib_numerix_SCVPolar (Radius, SIGLIB_TWO_PI * (StartFreq / SampleRate));
A_1 = siglib_numerix_SCVInverse (siglib_numerix_SCVPolarToRectangular (contour_start));

/* Calculate W0 values - W^(N^2/2) */
deltaomega = (EndFreq - StartFreq) / ((SFLOAT)OutLen);
deltasigma = Decay;

phinc = SIGLIB_TWO * (-SIGLIB_PI) * deltaomega / SampleRate;

w1inc = siglib_numerix_SDSExp (SIGLIB_PI * deltasigma / SampleRate);
w2inc = siglib_numerix_SDSSqrt (w1inc);

W1.real = w1inc * siglib_numerix_SDSCos (phinc); /* W */
W1.imag = w1inc * siglib_numerix_SDSSin (phinc);

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

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

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


SigLibErrorCode = siglib_numerix_SIFFft (pFFTCoeffs, pBitReverseAddressTable, FFTSize); /* Initialise FFT */
if (SigLibErrorCode != SIGLIB_NO_ERROR)
{
return (SigLibErrorCode);
}


siglib_numerix_SIFAwn (pAWNr, pAWNi, A_1, W1, W12, InLen); /* Gen. complex window coeffs */
siglib_numerix_SIFVl (pvLr, pvLi, W_1, W_12, InLen, OutLen, FFTSize); /* Gen. contour definition coeffs */
siglib_numerix_SIFWm (pWMr, pWMi, W1, W12, OutLen); /* Gen. weighting coeffs */


/* Performing this here will save time in the real time process */
siglib_numerix_SDACfft (pvLr, pvLi, pFFTCoeffs, SIGLIB_NULL_SFIX_PTR, FFTSize, Log2FFTSize); /* Contour coeffs. FFT */


#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);
getch ();
#endif

return (SIGLIB_NO_ERROR);

} /* End of siglib_numerix_SIFCzt () */

#endif /* End of conditional include */




/**/
/********************************************************
* Function : siglib_numerix_SIFAwn
*
* Parameters :
* SFLOAT * SIGLIB_PTR_DECL pAWNr - Real coefficient pointer
* SFLOAT * SIGLIB_PTR_DECL pAWNi - Imaginary coefficient pointer
* const Complex, - A ^ (-1)
* const Complex, - W
* const Complex, - W^(1/2)
* const SFIX SampleLength - Buffer length
*
* Return value :
* void
*
* Description :
* Generate complex window coeffs
*
********************************************************/

#if SIGLIB_INCLUDE_FN_siglib_numerix_SIFAwn

void SIGLIB_FUNC_DECL siglib_numerix_SIFAwn (SFLOAT * SIGLIB_PTR_DECL pAWNr,
SFLOAT * SIGLIB_PTR_DECL pAWNi,
const ComplexRect_s A_1,
const ComplexRect_s W1,
const ComplexRect_s W12,
const SFIX SampleLength)

{
register SFIX i;
ComplexRect_s C, D, Temp;

C.real = SIGLIB_ONE;
C.imag = SIGLIB_ZERO;

*pAWNr++ = SIGLIB_ONE;
*pAWNi++ = SIGLIB_ZERO;

D.real = (W12.real * A_1.real) - (W12.imag * A_1.imag);
D.imag = (W12.real * A_1.imag) + (W12.imag * A_1.real);

for (i = 1; i < SampleLength; i++) /* Iterate complex chirp (sin & cos) */
{
Temp.real = (C.real * D.real) - (C.imag * D.imag);
Temp.imag = (C.real * D.imag) + (C.imag * D.real);

C.real = Temp.real;
C.imag = Temp.imag;

*pAWNr++ = Temp.real;
*pAWNi++ = Temp.imag;

Temp.real = (D.real * W1.real) - (D.imag * W1.imag);
Temp.imag = (D.real * W1.imag) + (D.imag * W1.real);

D.real = Temp.real;
D.imag = Temp.imag;
}

} /* End of siglib_numerix_SIFAwn () */

#endif /* End of conditional include */



/**/
/********************************************************
* Function : siglib_numerix_SIFVl
*
* Parameters :
* SFLOAT * SIGLIB_PTR_DECL pvLr - Real coefficient pointer
* SFLOAT * SIGLIB_PTR_DECL pvLi - Imaginary coefficient pointer
* const Complex, - W^(-1)
* const Complex, - W^(-1/2)
* const SFIX InLen - Input buffer length
* const SFIX OutLen - Output buffer length
* const SFIX FFTLen - FFT buffer length
*
* Return value :
* void
*
* Description :
* Generate contour definition coeffs
*
********************************************************/

#if SIGLIB_INCLUDE_FN_siglib_numerix_SIFVl

void SIGLIB_FUNC_DECL siglib_numerix_SIFVl (SFLOAT * SIGLIB_PTR_DECL pvLr,
SFLOAT * SIGLIB_PTR_DECL pvLi,
const ComplexRect_s W_1,
const ComplexRect_s W_12,
const SFIX InLen,
const SFIX OutLen,
const SFIX FFTLen)

{
register SFIX i;
ComplexRect_s C, D, Temp;

siglib_numerix_SDAClear (pvLr, FFTLen); /* Clear output buffer */
siglib_numerix_SDAClear (pvLi, FFTLen);

C.real = SIGLIB_ONE;
C.imag = SIGLIB_ZERO;

pvLr[0] = SIGLIB_ONE;
pvLi[0] = SIGLIB_ZERO;

D.real = W_12.real;
D.imag = W_12.imag;

for (i = 1; i < OutLen; i++) /* Iterate complex chirp (sin & cos) */
{
Temp.real = (C.real * D.real) - (C.imag * D.imag);
Temp.imag = (C.real * D.imag) + (C.imag * D.real);

C.real = Temp.real;
C.imag = Temp.imag;

pvLr = Temp.real;
pvLi = Temp.imag;

Temp.real = (D.real * W_1.real) - (D.imag * W_1.imag);
Temp.imag = (D.real * W_1.imag) + (D.imag * W_1.real);

D.real = Temp.real;
D.imag = Temp.imag;
}


pvLr += (FFTLen - InLen); /* Set ptrs to start of 2nd chirp - don't save first result */
pvLi += (FFTLen - InLen);

C.real = SIGLIB_ONE;
C.imag = SIGLIB_ZERO;

D.real = W_12.real;
D.imag = W_12.imag;

for (i = 0; i < InLen; i++) /* Iterate complex chirp (sin & cos) */
{
Temp.real = (C.real * D.real) - (C.imag * D.imag);
Temp.imag = (C.real * D.imag) + (C.imag * D.real);

C.real = Temp.real;
C.imag = Temp.imag;

pvLr[InLen-i-1] = Temp.real;
pvLi[InLen-i-1] = Temp.imag;

Temp.real = (D.real * W_1.real) - (D.imag * W_1.imag);
Temp.imag = (D.real * W_1.imag) + (D.imag * W_1.real);

D.real = Temp.real;
D.imag = Temp.imag;
}

} /* End of siglib_numerix_SIFVl () */

#endif /* End of conditional include */




/**/
/********************************************************
* Function : siglib_numerix_SIFWm
*
* Parameters :
* SFLOAT * SIGLIB_PTR_DECL pWMr, - Real coefficient pointer
* SFLOAT * SIGLIB_PTR_DECL pWMi, - Imaginary coefficient pointer
* const Complex, - W
* const Complex, - W^(1/2)
* const SFIX SampleLength - Buffer length
*
* Return value :
* void
*
* Description :
* Generate weighting coeffs
*
********************************************************/

#if SIGLIB_INCLUDE_FN_siglib_numerix_SIFWm

void SIGLIB_FUNC_DECL siglib_numerix_SIFWm (SFLOAT * SIGLIB_PTR_DECL pWMr,
SFLOAT * SIGLIB_PTR_DECL pWMi,
const ComplexRect_s W1,
const ComplexRect_s W12,
const SFIX SampleLength)

{
register SFIX i;
ComplexRect_s C, D, Temp;

C.real = SIGLIB_ONE;
C.imag = SIGLIB_ZERO;

*pWMr++ = SIGLIB_ONE;
*pWMi++ = SIGLIB_ZERO;

D.real = W12.real;
D.imag = W12.imag;

for (i = 1; i < SampleLength; i++) /* Iterate complex chirp (sin & cos) */
{
Temp.real = (C.real * D.real) - (C.imag * D.imag);
Temp.imag = (C.real * D.imag) + (C.imag * D.real);

C.real = Temp.real;
C.imag = Temp.imag;

*pWMr++ = Temp.real;
*pWMi++ = Temp.imag;

Temp.real = (D.real * W1.real) - (D.imag * W1.imag);
Temp.imag = (D.real * W1.imag) + (D.imag * W1.real);

D.real = Temp.real;
D.imag = Temp.imag;
}
} /* End of siglib_numerix_SIFWm () */

#endif /* End of conditional include */