/* audioopmodule - Module to detect peak values in arrays */

#include "Python.h"

#if SIZEOF_INT == 4

typedef int Py_Int32;

typedef unsigned int Py_UInt32;

#else

#if SIZEOF_LONG == 4

typedef long Py_Int32;

typedef unsigned long Py_UInt32;

#else

#error "No 4-byte integral type"

#endif

#endif

typedef short PyInt16;

#if defined(__CHAR_UNSIGNED__)

#if defined(signed)

/* This module currently does not work on systems where only unsigned

   characters are available.  Take it out of Setup.  Sorry. */

#endif

#endif

/* Code shamelessly stolen from sox, 12.17.7, g711.c

** (c) Craig Reese, Joe Campbell and Jeff Poskanzer 1989 */

/* From g711.c:

 *

 * December 30, 1994:

 * Functions linear2alaw, linear2ulaw have been updated to correctly

 * convert unquantized 16 bit values.

 * Tables for direct u- to A-law and A- to u-law conversions have been

 * corrected.

 * Borge Lindberg, Center for PersonKommunikation, Aalborg University.

 * bli@cpk.auc.dk

 *

 */

#define BIAS 0x84   /* define the add-in bias for 16 bit samples */

#define CLIP 32635

#define SIGN_BIT        (0x80)          /* Sign bit for a A-law byte. */

#define QUANT_MASK      (0xf)           /* Quantization field mask. */

#define SEG_SHIFT       (4)             /* Left shift for segment number. */

#define SEG_MASK        (0x70)          /* Segment field mask. */

static PyInt16 seg_aend[8] = {0x1F, 0x3F, 0x7F, 0xFF,

                              0x1FF, 0x3FF, 0x7FF, 0xFFF};

static PyInt16 seg_uend[8] = {0x3F, 0x7F, 0xFF, 0x1FF,

                              0x3FF, 0x7FF, 0xFFF, 0x1FFF};

static PyInt16

search(PyInt16 val, PyInt16 *table, int size)

{

        int i;

        for (i = 0; i < size; i++) {

                if (val <= *table++)

                        return (i);

        }

        return (size);

}

#define st_ulaw2linear16(uc) (_st_ulaw2linear16[uc])

#define st_alaw2linear16(uc) (_st_alaw2linear16[uc])

static PyInt16 _st_ulaw2linear16[256] = {

    -32124,  -31100,  -30076,  -29052,  -28028,  -27004,  -25980,

    -24956,  -23932,  -22908,  -21884,  -20860,  -19836,  -18812,

    -17788,  -16764,  -15996,  -15484,  -14972,  -14460,  -13948,

    -13436,  -12924,  -12412,  -11900,  -11388,  -10876,  -10364,

     -9852,   -9340,   -8828,   -8316,   -7932,   -7676,   -7420,

     -7164,   -6908,   -6652,   -6396,   -6140,   -5884,   -5628,

     -5372,   -5116,   -4860,   -4604,   -4348,   -4092,   -3900,

     -3772,   -3644,   -3516,   -3388,   -3260,   -3132,   -3004,

     -2876,   -2748,   -2620,   -2492,   -2364,   -2236,   -2108,

     -1980,   -1884,   -1820,   -1756,   -1692,   -1628,   -1564,

     -1500,   -1436,   -1372,   -1308,   -1244,   -1180,   -1116,

     -1052,    -988,    -924,    -876,    -844,    -812,    -780,

      -748,    -716,    -684,    -652,    -620,    -588,    -556,

      -524,    -492,    -460,    -428,    -396,    -372,    -356,

      -340,    -324,    -308,    -292,    -276,    -260,    -244,

      -228,    -212,    -196,    -180,    -164,    -148,    -132,

      -120,    -112,    -104,     -96,     -88,     -80,     -72,

       -64,     -56,     -48,     -40,     -32,     -24,     -16,

        -8,       0,   32124,   31100,   30076,   29052,   28028,

     27004,   25980,   24956,   23932,   22908,   21884,   20860,

     19836,   18812,   17788,   16764,   15996,   15484,   14972,

     14460,   13948,   13436,   12924,   12412,   11900,   11388,

     10876,   10364,    9852,    9340,    8828,    8316,    7932,

      7676,    7420,    7164,    6908,    6652,    6396,    6140,

      5884,    5628,    5372,    5116,    4860,    4604,    4348,

      4092,    3900,    3772,    3644,    3516,    3388,    3260,

      3132,    3004,    2876,    2748,    2620,    2492,    2364,

      2236,    2108,    1980,    1884,    1820,    1756,    1692,

      1628,    1564,    1500,    1436,    1372,    1308,    1244,

      1180,    1116,    1052,     988,     924,     876,     844,

       812,     780,     748,     716,     684,     652,     620,

       588,     556,     524,     492,     460,     428,     396,

       372,     356,     340,     324,     308,     292,     276,

       260,     244,     228,     212,     196,     180,     164,

       148,     132,     120,     112,     104,      96,      88,

        80,      72,      64,      56,      48,      40,      32,

        24,      16,       8,       0

};

/*

 * linear2ulaw() accepts a 14-bit signed integer and encodes it as u-law data

 * stored in a unsigned char.  This function should only be called with

 * the data shifted such that it only contains information in the lower

 * 14-bits.

 *

 * In order to simplify the encoding process, the original linear magnitude

 * is biased by adding 33 which shifts the encoding range from (0 - 8158) to

 * (33 - 8191). The result can be seen in the following encoding table:

 *

 *      Biased Linear Input Code        Compressed Code

 *      ------------------------        ---------------

 *      00000001wxyza                   000wxyz

 *      0000001wxyzab                   001wxyz

 *      000001wxyzabc                   010wxyz

 *      00001wxyzabcd                   011wxyz

 *      0001wxyzabcde                   100wxyz

 *      001wxyzabcdef                   101wxyz

 *      01wxyzabcdefg                   110wxyz

 *      1wxyzabcdefgh                   111wxyz

 *

 * Each biased linear code has a leading 1 which identifies the segment

 * number. The value of the segment number is equal to 7 minus the number

 * of leading 0's. The quantization interval is directly available as the

 * four bits wxyz.  * The trailing bits (a - h) are ignored.

 *

 * Ordinarily the complement of the resulting code word is used for

 * transmission, and so the code word is complemented before it is returned.

 *

 * For further information see John C. Bellamy's Digital Telephony, 1982,

 * John Wiley & Sons, pps 98-111 and 472-476.

 */

static unsigned char

st_14linear2ulaw(PyInt16 pcm_val)	/* 2's complement (14-bit range) */

{

        PyInt16         mask;

        PyInt16         seg;

        unsigned char   uval;

        /* The original sox code does this in the calling function, not here */

        pcm_val = pcm_val >> 2;

        /* u-law inverts all bits */

        /* Get the sign and the magnitude of the value. */

        if (pcm_val < 0) {

                pcm_val = -pcm_val;

                mask = 0x7F;

        } else {

                mask = 0xFF;

        }

        if ( pcm_val > CLIP ) pcm_val = CLIP;           /* clip the magnitude */

        pcm_val += (BIAS >> 2);

        /* Convert the scaled magnitude to segment number. */

        seg = search(pcm_val, seg_uend, 8);

        /*

         * Combine the sign, segment, quantization bits;

         * and complement the code word.

         */

        if (seg >= 8)           /* out of range, return maximum value. */

                return (unsigned char) (0x7F ^ mask);

        else {

                uval = (unsigned char) (seg << 4) | ((pcm_val >> (seg + 1)) & 0xF);

                return (uval ^ mask);

        }

}

static PyInt16 _st_alaw2linear16[256] = {

     -5504,   -5248,   -6016,   -5760,   -4480,   -4224,   -4992,

     -4736,   -7552,   -7296,   -8064,   -7808,   -6528,   -6272,

     -7040,   -6784,   -2752,   -2624,   -3008,   -2880,   -2240,

     -2112,   -2496,   -2368,   -3776,   -3648,   -4032,   -3904,

     -3264,   -3136,   -3520,   -3392,  -22016,  -20992,  -24064,

    -23040,  -17920,  -16896,  -19968,  -18944,  -30208,  -29184,

    -32256,  -31232,  -26112,  -25088,  -28160,  -27136,  -11008,

    -10496,  -12032,  -11520,   -8960,   -8448,   -9984,   -9472,

    -15104,  -14592,  -16128,  -15616,  -13056,  -12544,  -14080,

    -13568,    -344,    -328,    -376,    -360,    -280,    -264,

      -312,    -296,    -472,    -456,    -504,    -488,    -408,

      -392,    -440,    -424,     -88,     -72,    -120,    -104,

       -24,      -8,     -56,     -40,    -216,    -200,    -248,

      -232,    -152,    -136,    -184,    -168,   -1376,   -1312,

     -1504,   -1440,   -1120,   -1056,   -1248,   -1184,   -1888,

     -1824,   -2016,   -1952,   -1632,   -1568,   -1760,   -1696,

      -688,    -656,    -752,    -720,    -560,    -528,    -624,

      -592,    -944,    -912,   -1008,    -976,    -816,    -784,

      -880,    -848,    5504,    5248,    6016,    5760,    4480,

      4224,    4992,    4736,    7552,    7296,    8064,    7808,

      6528,    6272,    7040,    6784,    2752,    2624,    3008,

      2880,    2240,    2112,    2496,    2368,    3776,    3648,

      4032,    3904,    3264,    3136,    3520,    3392,   22016,

     20992,   24064,   23040,   17920,   16896,   19968,   18944,

     30208,   29184,   32256,   31232,   26112,   25088,   28160,

     27136,   11008,   10496,   12032,   11520,    8960,    8448,

      9984,    9472,   15104,   14592,   16128,   15616,   13056,

     12544,   14080,   13568,     344,     328,     376,     360,

       280,     264,     312,     296,     472,     456,     504,

       488,     408,     392,     440,     424,      88,      72,

       120,     104,      24,       8,      56,      40,     216,

       200,     248,     232,     152,     136,     184,     168,

      1376,    1312,    1504,    1440,    1120,    1056,    1248,

      1184,    1888,    1824,    2016,    1952,    1632,    1568,

      1760,    1696,     688,     656,     752,     720,     560,

       528,     624,     592,     944,     912,    1008,     976,

       816,     784,     880,     848

};

/*

 * linear2alaw() accepts an 13-bit signed integer and encodes it as A-law data

 * stored in a unsigned char.  This function should only be called with

 * the data shifted such that it only contains information in the lower

 * 13-bits.

 *

 *              Linear Input Code       Compressed Code

 *      ------------------------        ---------------

 *      0000000wxyza                    000wxyz

 *      0000001wxyza                    001wxyz

 *      000001wxyzab                    010wxyz

 *      00001wxyzabc                    011wxyz

 *      0001wxyzabcd                    100wxyz

 *      001wxyzabcde                    101wxyz

 *      01wxyzabcdef                    110wxyz

 *      1wxyzabcdefg                    111wxyz

 *

 * For further information see John C. Bellamy's Digital Telephony, 1982,

 * John Wiley & Sons, pps 98-111 and 472-476.

 */

static unsigned char

st_linear2alaw(PyInt16 pcm_val)	/* 2's complement (13-bit range) */

{

        PyInt16         mask;

        short           seg;

        unsigned char   aval;

        /* The original sox code does this in the calling function, not here */

        pcm_val = pcm_val >> 3;

        /* A-law using even bit inversion */

        if (pcm_val >= 0) {

                mask = 0xD5;            /* sign (7th) bit = 1 */

        } else {

                mask = 0x55;            /* sign bit = 0 */

                pcm_val = -pcm_val - 1;

        }

        /* Convert the scaled magnitude to segment number. */

        seg = search(pcm_val, seg_aend, 8);

        /* Combine the sign, segment, and quantization bits. */

        if (seg >= 8)           /* out of range, return maximum value. */

                return (unsigned char) (0x7F ^ mask);

        else {

                aval = (unsigned char) seg << SEG_SHIFT;

                if (seg < 2)

                        aval |= (pcm_val >> 1) & QUANT_MASK;

                else

                        aval |= (pcm_val >> seg) & QUANT_MASK;

                return (aval ^ mask);

        }

}

/* End of code taken from sox */

/* Intel ADPCM step variation table */

static int indexTable[16] = {

        -1, -1, -1, -1, 2, 4, 6, 8,

        -1, -1, -1, -1, 2, 4, 6, 8,

};

static int stepsizeTable[89] = {

        7, 8, 9, 10, 11, 12, 13, 14, 16, 17,

        19, 21, 23, 25, 28, 31, 34, 37, 41, 45,

        50, 55, 60, 66, 73, 80, 88, 97, 107, 118,

        130, 143, 157, 173, 190, 209, 230, 253, 279, 307,

        337, 371, 408, 449, 494, 544, 598, 658, 724, 796,

        876, 963, 1060, 1166, 1282, 1411, 1552, 1707, 1878, 2066,

        2272, 2499, 2749, 3024, 3327, 3660, 4026, 4428, 4871, 5358,

        5894, 6484, 7132, 7845, 8630, 9493, 10442, 11487, 12635, 13899,

        15289, 16818, 18500, 20350, 22385, 24623, 27086, 29794, 32767

};

#define CHARP(cp, i) ((signed char *)(cp+i))

#define SHORTP(cp, i) ((short *)(cp+i))

#define LONGP(cp, i) ((Py_Int32 *)(cp+i))

static PyObject *AudioopError;

static PyObject *

audioop_getsample(PyObject *self, PyObject *args)

{

        signed char *cp;

        int len, size, val = 0;

        int i;

        if ( !PyArg_ParseTuple(args, "s#ii:getsample", &cp, &len, &size, &i) )

                return 0;

        if ( size != 1 && size != 2 && size != 4 ) {

                PyErr_SetString(AudioopError, "Size should be 1, 2 or 4");

                return 0;

        }

        if ( i < 0 || i >= len/size ) {

                PyErr_SetString(AudioopError, "Index out of range");

                return 0;

        }

        if ( size == 1 )      val = (int)*CHARP(cp, i);

        else if ( size == 2 ) val = (int)*SHORTP(cp, i*2);

        else if ( size == 4 ) val = (int)*LONGP(cp, i*4);

        return PyInt_FromLong(val);

}

static PyObject *

audioop_max(PyObject *self, PyObject *args)

{

        signed char *cp;

        int len, size, val = 0;

        int i;

        int max = 0;

        if ( !PyArg_ParseTuple(args, "s#i:max", &cp, &len, &size) )

                return 0;

        if ( size != 1 && size != 2 && size != 4 ) {

                PyErr_SetString(AudioopError, "Size should be 1, 2 or 4");

                return 0;

        }

        for ( i=0; i<len; i+= size) {

                if ( size == 1 )      val = (int)*CHARP(cp, i);

                else if ( size == 2 ) val = (int)*SHORTP(cp, i);

                else if ( size == 4 ) val = (int)*LONGP(cp, i);

                if ( val < 0 ) val = (-val);

                if ( val > max ) max = val;

        }

        return PyInt_FromLong(max);

}

static PyObject *

audioop_minmax(PyObject *self, PyObject *args)

{

        signed char *cp;

        int len, size, val = 0;

        int i;

        int min = 0x7fffffff, max = -0x7fffffff;

        if (!PyArg_ParseTuple(args, "s#i:minmax", &cp, &len, &size))

                return NULL;

        if (size != 1 && size != 2 && size != 4) {

                PyErr_SetString(AudioopError, "Size should be 1, 2 or 4");

                return NULL;

        }

        for (i = 0; i < len; i += size) {

                if (size == 1) val = (int) *CHARP(cp, i);

                else if (size == 2) val = (int) *SHORTP(cp, i);

                else if (size == 4) val = (int) *LONGP(cp, i);

                if (val > max) max = val;

                if (val < min) min = val;

        }

        return Py_BuildValue("(ii)", min, max);

}

static PyObject *

audioop_avg(PyObject *self, PyObject *args)

{

        signed char *cp;

        int len, size, val = 0;

        int i;

        double avg = 0.0;

        if ( !PyArg_ParseTuple(args, "s#i:avg", &cp, &len, &size) )

                return 0;

        if ( size != 1 && size != 2 && size != 4 ) {

                PyErr_SetString(AudioopError, "Size should be 1, 2 or 4");

                return 0;

        }

        for ( i=0; i<len; i+= size) {

                if ( size == 1 )      val = (int)*CHARP(cp, i);

                else if ( size == 2 ) val = (int)*SHORTP(cp, i);

                else if ( size == 4 ) val = (int)*LONGP(cp, i);

                avg += val;

        }

        if ( len == 0 )

                val = 0;

        else

                val = (int)(avg / (double)(len/size));

        return PyInt_FromLong(val);

}

static PyObject *

audioop_rms(PyObject *self, PyObject *args)

{

        signed char *cp;

        int len, size, val = 0;

        int i;

        double sum_squares = 0.0;

        if ( !PyArg_ParseTuple(args, "s#i:rms", &cp, &len, &size) )

                return 0;

        if ( size != 1 && size != 2 && size != 4 ) {

                PyErr_SetString(AudioopError, "Size should be 1, 2 or 4");

                return 0;

        }

        for ( i=0; i<len; i+= size) {

                if ( size == 1 )      val = (int)*CHARP(cp, i);

                else if ( size == 2 ) val = (int)*SHORTP(cp, i);

                else if ( size == 4 ) val = (int)*LONGP(cp, i);

                sum_squares += (double)val*(double)val;

        }

        if ( len == 0 )

                val = 0;

        else

                val = (int)sqrt(sum_squares / (double)(len/size));

        return PyInt_FromLong(val);

}

static double _sum2(short *a, short *b, int len)

{

        int i;

        double sum = 0.0;

        for( i=0; i<len; i++) {

                sum = sum + (double)a[i]*(double)b[i];

        }

        return sum;

}

/*

** Findfit tries to locate a sample within another sample. Its main use

** is in echo-cancellation (to find the feedback of the output signal in

** the input signal).

** The method used is as follows:

**

** let R be the reference signal (length n) and A the input signal (length N)

** with N > n, and let all sums be over i from 0 to n-1.

**

** Now, for each j in {0..N-n} we compute a factor fj so that -fj*R matches A

** as good as possible, i.e. sum( (A[j+i]+fj*R[i])^2 ) is minimal. This

** equation gives fj = sum( A[j+i]R[i] ) / sum(R[i]^2).

**

** Next, we compute the relative distance between the original signal and

** the modified signal and minimize that over j:

** vj = sum( (A[j+i]-fj*R[i])^2 ) / sum( A[j+i]^2 )  =>

** vj = ( sum(A[j+i]^2)*sum(R[i]^2) - sum(A[j+i]R[i])^2 ) / sum( A[j+i]^2 )

**

** In the code variables correspond as follows:

** cp1          A

** cp2          R

** len1         N

** len2         n

** aj_m1        A[j-1]

** aj_lm1       A[j+n-1]

** sum_ri_2     sum(R[i]^2)

** sum_aij_2    sum(A[i+j]^2)

** sum_aij_ri   sum(A[i+j]R[i])

**

** sum_ri is calculated once, sum_aij_2 is updated each step and sum_aij_ri

** is completely recalculated each step.

*/

static PyObject *

audioop_findfit(PyObject *self, PyObject *args)

{

        short *cp1, *cp2;

        int len1, len2;

        int j, best_j;

        double aj_m1, aj_lm1;

        double sum_ri_2, sum_aij_2, sum_aij_ri, result, best_result, factor;

	/* Passing a short** for an 's' argument is correct only

	   if the string contents is aligned for interpretation

	   as short[]. Due to the definition of PyStringObject,

	   this is currently (Python 2.6) the case. */

        if ( !PyArg_ParseTuple(args, "s#s#:findfit",

	                       (char**)&cp1, &len1, (char**)&cp2, &len2) )

                return 0;

        if ( len1 & 1 || len2 & 1 ) {

                PyErr_SetString(AudioopError, "Strings should be even-sized");

                return 0;

        }

        len1 >>= 1;

        len2 >>= 1;

        if ( len1 < len2 ) {

                PyErr_SetString(AudioopError, "First sample should be longer");

                return 0;

        }

        sum_ri_2 = _sum2(cp2, cp2, len2);

        sum_aij_2 = _sum2(cp1, cp1, len2);

        sum_aij_ri = _sum2(cp1, cp2, len2);

        result = (sum_ri_2*sum_aij_2 - sum_aij_ri*sum_aij_ri) / sum_aij_2;

        best_result = result;

        best_j = 0;

        j = 0;

        for ( j=1; j<=len1-len2; j++) {

                aj_m1 = (double)cp1[j-1];

                aj_lm1 = (double)cp1[j+len2-1];

                sum_aij_2 = sum_aij_2 + aj_lm1*aj_lm1 - aj_m1*aj_m1;

                sum_aij_ri = _sum2(cp1+j, cp2, len2);

                result = (sum_ri_2*sum_aij_2 - sum_aij_ri*sum_aij_ri)

                        / sum_aij_2;

                if ( result < best_result ) {

                        best_result = result;

                        best_j = j;

                }

        }

        factor = _sum2(cp1+best_j, cp2, len2) / sum_ri_2;

        return Py_BuildValue("(if)", best_j, factor);

}

/*

** findfactor finds a factor f so that the energy in A-fB is minimal.

** See the comment for findfit for details.

*/