3
$\begingroup$

Im trying to calculate the autocorrelation of soundwaves when I noticed that I get different results with scipys FFT based and with numpys methods.

The 4 functions Im using:

def c1(x):
    return np.correlate(x,x,'full')
def c2(x):
    return np.convolve(x,x[::-1])
def c3(x):
    return np.int16(np.rint(signal.fftconvolve(x,x[::-1])))
def c4(x):
    return np.int16(np.rint(signal.correlate(x,x,method='fft')))

Its all good for simple arrays, I even created a method that checks for large random numpy arrays, they give back the same result: 

for i in range(1000):
    x = np.random.randint(-200,200,100000)
    r1,r2,r3,r4 = c1(x),c2(x),c3(x),c4(x)
    if not (np.allclose(r1,r2) and np.allclose(r1,r3) and np.allclose(r1,r4)):
        return x

It hasn't failed once. However when I try them on real sound data (read from a wav by scipy) it turns out that c1(x)=c2(x)!=c3(x)=c4(x). I also noticed that c3(x) and c4(x) has the same exact values as c1(x) and c2(x), except that they have large chunks of 0-s in many places (especially in the center) where c1 and c2 dont. Apart from those missing values they are the same.

Can someone tell me what am I doing wrong? Why are there missing values (0-s) in c3 and c4, and why does this happen only when Im processing real wavs and not random arrays? Thanks.

EDIT

Here is an example array (600 in length) for which the functions' results dont match:

array([ 1048,  1052,  1122,  1066,   972,   992,  1086,  1072,  1099,
        1308,  1373,  1388,  1581,  1749,  1781,  1928,  2122,  2158,
        2308,  2539,  2654,  2787,  3071,  3193,  3196,  3328,  3154,
        2896,  2771,  2529,  2231,  1952,  1616,  1157,   760,   406,
         -31,  -518,  -897, -1229, -1512, -1769, -2099, -2420, -2755,
       -3124, -3486, -3697, -3821, -3965, -3970, -3914, -3932, -3935,
       -3952, -3971, -3974, -3913, -3826, -3741, -3621, -3529, -3465,
       -3381, -3257, -3097, -2875, -2562, -2299, -2069, -1805, -1598,
       -1388, -1197,  -931,  -645,  -355,   -10,   307,   610,   840,
         989,  1153,  1356,  1471,  1567,  1719,  1809,  1760,  1700,
        1649,  1552,  1476,  1407,  1292,  1187,  1049,   886,   741,
         529,   299,    94,  -107,  -302,  -454,  -613,  -752,  -909,
       -1129, -1288, -1442, -1604, -1689, -1806, -1959, -2042, -2115,
       -2223, -2271, -2334, -2380, -2383, -2335, -2243, -2170, -2114,
       -2063, -1985, -1927, -1795, -1689, -1567, -1402, -1235, -1034,
        -858,  -720,  -639,  -528,  -396,  -249,   -92,    36,   153,
         253,   323,   363,   456,   520,   572,   667,   741,   825,
         818,   821,   880,   876,   884,   928,   958,   985,   990,
         963,   902,   858,   776,   723,   683,   627,   582,   509,
         466,   383,   305,   219,   142,   116,    85,    54,    17,
         -17,   -72,   -81,   -46,   -39,     3,    84,   113,   130,
         179,   222,   261,   296,   310,   356,   432,   474,   558,
         623,   633,   689,   787,   846,   882,   957,  1027,  1098,
        1169,  1202,  1233,  1271,  1291,  1376,  1378,  1370,  1419,
        1377,  1359,  1365,  1332,  1301,  1285,  1177,  1128,  1108,
        1009,   944,   909,   844,   759,   722,   655,   616,   564,
         520,   510,   514,   467,   410,   428,   396,   410,   475,
         531,   567,   609,   645,   715,   787,   865,   943,   993,
        1111,  1215,  1333,  1416,  1412,  1450,  1669,  1880,  1933,
        2038,  2200,  2325,  2417,  2528,  2683,  2854,  3052,  3276,
        3504,  3697,  3917,  4163,  4338,  4313,  4098,  3846,  3571,
        3100,  2475,  2011,  1621,  1125,   763,   598,   339,   -62,
        -528, -1038, -1606, -2171, -2748, -3271, -3671, -4002, -4294,
       -4436, -4478, -4600, -4672, -4583, -4451, -4454, -4449, -4383,
       -4371, -4416, -4513, -4553, -4499, -4378, -4187, -3933, -3642,
       -3432, -3192, -2966, -2804, -2696, -2557, -2282, -1996, -1720,
       -1396, -1013,  -680,  -362,    64,   518,   918,  1301,  1580,
        1798,  1987,  2071,  2131,  2220,  2305,  2348,  2394,  2472,
        2473,  2403,  2327,  2222,  2009,  1745,  1536,  1310,  1068,
         813,   540,   289,    87,   -87,  -261,  -439,  -597,  -747,
        -880, -1055, -1246, -1407, -1595, -1748, -1868, -1967, -2043,
       -2078, -2112, -2206, -2272, -2316, -2365, -2440, -2486, -2516,
       -2531, -2534, -2491, -2393, -2309, -2163, -1981, -1814, -1638,
       -1433, -1236, -1102,  -976,  -824,  -656,  -472,  -253,   -34,
         111,   297,   456,   543,   620,   748,   834,   859,   890,
         915,   925,   942,   962,   972,   982,  1003,  1029,   999,
         968,   904,   848,   757,   677,   602,   536,   470,   418,
         375,   282,   209,   138,    55,   -47,  -157,  -243,  -330,
        -429,  -468,  -477,  -493,  -517,  -509,  -510,  -467,  -357,
        -289,  -251,  -201,  -125,   -59,    -1,    97,   206,   334,
         463,   597,   700,   783,   897,   981,  1030,  1030,  1032,
        1094,  1153,  1176,  1193,  1187,  1173,  1191,  1223,  1228,
        1220,  1222,  1189,  1102,  1046,  1042,  1049,  1007,   982,
         952,   885,   877,   907,   839,   746,   716,   652,   638,
         628,   540,   483,   533,   474,   434,   461,   448,   469,
         420,   410,   514,   527,   444,   450,   456,   461,   546,
         653,   778,   962,  1115,  1284,  1369,  1365,  1465,  1588,
        1632,  1657,  1737,  1804,  1973,  2118,  2157,  2257,  2365,
        2465,  2497,  2626,  2750,  2765,  2696,  2755,  2960,  3020,
        3094,  3315,  3632,  3837,  3980,  4042,  3927,  3592,  3067,
        2474,  1995,  1582,  1225,  1047,   905,   778,   492,    30,
        -486, -1079, -1802, -2394, -2845, -3289, -3617, -3865, -4141,
       -4392, -4595, -4756, -4711, -4586, -4440, -4328, -4287, -4335,
       -4466, -4543, -4581, -4571, -4445, -4179, -3824, -3514, -3265,
       -3093, -2956, -2859, -2796, -2668, -2473, -2224, -2034, -1852,
       -1645, -1415, -1166,  -778,  -338,    36,   443,   788,  1051,
        1242,  1355,  1391,  1486,  1661,  1869], dtype=int16)

Running a simple checking on it we get:

i:  594 c1(x)[i]:  -11191 c4(x)[i]:  0
i:  595 c1(x)[i]:  8216 c4(x)[i]:  0
i:  596 c1(x)[i]:  -32326 c4(x)[i]:  0
i:  597 c1(x)[i]:  -30734 c4(x)[i]:  0
i:  598 c1(x)[i]:  -23815 c4(x)[i]:  0
i:  599 c1(x)[i]:  4296 c4(x)[i]:  0
i:  600 c1(x)[i]:  -23815 c4(x)[i]:  0
i:  601 c1(x)[i]:  -30734 c4(x)[i]:  0
i:  602 c1(x)[i]:  -32326 c4(x)[i]:  0
i:  603 c1(x)[i]:  8216 c4(x)[i]:  0
i:  604 c1(x)[i]:  -11191 c4(x)[i]:  0

As you can see in the middle of the autocorrelation in symmetric pattern there are zeros where there shouldnt be.

Improve this question
$\endgroup$
2
  • $\begingroup$ you're possibly not doing something wrong. tell us more about the x for which things aren't equal, and what the difference is to those where they are. $\endgroup$
    – Marcus Müller
    Commented Jul 17, 2017 at 23:07
  • $\begingroup$ @MarcusMüller Thanks for asking, I've added an example array. $\endgroup$
    – v.s.
    Commented Jul 18, 2017 at 11:59

1 Answer 1

Reset to default
3
$\begingroup$

By converting to np.int16, you're hitting integer overflow and wrapping the values around. I assume this is not intentional, since the output looks much more sane without it.

int16 is limited to values from −32,768 through +32,767, but the output of the convolution goes from −1,255,891,038 to +2,459,046,088:

np.int16:

wrapped

np.int64:

not wrapped

You can make your test array float or int64 format to handle the required level, and remove the np.int16 conversion after convolution, then all 4 methods are identical.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.