/*
   nasm -f elf wav_korr_asm.asm; gcc -O2 -s -o wav_korr wav_korr.c wav_korr_asm.o -lm
     or:
   gcc -DNONASM -O2 -s -o wav_korr wav_korr.c -lm
*/

#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <memory.h>

typedef signed short stereo [2];
typedef signed short mono;
typedef struct {
    unsigned long long  n;
    long double         x;
    long double         x2;
    long double         y;
    long double         y2;
    long double         xy;
} korr_t;

#ifdef NONASM

void analyze_stereo ( const stereo* p, size_t len, korr_t* k )
{
    long double  _x = 0, _x2 = 0, _y = 0, _y2 = 0, _xy = 0;
    double       t1;
    double       t2;

    k -> n  += len;

    for ( ; len--; p++ ) {
        _x  += (t1 = (*p)[0]); _x2 += t1 * t1;
        _y  += (t2 = (*p)[1]); _y2 += t2 * t2;
                               _xy += t1 * t2;
    }

    k -> x  += _x ;
    k -> x2 += _x2;
    k -> y  += _y ;
    k -> y2 += _y2;
    k -> xy += _xy;
}

void analyze_dstereo ( const stereo* p, size_t len, korr_t* k )
{
    static double l0 = 0;
    static double l1 = 0;
    long double   _x = 0, _x2 = 0, _y = 0, _y2 = 0, _xy = 0;
    double        t1;
    double        t2;

    k -> n  += len;

    for ( ; len--; p++ ) {
        _x  += (t1 = (*p)[0] - l0);  _x2 += t1 * t1;
        _y  += (t2 = (*p)[1] - l1);  _y2 += t2 * t2;
                                     _xy += t1 * t2;
        l0   = (*p)[0];
        l1   = (*p)[1];
    }

    k -> x  += _x ;
    k -> x2 += _x2;
    k -> y  += _y ;
    k -> y2 += _y2;
    k -> xy += _xy;
}


void analyze_mono   ( const mono* p, size_t len, korr_t* k )
{
    long double   _x = 0, _x2 = 0;
    double        t1;

    k -> n  += len;

    for ( ; len--; p++ ) {
        _x  += (t1 = (*p)); _x2 += t1 * t1;
    }

    k -> x  += _x ;
    k -> x2 += _x2;
    k -> y  += _x ;
    k -> y2 += _x2;
    k -> xy += _x2;
}

void analyze_dmono   ( const mono* p, size_t len, korr_t* k )
{
    static double l0 = 0;
    long double   _x = 0, _x2 = 0;
    double        t1;

    k -> n  += len;

    for ( ; len--; p++ ) {
        _x  += (t1 = (*p) - l0); _x2 += t1 * t1;
        l0   = *p;
    }

    k -> x  += _x ;
    k -> x2 += _x2;
    k -> y  += _x ;
    k -> y2 += _x2;
    k -> xy += _x2;
}

#else

extern void __analyze_stereo ( const stereo* p, size_t len, korr_t* dst );
extern void __analyze_mono   ( const mono*   p, size_t len, korr_t* dst );

#define analyze_stereo(ptr,len,k)       __analyze_stereo ( ptr, len, k )
#define analyze_mono(ptr,len,k)         __analyze_mono   ( ptr, len, k )

#endif


void report_init ( void )
{
    printf ( " x[AC]     y[AC]     r         type        x[DC]     y[DC]   sy/sx  File\n");
}

int sgn ( long double x )
{
    if ( x > 0 ) return +1;
    if ( x < 0 ) return -1;
    return 0;
}


void report ( int channels, korr_t* k, int supress_DC_report )
{
    long double  scale = sqrt ( 1.e5 / (1<<29) ); // Sine Full Scale is +10 dB, 7327 = 100%
    long double  r;
    long double  sx;
    long double  sy;
    long double  x;
    long double  y;
    long double  b;

    // printf ("n=%Lu (7036596)\nx=%Lf (-1136345)\ny=%Lf (-783749)\nxy=%Lf (103252784558988)\nx²=%Lf (182857029624921)\ny²=%Lf (201045032621475)\n\n",k.n,k.x,k.y,k.xy,k.x2,k.y2);

    r  = (k->x2*k->n - k->x*k->x) * (k->y2*k->n - k->y*k->y);
    r  = r  > 0.l  ?  (k->xy*k->n - k->x*k->y) / sqrt (r)  :  1.l;
    sx = k->n > 1  ?  sqrt ( (k->x2 - k->x*k->x/k->n) / (k->n - 1) )  :  0.l;
    sy = k->n > 1  ?  sqrt ( (k->y2 - k->y*k->y/k->n) / (k->n - 1) )  :  0.l;
    x  = k->n > 0  ?  k->x/k->n  :  0.l;
    y  = k->n > 0  ?  k->y/k->n  :  0.l;

    b  = sx != 0   ?  sy/sx * sgn(r)  :  0.l;

    printf ( "%7.3Lf%%%9.3Lf%%%9.3Lf%%   ",
             sx * scale, sy * scale, 100. * r );
    if ( supress_DC_report )
        printf ( "                  " );
    else
        printf ( "%7.3Lf%%%9.3Lf%%%", x * scale, y * scale );
    printf ( "  %6.3Lf\n", b );
    fflush ( stdout );
}


void readfile ( const char* name, int fd )
{
    unsigned short  header [22];
    stereo          s [1152];
    mono            m [1152];
    size_t          samples;
    korr_t          k0;
    korr_t          k1;
    korr_t          kd;


    memset ( &k0, 0, sizeof(k0) );
    memset ( &k1, 0, sizeof(k1) );

    read ( fd, header, sizeof(header) );

    switch ( header[11] ) {
    case 1:
        printf ("\n%s\n", name);
        while  ( ( samples = read (fd, m, sizeof(m)) ) > 0 ) {
            analyze_mono    ( m, samples / sizeof (*m), &k0 );
            analyze_dmono   ( m, samples / sizeof (*m), &k1 );
        }
        report ( header[11], &k0, 0 );
        report ( header[11], &k1, 1 );
        break;

    case 2:
        printf ("\n%s\n", name);
        while  ( ( samples = read (fd, s, sizeof(s)) ) > 0 ) {
            analyze_stereo  ( s, samples / sizeof (*s), &k0 );
            analyze_dstereo ( s, samples / sizeof (*s), &k1 );
            memset ( &kd, 0, sizeof(kd) );
            analyze_dstereo ( s, samples / sizeof (*s), &kd );
        }
        report ( header[11], &k0, 0 );
        report ( header[11], &k1, 1 );
        break;

    default:
        fprintf ( stderr, "%u Channels not supported: %s\n", header[11], name );
        break;
    }
}

int main ( int argc, char** argv )
{
    char*  name;
    int    fd;

    report_init ();

    if (argc < 2)
        readfile ( "<stdin>", 0 );
    else
        while ( (name = *++argv) != NULL ) {
            if ( (fd = open ( name, O_RDONLY )) >= 0 ) {
                readfile ( name, fd );
                close ( fd );
            } else {
                fprintf ( stderr, "Can't open: %s\n", name );
            }
        }

    return 0;
}