/*
 *  Generate graph with palette
 */

#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <memory.h>
#include <math.h>

#define FS          44100

#ifndef MAX
# define MAX        2048                // max. elements
#endif
#ifndef M_PI
# define M_PI       3.1415926535897932384626433832795029L
#endif


#define FFT_NORM        0                   // y(t)     -> y(f)
#define FFT_INVS        1                   // y(f)     -> y(t)

#define FFT_ERR_OK      0                   // no error
#define FFT_ERR_LD      1                   // len is not a power of 2
#define FFT_ERR_MAX     2                   // len too large


typedef float   compl   [2];
compl           root    [MAX >> 1];                 // Sine-/cosine-table
size_t          shuffle [MAX >> 1] [2];             // Shuffle-table
size_t          shuffle_len;


static long double
sinus ( long double x )
{
    x -= floor (x);

    switch ( (int)(8 * x) ) {
    case 0: return +sin (2*M_PI*      x );
    case 1: return +cos (2*M_PI*(0.25-x));
    case 2: return +cos (2*M_PI*(x-0.25));
    case 3: return +sin (2*M_PI*(0.50-x));
    case 4: return -sin (2*M_PI*(x-0.50));
    case 5: return -cos (2*M_PI*(0.75-x));
    case 6: return -cos (2*M_PI*(x-0.75));
    case 7: return -sin (2*M_PI*(1.00-x));
    }
}

static long double
cosinus ( long double x )
{
    x -= floor (x);

    switch ( (int)(8 * x) ) {
    case 0: return +cos (2*M_PI*      x );
    case 1: return +sin (2*M_PI*(0.25-x));
    case 2: return -sin (2*M_PI*(x-0.25));
    case 3: return -cos (2*M_PI*(0.50-x));
    case 4: return -cos (2*M_PI*(x-0.50));
    case 5: return -sin (2*M_PI*(0.75-x));
    case 6: return +sin (2*M_PI*(x-0.75));
    case 7: return +cos (2*M_PI*(1.00-x));
    }
}

// Bitinversion

static size_t
swap ( size_t number, int bits )
{
    size_t  ret;
    for ( ret = 0; bits--; number >>= 1 ) {
        ret = ret + ret + (number & 1);
    }
    return ret;
}


// Determine the logarithmus dualis

static int
ld ( size_t number )
{
    int i;

    for ( i = 0; i < (int) (sizeof(size_t) * CHAR_BIT); i++ )
        if ( ((size_t)1 << i) == number )
            return i;

    return -1;
}

// The actual FFT

int
fft ( compl* fn, const size_t newlen, const int direction )
{
    static size_t    len  = 0;
    static int       bits = 0;
    static int       last = 0;
    register size_t  i;
    register size_t  j;
    register size_t  k;
    float*           p;
    size_t           pp;

    /* initialize tables */

    if ( newlen != len ) {
        len  = newlen;

        last = FFT_INVS;

        if ( (bits = ld(len)) == -1 )
            return FFT_ERR_LD;

        for ( i = 0; i < len; i++ ) {
            j = swap ( i, bits );
            if ( i < j ) {
                shuffle [shuffle_len] [0] = i;
                shuffle [shuffle_len] [1] = j;
                shuffle_len++;
            }
        }

        for ( i = 0; i < (len>>1); i++ ) {
            long double  x = (long double) swap ( i+i, bits ) / len;
            root [i] [0] = cosinus (x);
            root [i] [1] = sinus   (x);
        }
    }

    if ( last != (direction & FFT_INVS) ) {
        last = direction & FFT_INVS;

        p = (float*) (&root[0][1]);

        for ( i = len>>1; i--; p += 2 )
            *p = -*p;
    }

    /* Actual transformation */

    pp = len >> 1;
    do {
        float*  bp = (float*) root;
        float*  si = (float*) fn;
        float*  di = (float*) fn+pp+pp;

        do {
            k = pp;
            do {
                float  mulr = bp[0]*di[0] - bp[1]*di[1];
                float  muli = bp[1]*di[0] + bp[0]*di[1];
                float  addr = si[0];
                float  addi = si[1];

                si [0] = addr + mulr;
                si [1] = addi + muli;
                di [0] = addr - mulr;
                di [1] = addi - muli;

                si += 2, di += 2;
            } while ( --k );
            si += pp+pp, di += pp+pp, bp += 2;
        } while ( si < &fn[len][0] );
    } while ( pp >>= 1 );

    /* Bitinversion */

    for ( k = 0; k < shuffle_len; k++ ) {
        float  tmp;
        i   = shuffle [k] [0];
        j   = shuffle [k] [1];
        tmp = fn [i][0]; fn [i][0] = fn [j][0]; fn [j][0] = tmp;
        tmp = fn [i][1]; fn [i][1] = fn [j][1]; fn [j][1] = tmp;
    }

    return FFT_ERR_OK;
}


#define S             8
#if   FS == 44100
# define USED_BANDS  203
#elif FS == 48000
# define USED_BANDS  205
#else
# error
#endif

float  bands [] = {
       0,  100,  200,   300,   400,   510,   630,  770,  920, 1080,
    1270, 1480, 1720,  2000,  2320,  2700,  3150, 3700, 4400, 5300,
    6400, 7700, 9500, 12000, 15500, 20500, 27000
};


void
procedure ( float* _A, FILE* fp )
{
    static int    init = 0;
    static float  tab [USED_BANDS] [1024];
    float         A [2048];
    float         B [2048] [2];
    float         C [2048] [2];
    float         D [2048];
    int           i;
    unsigned int  band;
    double        Sum;
    double        Diff;

    if ( init == 0 ) {  // +/- 0.5 bark = -6 dB, +/- 1 bark = -96 dB

        // Center the bands differently, so that the 0. Band starts at the bottom?
        for ( band = 0; band < USED_BANDS; band++ ) {
            double  f1  = ( (S - band%S)*bands [band/S+0] + (band%S)*bands [band/S+1] ) * (2048. / FS / S);
            double  f2  = ( (S - band%S)*bands [band/S+1] + (band%S)*bands [band/S+2] ) * (2048. / FS / S);
            double  tmp = 2. / ( f2 - f1 );

            for ( i = 0; i < 1024; i++ ) {
                double  w    = (i - f1) * tmp - 1;
                double  mult = exp ( -0.69314718055994530941723212145818 * w * w * w * w);
                tab [band] [i] = mult;
            }
        }
        init = 1;
    }


    for ( i = 0; i < 2048; i++ ) {
        A [i]     = _A [i];
        B [i] [0] = A[i];
        B [i] [1] = 0.;
    }

    fft ( B, 2048, FFT_NORM );

    for ( band = 0; band < USED_BANDS; band++ ) {

        memset ( C [1024], 0, 1024 * sizeof(*C) );
        for ( i = 0; i < 1024; i++ ) {
            C [i] [0] = B [i] [0] * tab [band] [i];
            C [i] [1] = B [i] [1] * tab [band] [i];
        }
        fft ( C, 2048, FFT_INVS );

        for ( i = 0; i < 2048; i++ ) {
            D [i] = sqrt (C[i][0] * C[i][0] + C[i][1] * C[i][1]);
        }

        Sum = 1.e-70;
        for ( i = 0; i < 2048; i++ ) {
            Sum += D[i];
        }
        Sum /= 2048.;

        Diff = 1.e-70;
        for ( i = 0; i < 2048; i++ ) {
            D[i]  = D[i] / Sum - 1;
            Diff += D[i] * D[i];
        }
        Diff = 1 - 2 * sqrt ( Diff / 2048.);        // 1 for sine signals, ~0 for noise

        printf ("%4.0f", Diff*100 );

        Diff = 170 * Diff + 85;
        if ( Diff >= 0 )
            putc ( (int)Diff, fp );
        else
            putc ( 0, fp );
    }

    printf ("\n");
    return;
}


int
main ( int argc, char** argv )
{
    static char          buff [1 << 18];
    static float         A [4000000];
    static signed short  b [35 * FS] [2];
    float*               p = A;
    int                  i;
    int                  len;
    long double          w = 0.L;
    FILE*                graph;
    FILE*                fp;
    char                 name [256];

    freopen ( "report", "w", stdout );
    setvbuf ( stdout, buff, _IOFBF, sizeof buff );

    // Calculate tone
    while ( *++argv ) {

        p = A;

        if ( 0 ) {
            for ( i = 0; i < 10000; i++ )
                *p++ = (rand () + rand ()) / (double) RAND_MAX - 1;

            for ( i = 0; i < 200000; i++ ) {
                w   +=  i / 200000. * M_PI;
                *p++ = sin (w);
            }
            for ( i = 0; i < 20000; i++ ) {
                w   +=  i / 20000. * M_PI;
                *p++ = sin (w);
            }
            for ( i = 0; i < 2000; i++ ) {
                w   +=  i / 2000. * M_PI;
                *p++ = sin (w);
            }

            for ( i = 0; i < 200000; i++ ) {
                w   +=  i / 200000. * M_PI;
                *p++ = sin (w) + (rand () + rand ()) / (double) RAND_MAX - 1;
            }
            for ( i = 0; i < 20000; i++ ) {
                w   +=  i / 20000. * M_PI;
                *p++ = sin (w) + (rand () + rand ()) / (double) RAND_MAX - 1;
            }
            for ( i = 0; i < 2000; i++ ) {
                w   +=  i / 2000. * M_PI;
                *p++ = sin (w) + (rand () + rand ()) / (double) RAND_MAX - 1;
            }
        }
        else {

            fp = fopen ( *argv, "rb" );
            fread ( b, 1, 44, fp );
            len = fread ( b, sizeof(*b), sizeof(b)/sizeof(*b), fp );
            fclose (fp);
            fprintf ( stderr, "Add WAV file: %u\n", len );
            for ( i = 0; i < len; i++ )
                *p++ = b[i][0] * 1.e-4;
        }

        fprintf ( stderr, "Total length: %u\n", p - A );
        sprintf ( name, "%s.ppm", *argv );
        graph = fopen ( name, "wb" );
        fprintf ( graph, "P5\n%u %u\n255\n", USED_BANDS, (p - A - MAX) / 256 );

        // Analyze tone
        for ( i = 0; i < p - A - MAX; i += 256 ) {
            fprintf ( stderr, "Proceed: %7u %3.0f%%\r", i, 100.*i/(p - A - MAX) );
            printf ( "%7u: ", i );
            procedure ( A + i, graph );
        }

        fclose (graph);

    }

    return 0;
}