/************************************************************************
  fftbench
 ------------------------------------------------------------------------
  Function:
      The FFT routine under test is characterized in terms of speed and
      accuracy. The routine is tested at different FFT lengths and the
      test results are saved in the file FFTBENCH.TXT.
      
      The FFT benchmark program is distributed with the mixed radix FFT
      routine. If you want to test another FFT routine, you probably anly
      have to modify the make file, the FFT routine declaration and the 
      fftio routine below.
 ------------------------------------------------------------------------
  Compile options:
      The huge memory model is used.
 ------------------------------------------------------------------------
  Benchmarks:                   
      As an example the mixed radix FFT routine was tested and the 
      results are shown below. The Microsoft Visual C++ compiler was used 
      with the following compile options:
      /nologo /Gs /G2 /W4 /AH /Ox /D "NDEBUG" /D "_DOS" /FR
      and the test executed on a 50MHz 486DX :
      
      Length  Time [s]  Accuracy [dB]

         128   0.0054     -314.8   
         256   0.0116     -309.8   
         512   0.0251     -290.8   
        1024   0.0567     -313.6   
        2048   0.1203     -306.4   
        4096   0.2600     -291.8   
        8192   0.5800     -305.1   
         100   0.0040     -278.5   
         200   0.0099     -280.3   
         500   0.0256     -278.5   
        1000   0.0540     -278.5   
        2000   0.1294     -280.6   
        5000   0.3300     -278.4   
       10000   0.7133     -278.5   
 ------------------------------------------------------------------------
  Author:
      Jens Joergen Nielsen            For non-commercial use only.
      Bakkehusene 54                  A $100 fee must be paid if used
      DK-2970 Hoersholm               commercially. Please contact.
      DENMARK

      E-mail : jjn@get2net.dk   All rights reserved. October 2000.
      Homepage : http://home.get2net.dk/jjn
*************************************************************************/

#include <math.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>
#include <sys\timeb.h>

#define  maxIndex  10000L


/* Declaration of FFT routine under test.
 */
void fft(int n,double are[],double aim[],double bre[],double bim[]);

/* The purpose of this routine, which is called from the benchmark  
 * program, is to convert between the data formats used in the      
 * benchmark routine and the formats used in the FFT routine under  
 * test.                                                            
 */
int fftio(long n,double are[],double aim[],double bre[],double bim[])
{
    fft((int)n, are, aim, bre, bim);
    return(0);
}

/* Calculate time in seconds from start. Resolution 1/1000 second.
 */
double timeUsed(struct _timeb *start)
{                               
    double seconds;
    struct _timeb   finish;

    _ftime(&finish);
    seconds = (double) finish.millitm;
    seconds = seconds - start->millitm;
    seconds = seconds / 1000;
    seconds = seconds + finish.time - start->time;
    
    return(seconds);
}


/* Measure the time taken to perform an FFT. The test signal is a 
 * linear sweep from 0 to fs/2. At least 2 seconds of test time is
 * used and output is the average time for one FFT in seconds.
 */
double fftTime(long fftLength, 
               double xre[], double xim[], double yre[], double yim[])
{
    long            i, k, n, m, q;  
    double          w, pi;
    double          elapsed_time;
    struct _timeb   start;
                             
    n = fftLength;                             
    pi = 4 * atan(1); 
    k  = 1;
    w  = k*pi/n/n;
    for (i=0; i<n; i++) {
        xre[i] = cos(i*i*w);
        xim[i] = sin(i*i*w);
    }                    
    
    m = 0;
    q = 1;            
    _ftime(&start);
    while (timeUsed(&start) < 2.0)
    {
        for (i=0; i<q; i++)
            fftio(n, xre, xim, yre, yim);
        m = m + q;
        q = q + q;
    }
    elapsed_time = timeUsed(&start)/(double)m;
    return(elapsed_time);
}

/* Measure the difference between the theoretical FFT output and the
 * actual FFT output. The power leak is measured in terms of error 
 * to signal power ratio. The value is returned in dB.
 */
double fftLeak(long fftLength, 
               double xre[], double xim[], double yre[], double yim[])
{
    long            i, k, n;  
    double          w, pi;
    double          sqr, sqrsum, sqrmax;
                             
    n = fftLength;                             
    pi = 4 * atan(1); 
    k = 3;
    w = 2*k*pi/n;
    for (i=0; i<n; i++) {
        xre[i] = cos(i*w);
        xim[i] = sin(i*w);
    }                    
                     
    fftio(n, xre, xim, yre, yim);

    sqrsum = (yre[k]-n)*(yre[k]-n) + yim[k]*yim[k];
    sqrmax = 0;         
    for (i=0; i<n; i++)
        if (i != k)
        {
            sqr = yre[i]*yre[i] + yim[i]*yim[i];
            sqrsum = sqrsum + sqr;
            if (sqr > sqrmax) sqrmax = sqr;
        }
    sqrsum = sqrsum/n;    /* division by n*n in two stages to avoid */
    sqrsum = sqrsum/n;    /* integer overflow !!!!                  */
    sqrmax = sqrmax/n;
    sqrmax = sqrmax/n;
             
    if (sqrsum > 0)
        sqrsum = 10*log10(sqrsum);
    else                                
        sqrsum = -1000;
    if (sqrmax > 0)
        sqrmax = 10*log10(sqrmax);
    else
        sqrmax = -1000;
        
    return(sqrsum);
}

/* FFT benchmark main program.
 */
void main()
{                
    double          t, p;
    long            i, n;
    long            nmax, nfft[100];  

    double   *are, *aim, *bre, *bim;             
    
    FILE *txtFile;

    are = (double  *) calloc(maxIndex, sizeof(double));                        
    aim = (double  *) calloc(maxIndex, sizeof(double));                        
    bre = (double  *) calloc(maxIndex, sizeof(double));
    bim = (double  *) calloc(maxIndex, sizeof(double));

    if( bim == NULL )
    {
       printf( "Insufficient memory available\n" );
       exit(2);
    }                 
      
    printf("\nFFT Benchmark Test");        
    txtFile = fopen("fftbench.txt", "w");
    fprintf(txtFile,"\nResults from FFT benchmark test:");
    fprintf(txtFile,"\n\nLength  Time [s]  Accuracy [dB]\n");
              
    do {        
        printf("\n\n");
        printf("(0)  End FFT benchmark test.\n");
        printf("(1)  FFT lengths : 128, 256, 512, ... , 8192.\n");
        printf("(2)  FFT lengths : 100, 200, 500, ... , 10000.\n");
        printf("(3)  User selectable FFT length.\n\n");
        printf("Enter selection, 0-3 : ");
        scanf( "%ld", &n );
        nmax = 0;
        switch(n) {
            case 0 : nmax = 0; break;
            case 1 : 
                nmax = 7;
                nfft[0]  =  128;
                nfft[1]  =  256;
                nfft[2]  =  512;
                nfft[3]  = 1024;
                nfft[4]  = 2048;
                nfft[5]  = 4096;
                nfft[6]  = 8192;
                break;
            case 2 : 
                nmax = 7;
                nfft[0]  =   100;
                nfft[1]  =   200;
                nfft[2]  =   500;
                nfft[3]  =  1000;
                nfft[4]  =  2000;
                nfft[5]  =  5000;
                nfft[6]  = 10000;
                break;
            case 3 : 
                printf(     "\nEnter FFT length, N<=%5ld, maxPrime <= 37     : ", maxIndex);
                scanf( "%ld", &n );
                if ((n > 0) && (n <= maxIndex)) {
                    nfft[0] = n;
                    nmax = 1;
                }
                break;
        }
                
        for (i=0; i<nmax; i++) {
            t = fftTime( nfft[i], are, aim, bre, bim);
            p = fftLeak( nfft[i], are, aim, bre, bim);
            printf(         "\n %5ld %8.4f s %8.1f dB", nfft[i], t, p);
            fprintf(txtFile,"\n %5ld %8.4f   %8.1f   ", nfft[i], t, p);
        }                    
    } while(n > 0);

    fclose(txtFile);   

    free(bim);
    free(bre);
    free(aim);
    free(are);
    printf("\n");
    printf("\nFFT test terminated.");
    printf("\nCopyright 2000, Jens Joergen Nielsen, jjn@get2net.dk.");
    printf("\n");
    
 }