1/*
2 *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
3 *
4 *  Use of this source code is governed by a BSD-style license
5 *  that can be found in the LICENSE file in the root of the source
6 *  tree. An additional intellectual property rights grant can be found
7 *  in the file PATENTS.  All contributing project authors may
8 *  be found in the AUTHORS file in the root of the source tree.
9 */
10
11#include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
12
13#include <assert.h>
14
15size_t WebRtcSpl_AutoCorrelation(const int16_t* in_vector,
16                                 size_t in_vector_length,
17                                 size_t order,
18                                 int32_t* result,
19                                 int* scale) {
20  int32_t sum = 0;
21  size_t i = 0, j = 0;
22  int16_t smax = 0;
23  int scaling = 0;
24
25  assert(order <= in_vector_length);
26
27  // Find the maximum absolute value of the samples.
28  smax = WebRtcSpl_MaxAbsValueW16(in_vector, in_vector_length);
29
30  // In order to avoid overflow when computing the sum we should scale the
31  // samples so that (in_vector_length * smax * smax) will not overflow.
32  if (smax == 0) {
33    scaling = 0;
34  } else {
35    // Number of bits in the sum loop.
36    int nbits = WebRtcSpl_GetSizeInBits((uint32_t)in_vector_length);
37    // Number of bits to normalize smax.
38    int t = WebRtcSpl_NormW32(WEBRTC_SPL_MUL(smax, smax));
39
40    if (t > nbits) {
41      scaling = 0;
42    } else {
43      scaling = nbits - t;
44    }
45  }
46
47  // Perform the actual correlation calculation.
48  for (i = 0; i < order + 1; i++) {
49    sum = 0;
50    /* Unroll the loop to improve performance. */
51    for (j = 0; i + j + 3 < in_vector_length; j += 4) {
52      sum += (in_vector[j + 0] * in_vector[i + j + 0]) >> scaling;
53      sum += (in_vector[j + 1] * in_vector[i + j + 1]) >> scaling;
54      sum += (in_vector[j + 2] * in_vector[i + j + 2]) >> scaling;
55      sum += (in_vector[j + 3] * in_vector[i + j + 3]) >> scaling;
56    }
57    for (; j < in_vector_length - i; j++) {
58      sum += (in_vector[j] * in_vector[i + j]) >> scaling;
59    }
60    *result++ = sum;
61  }
62
63  *scale = scaling;
64  return order + 1;
65}
66