/******************************************************************************
*                             SoundTrack, v. 3.20                             *
*                                                                             *
*                                     by                                      *
*                                                                             *
*             L. Padilla (e-mail: padilla at domain "gae ucm es")             *
*                                                                             *
*                        Madrid, July 2003 - June 2012                        *
*                                                                             *
* Compile with: cc soundtrack.c -o soundtrack -lm                             *
* Latest version at http://www.gae.ucm.es/~padilla/extrawork/soundtrack.c     *
******************************************************************************/


#include   <math.h>
#include  <stdio.h>
#include <stdlib.h>
#include <string.h>


#define TRUE   1
#define FALSE  0

#define BITBCD    5  /*  Bits in a  BCD   byte  */
#define BITALPHA  7  /*  Bits in an ALPHA byte  */

#define SSBCD    0x0B  /*  Start sentinel in ANSI BCD    */
#define SSALPHA  0x05  /*  Start sentinel in ANSI ALPHA  */
#define ESBCD    0x0F  /*  End   sentinel in ANSI BCD    */
#define ESALPHA  0x1F  /*  End   sentinel in ANSI ALPHA  */

#define OFFBCD    48  /*  Start of BCD   characters in ASCII table  */
#define OFFALPHA  32  /*  Start of ALPHA characters in ASCII table  */

#define CLKLOW   20  /*  Max. clocking bits to use in low  density format  */
#define CLKHIGH  60  /*  Max. clocking bits to use in high density format  */
#define THDLOW    8  /*  Minimum points to be a good peak in low  density  */
#define THDHIGH   4  /*  Minimum points to be a good peak in high density  */

#define NUM  2048  /*  Maximum number of bits to record         */
#define NUMS  500  /*  Points under threshold to stop read out  */

#define SLOPE  0  /*  Detection methods  */
#define PEAK   1
#define ZERO   2

#define ABS(x)    (((double) (x) < 0.0) ? (-(x)) : (x))
#define ROUND(x)  ((int) ((float) (x) + 0.5))
#define SQR(x)    ((x) * (x))

#define STATPRINT() \
        mean /= (double) count; \
        sigma = sqrt (sigma/(double) count - SQR(mean)); \
        printf ("Swipe time: %.3f s\n", t - t0); \
        printf ("Data max./mean/rms amplitude: %.3f/%.3lf/%.3lf\n", max, mean, \
                sigma); \
        printf ("Bits: %i, individual duration: %.3f ms\n", i, 1000. * \
                (t - t0)/(float) i);


main (int argc, char * argv[])
{
    int i, j, track = 2, printdata = FALSE, printbits = FALSE, count = 0;
    int seq = 0, num = 0, nums = -1, num0 = 0, thd = 0, nostrict = TRUE;
    int invclk = FALSE, data = FALSE, seen1 = FALSE, tclk = FALSE, bits[NUM];
    int bit, off, ss, es, clk = -1, parity, CheckLRC, LRCBit[BITALPHA];
    int byte[NUM/BITBCD], LRC, nostop = FALSE, p = 0, method = -1;
    float t, s, n, max = 0., threshold = -1., t0, ot, os;
    float m = 1.e-10, om = 0., pm = 0., pt, hpt, dur, hdur, dur0, dur1;
    double mean = 0., sigma = 0.;

    printf ("\n");


/*  Processing command line arguments  */

    for (i = 1; i < argc; i++)
    {
        if (strcmp (argv[i], "-d") == 0)
        {
            printdata = TRUE;
        }
        else if (strcmp (argv[i], "-b") == 0)
        {
            printbits = TRUE;
        }
        else if (strcmp (argv[i], "-n") == 0)
        {
            nostop = TRUE;
        }
        else if (strcmp (argv[i], "-i") == 0)
        {
            invclk = TRUE;
        }
        else if (strcmp (argv[i], "-x") == 0)
        {
            nostrict = FALSE;
        }
        else if (strcmp (argv[i], "-t") == 0)
        {
            i++;
            if (i < argc)
            {
                track = atoi (argv[i]);
            }
            else
            {
                printf ("Parameter -t error!\n\n");
                return -1;
            }
        }
        else if (strcmp (argv[i], "-c") == 0)
        {
            i++;
            if (i < argc)
            {
                clk = atoi (argv[i]);
            }
            else
            {
                printf ("Parameter -c error!\n\n");
                return -1;
            }
        }
        else if (strcmp (argv[i], "-p") == 0)
        {
            i++;
            if (i < argc)
            {
                thd = atoi (argv[i]);
            }
            else
            {
                printf ("Parameter -p error!\n\n");
                return -1;
            }
        }
        else if (strcmp (argv[i], "-s") == 0)
        {
            i++;
            if (i < argc)
            {
                nums = atoi (argv[i]);
            }
            else
            {
                printf ("Parameter -s error!\n\n");
                return -1;
            }
        }
        else if (strcmp (argv[i], "-l") == 0)
        {
            i++;
            if (i < argc)
            {
                threshold = (float) atof (argv[i]);
            }
            else
            {
                printf ("Parameter -l error!\n\n");
                return -1;
            }
        }
        else if (strcmp (argv[i], "-m") == 0)
        {
            i++;
            if (i < argc)
            {
                if (strcmp (argv[i], "s") == 0)
                    method = SLOPE;
                else if (strcmp (argv[i], "p") == 0)
                    method = PEAK;
                else if (strcmp (argv[i], "z") == 0)
                    method = ZERO;
                else
                {
                    printf ("Parameter -m error!\n\n");
                    return -1;
                }
            }
            else
            {
                printf ("Parameter -m error!\n\n");
                return -1;
            }
        }
        else
        {
	    printf ("Usage: soundtrack [-d] [-b] [-n] [-i] [-x] ");
            printf ("[-t <track>] [-c <clk_bits>]\n");
            printf ("                  [-p <peak_points>] [-s <stop_points>]");
            printf (" [-l <threshold>]\n");
            printf ("                  [-m <s|p|z>]\n\n");
            printf ("Options:\n\n");
            printf ("-d: Print raw data to stdout.\n");
            printf ("-b: Print bits to stdout.\n");
            printf ("-n: No stop after LRC.\n");
            printf ("-i: Inverse clocking bits, 1s instead of 0s.\n");
            printf ("-x: Strict detection, all peaks above threshold.\n\n");
            printf ("-t <track>      : Track number.\n");
            printf ("-c <clk_bits>   : Number of clocking bits.\n");
            printf ("-p <peak_points>: Minimum peak duration.\n");
            printf ("-s <stop_points>: Pause in data before stopping.\n");
            printf ("-l <threshold>  : Minimum level to be data.\n");
            printf ("-m <s|p|z>      : Detection method:\n");
            printf ("                      s = slope\n");
            printf ("                      p = peak (default)\n");
            printf ("                      z = zero crossing\n\n");


            return 2;
        }
    }

    if (track == 1)  /*  Track 1  */
    {
        bit = BITALPHA;
        off = OFFALPHA;
        ss  = SSALPHA;
        es  = ESALPHA;
        if (clk < 0)
        {
            clk = CLKHIGH;
            printf ("Bad or no clocking bits parameter, using default: %i\n",
                    clk);
        }
        if (thd < 1)
        {
            if (method == ZERO)
                thd = 1;
            else
                thd = THDHIGH;

            printf ("Bad or no peak duration parameter, using default: %i\n",
                    thd);
        }
    }
    else if (track == 2)  /*  Track 2  */
    {
        bit = BITBCD;
        off = OFFBCD;
        ss  = SSBCD;
        es  = ESBCD;
        if (clk < 0)
        {
            clk = CLKLOW;
            printf ("Bad or no clocking bits parameter, using default: %i\n",
                    clk);
        }
        if (thd < 1)
        {
            if (method == ZERO)
                thd = 1;
            else
                thd = THDLOW;

            printf ("Bad or no peak duration parameter, using default: %i\n",
                    thd);
        }
    }
    else if (track == 3)  /*  Track 3  */
    {
        bit = BITBCD;
        off = OFFBCD;
        ss  = SSBCD;
        es  = ESBCD;
        if (clk < 0)
        {
            clk = CLKHIGH;
            printf ("Bad or no clocking bits parameter, using default: %i\n",
                    clk);
        }
        if (thd < 1)
        {
            if (method == ZERO)
                thd = 1;
            else
                thd = THDHIGH;

            printf ("Bad or no peak duration parameter, using default: %i\n",
                    thd);
        }
    }
    else if (track == 4)  /*  For non standard tracks  */
    {
        bit = BITBCD;
        off = OFFBCD;
        ss  = ESBCD + 1;  /*  Non existent  */
        es  = ESBCD + 1;  /*  Non existent  */
        if (clk < 0)
        {
            clk = CLKLOW;
            printf ("Bad or no clocking bits parameter, using default: %i\n",
                    clk);
        }
        if (thd < 1)
        {
            if (method == ZERO)
                thd = 1;
            else
                thd = THDHIGH;

            printf ("Bad or no peak duration parameter, using default: %i\n",
                    thd);
        }
    }
    else
    {
        printf ("Bad track number, exiting.\n\n");
        return 3;
    }

    if (nums < 0)
    {
        nums = NUMS;
        printf ("Bad or no stop duration parameter, using default: %i\n",
                nums);
    }

    if (method < 0)
    {
        method = PEAK;
        printf ("No method parameter, using default: peak\n");
    }

    printf ("Reading track %i\n", track);


/*  Input data read-out  */

/*  Skip some inconvenient initial data  */

    scanf ("; %*s %*s %*i\n");
    for (i = 0; i < 500 && scanf ("%f %f\n", & t, & s) == 2; i++)
        ;

/*  Initial data gives noise level  */

    if (threshold < 0.)
    {
        printf ("Bad or no threshold level parameter, using data...\n");

        for (i = 0; i < 5000 && scanf ("%f %f\n", & t, & s) == 2; i++)
        {
            if (ABS(s) > max)
                max = ABS(s);
            mean  += (double) s;
            sigma += (double) SQR(s);
            count++;
        }

        mean /= (double) count;
        sigma = sqrt (sigma/(double) count - SQR(mean));
        threshold = 10. * (float) sigma;
        printf ("Noise max./mean/rms amplitude: %.3f/%.3lf/%.3lf\n", max, mean,
                sigma);
        printf ("Threshold: %.3f\n", threshold);
    }

    n = (float) mean;

/*  Skipping noise  */

    while (scanf ("%f %f\n", & t, & s) == 2 && ABS(s - n) < threshold)
    {
        ot = t;
        pt = t;
        os = s - n;
    }

/*  Real data detected  */

    t0 = t;
    i = 0;
    max   = 0.;
    mean  = 0.;
    sigma = 0.;
    count = 0;

    while (scanf ("%f %f\n", & t, & s) == 2)
    {
        s = s - n;

        if (printdata)  /*  Printing raw data if option given  */
            printf ("%f %f\n", t, s);

        if (ABS(s) > max)
            max = ABS(s);
        mean  += (double) s;
        sigma += (double) SQR(s);
        count++;

        if (ABS(s) < threshold)  /*  Detecting end of data  */
            num++;
        else
            num = 0;

        if (num > nums)  /*  Stop if end of data  */
        {
            STATPRINT()
            break;
        }


        if (ABS(s) < 1.e-10) s = 1.e-10;

        if (method == SLOPE)
        {
            /*  Determining slope  */

            m = s - os;
            if (m == 0.)
                m = om;
            else
                m /= ABS(m);

            if (m != pm && seq == 0 && (nostrict || ABS(s) > threshold))
                                                          /*  Peak detected  */
            {
                seq = 1;
            }
            else if (m == om && seq > 0)  /*  Peak looks good  */
            {
                seq += 1;

                if (seq >= thd)  /*  Peak confirmed  */
                {
                    if (p == 0)  /*  First peak  */
                    {
                        hdur = t - ot;
                        p++;
                    }
                    else  /*  Second peak  */
                    {
                        dur = hdur + t - ot;
                        p = 0;
                        data = TRUE;
                    }
                    ot = t;
                    seq = 0;
                    pm = m;
                }
            }
            else  /*  No good peak yet  */
            {
                seq = 0;
            }
            os = s;
            om = m;
        }
        else if (method == PEAK)
        {
            if (ABS(s) > ABS(m))
            {
                m = s;
                pt = t;
            }

            seq += 1;

            if (s/ABS(s) != m/ABS(m) && seq >= thd && (nostrict || ABS(m) >
                                             threshold))  /*  Peak detected  */
            {
                if (p == 0)  /*  First peak  */
                {
                    hdur = pt - ot;
                    p++;
                }
                else  /*  Second peak  */
                {
                    dur = hdur + pt - ot;
                    p = 0;
                    data = TRUE;
                }
                m = 1.e-10;
                ot = pt;
                seq = 0;
            }
        }
        else if (method == ZERO)
        {
            if (ABS(os) < 1.e-10) os = 1.e-10;

            seq += 1;

            if (s/ABS(s) != os/ABS(os) && seq >= thd)  /*  Zero crossing  */
            {
                if (p == 0)  /*  First peak  */
                {
                    hdur = t - (t - ot) * s/(s - os) - pt;
                    hpt = hdur + pt;
                    p++;
                }
                else  /*  Second peak  */
                {
                    dur = t - (t - ot) * s/(s - os) - pt;
                    pt = dur + pt;
                    p = 0;
                    data = TRUE;
                }
                seq = 0;
            }
            ot = t;
            os = s;
        }

        if (data)
        {
            data = FALSE;

            /*  Calculating bit durations  */

            if (i == 0)  /*  Discard first two peaks  */
            {
                if (! invclk)
                {
                    dur0 = 0.0007;
                    dur1 = dur0 * 0.60;
                    bits[i] = 0;
                    i += 1;
                    bits[i] = 0;
                }
                else
                {
                    dur1 = 0.0004;
                    dur0 = dur1 / 0.60;
                    bits[i] = 1;
                    i += 1;
                    bits[i] = 1;
                }
                i += 1;
            }
            else if (i <= 2 && i < clk)  /*  First two clocking bits  */
            {
                if (! invclk)
                {
                    dur0 = dur/2.;
                    dur1 = dur0 * 0.60;
                    bits[i] = 0;
                    i += 1;
                    bits[i] = 0;
                }
                else
                {
                    dur1 = dur/2.;
                    dur0 = dur1 / 0.60;
                    bits[i] = 1;
                    i += 1;
                    bits[i] = 1;
                }
                i += 1;
            }
            else if (i > 1 && i < clk)  /*  Next clocking bits  */
            {
                if (! invclk)
                {
                    dur0 = (dur0 + dur/2.)/2.;
                    dur1 = dur0 * 0.60;
                    bits[i] = 0;
                    i += 1;
                    bits[i] = 0;
                }
                else
                {
                    dur1 = (dur1 + dur/2.)/2.;
                    dur0 = dur1 / 0.60;
                    bits[i] = 1;
                    i += 1;
                    bits[i] = 1;
                }
                i += 1;
            }
            else  /*  Data bits  */
            {
                if (dur > (3. * dur0 + dur1)/2.)  /*  Bits = 00  */
                {
                    /*  Update duration if good value  */

                    if (ABS(dur/2. - dur0) < dur0 - dur1)

                        dur0 = (dur0 + dur/2.)/2.;

                    if (invclk)
                    {
                        seen1 = TRUE;
                        num0 = 0;
                    }
                    bits[i] = 0;
                    i += 1;
                    bits[i] = 0;
                    i += 1;

                    /*  Detection of trailing clocking bits  */

                    if (seen1 == TRUE && ! invclk)
                    {
                        num0 += 1;
                        if (num0 > 12) tclk = TRUE;
                    }
                }
                else if (dur < (dur0 + 3. * dur1)/2.)  /*  Bits = 11  */
                {
                    /*  Update duration if good value  */

                    if (ABS(dur/2. - dur1) < dur0 - dur1)

                        dur1 = (dur1 + dur/2.)/2.;

                    if (! invclk)
                    {
                        seen1 = TRUE;
                        num0 = 0;
                    }
                    bits[i] = 1;
                    i += 1;
                    if (track == 4 || (invclk && ! seen1))
                    {
                        bits[i] = 1;
                        i += 1;
                    }

                    /*  Detection of inverse trailing clocking bits  */

                    if (seen1 == TRUE && invclk)
                    {
                        num0 += 1;
                        if (num0 > 12) tclk = TRUE;
                    }
                }
                else if (track == 4 || (invclk && ! seen1))  /*  Bits = 01X  */
                                                             /*      or 10   */
                {
                    seen1 = TRUE;
                    num0 = 0;

                    if (hdur > dur - hdur)  /*  Bits = 01X  */
                    {
                        /*  Update duration if good value  */

                        if (ABS(hdur - dur0) < dur0 - dur1)

                            dur0 = (dur0 + hdur)/2.;

                        if (ABS(dur - hdur - dur1) < dur0 - dur1)

                            dur1 = (dur1 + dur - hdur)/2.;

                        bits[i] = 0;
                        i += 1;
                        bits[i] = 1;
                    }
                    else  /*  Bits = 10  */
                    {
                        /*  Update duration if good value  */

                        if (ABS(hdur - dur1) < dur0 - dur1)

                            dur1 = (dur1 + hdur)/2.;

                        if (ABS(dur - hdur - dur0) < dur0 - dur1)

                            dur0 = (dur0 + dur - hdur)/2.;

                        bits[i] = 1;
                        i += 1;
                        bits[i] = 0;
                    }

                    i += 1;

                }
                else  /*  Bits = 011  */
                {
                    if (hdur > (dur0 + dur1)/2.)
                    {
                        /*  Update duration if good value  */

                        if (ABS(hdur - dur0) < dur0 - dur1)

                            dur0 = (dur0 + hdur)/2.;

                        if (invclk)
                        {
                            seen1 = TRUE;
                            num0 = 0;
                        }
                        bits[i] = 0;
                        i += 1;

                        /*  Detection of trailing clocking bits  */

                        if (seen1 == TRUE && ! invclk)
                        {
                            num0 += 1;
                            if (num0 > 12) tclk = TRUE;
                        }

                        p = 1;
                        hdur = dur - hdur;
                        if (method == ZERO) pt = hpt;
                    }
                    else  /*  Inconsistency  */
                    {
                        printf ("Recognition error of bit %i, stopping.\n\n",
                                                                        i + 1);
                        STATPRINT()
                        break;
                    }
                }
            }
        }
    }

    num = i;


/*  Printing bits if option given  */

    if (printbits)
    {
        for (i = 0; i < num; i++)
            printf ("%i", bits[i]);

        printf ("\n");
    }


/*  Grouping bits into bytes and checking parity and LRC  */

    seq = 0;
    seen1 = FALSE;
    CheckLRC = FALSE;
    for (j = 0; j < bit; j++)
        LRCBit[j] = FALSE;

    for (i = 0; i < num; i++)  /*  Loop for all bits  */

        if ((bits[i] && ! invclk) || (! bits[i] && invclk) || seen1)
                                         /*  Skipping leading clocking bits  */
        {
            seen1 = TRUE;
            parity = FALSE;
            byte[seq] = 0;

            /*  Calculating byte and parity  */

            for (j = 0; j < (bit - 1) && i < num; j++, i++)
            {
                byte[seq] += bits[i] * ROUND(pow (2., (float) j));

                if (! CheckLRC)  /*  LRC byte does not contribute to LRC  */
                    LRCBit[j] ^= bits[i];

                parity ^= bits[i];
            }

            /*  Checking parity (odd) in standard tracks  */

            if (bits[i] == parity && track != 4)
            {
                printf ("\nParity error in byte %i, exiting.\n\n", seq + 1);
                return -2;
            }

            if (byte[seq] == ss)
            {
                printf ("Start sentinel found\n");
                printf ("%c", (char) (byte[seq] + off));
            }
            else if (byte[seq] == es)
            {
                printf ("%c", (char) (byte[seq] + off));
                printf ("\nEnd sentinel found\n");
            }
            else if (! CheckLRC)
                printf ("%c", (char) (byte[seq] + off));

            if (CheckLRC)  /*  This is LRC byte, check it  */
            {
                LRC = 0;
                for (j = 0; j < bit - 1; j++)
                    LRC += LRCBit[j] * ROUND(pow (2., (float) j));

                if (LRC == byte[seq])  /*  LRC bits are even  */
                {
                    printf ("LRC OK.\n\n");
                    seen1 = FALSE;
                    CheckLRC = FALSE;
                    for (j = 0; j < bit; j++)
                        LRCBit[j] = FALSE;
                    if (! nostop)
                        return 0;
                }
                else
                {
                    printf ("LRC error.\n\n");
                    return -2;
                }
            }

            if (byte[seq] == es)
                CheckLRC = TRUE;

            seq++;

        }

    printf ("\n");


    printf ("\n");


    return 1;
}