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

	RUNFFTCPLX.C
   Z-World, 2000

	Sample test program for fftcplx()

	First construct a floating-point array comprising a dc
	value and up to three sinusoids with configurable amplitude,
	frequency, and phase. Round the 	floating values to
	integers, scaling if necessary to fit in a 16-bit word,
	store them in the integer array x[], and print x[].  Then
	transform x[] via fftcplx() and print the results. Finally
	inverse transform x[] via fftcplxinv() and print it a third
	time.

************************************************************/
#use fft.lib

int k, fsize, scale, blockexp, x[2048];
unsigned long int timer;

float fxmax, fscale, fx[2048];
float DC,
      A1, F1, P1,
      A2, F2, P2,
      A3, F3, P3;

#define     max(a,b)    ( a > b ? a : b )
#define     twoPI       (2*PI)

void dumpx( char *title );
void printsigparams( void );

void main ( void ) {

    blockexp = 0;

//  First construct a complex floating-point array comprising a dc value and up
//  to three sinusoids with configurable amplitude, frequency, and phase.  Round
//  the floating values to integers, scaling if necessary to fit in a 16-bit
//  word, store them in the integer array x[], and print x[].  Then transform
//  x[] via fftcplx() and print the results.  Finally inverse transform x[] via
//  fftcplxinv() and print it a third time.

	fsize = 16;

    DC =  2000.0;
    A1 =  4000.0;       F1 = 1.0/16;        P1 =  0.0;
    A2 =  8000.0;       F2 = 1.0/8;         P2 = PI/4;
    A3 = 16000.0;       F3 = 1.0/2;         P3 = PI/2;
    // Amplitude        // Frequency        // Phase

    // first calculate floating point values and note maximum value

    fxmax = 0;
    for (k = 0; k < fsize; k++) {
        fx[  2*k  ] =  DC
                   + A1 * sin(twoPI * F1 * k + P1)
                   + A2 * sin(twoPI * F2 * k + P2)
                   + A3 * sin(twoPI * F3 * k + P3);
        fx[2*k + 1] = 0.0;
        fxmax = max( fxmax, fabs( fx[2*k] ) );
    }
    printsigparams();

    // scale floating point values to be less than 32767 in
    // absolute value and load integer x[k] array.

    fscale = max( 1.0, fxmax/32767.0 );
    for (k = 0; k < 2*fsize; k++)
        x[k] = (int) ( fx[k] / fscale + 0.5);  // 0.5 rounds

    if (fscale > 1.0 ) printf("input divided by %7.5f\n"\
                              "to fit in 16-bit integer\n"\
                              "------------------------\n", fscale);
    dumpx( "before fftcplx()" );

	// compute the FFT

    timer = MS_TIMER;
    fftcplx(x, fsize, &blockexp);
    timer = MS_TIMER-timer;
    printf("%4d-point fftcplx() took %ld ms\n\n", fsize, timer);
    dumpx( "after fftcplx()" );

	// compute the inverse FFT

    fftcplxinv(x, fsize, &blockexp);
	dumpx("after fftcplxinv()");

} // end of main

//-------------------- Utility Functions ------------------------

void dumpx( char *title ) {
    int k;
    printf("x[n] %s\n", title);
    for (k = 0; k < fsize; k++)
        printf("%4d  %6d %6d\n", k, x[2*k], x[2*k + 1]);
    printf("blockexp is %d\n\n", blockexp);
}

void printsigparams( void ) {
    printf("---------- Signal Parameters ----------\n");
    printf("%7.f (dc)\n", DC);
    printf("%7.f sin(%6.4f * twoPI * k + %4.2f)\n", A1, F1, P1);
    printf("%7.f sin(%6.4f * twoPI * k + %4.2f)\n", A2, F2, P2);
    printf("%7.f sin(%6.4f * twoPI * k + %4.2f)\n", A3, F3, P3);
    printf("---------------------------------------\n");
}
