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
WebRtcSpl_AutoCorrelation(const int16_t * in_vector,int in_vector_length,int order,int32_t * result,int * scale)13 int WebRtcSpl_AutoCorrelation(const int16_t* in_vector,
14 int in_vector_length,
15 int order,
16 int32_t* result,
17 int* scale) {
18 int32_t sum = 0;
19 int i = 0, j = 0;
20 int16_t smax = 0;
21 int scaling = 0;
22
23 if (order > in_vector_length) {
24 /* Undefined */
25 return -1;
26 } else if (order < 0) {
27 order = in_vector_length;
28 }
29
30 // Find the maximum absolute value of the samples.
31 smax = WebRtcSpl_MaxAbsValueW16(in_vector, in_vector_length);
32
33 // In order to avoid overflow when computing the sum we should scale the
34 // samples so that (in_vector_length * smax * smax) will not overflow.
35 if (smax == 0) {
36 scaling = 0;
37 } else {
38 // Number of bits in the sum loop.
39 int nbits = WebRtcSpl_GetSizeInBits(in_vector_length);
40 // Number of bits to normalize smax.
41 int t = WebRtcSpl_NormW32(WEBRTC_SPL_MUL(smax, smax));
42
43 if (t > nbits) {
44 scaling = 0;
45 } else {
46 scaling = nbits - t;
47 }
48 }
49
50 // Perform the actual correlation calculation.
51 for (i = 0; i < order + 1; i++) {
52 sum = 0;
53 /* Unroll the loop to improve performance. */
54 for (j = 0; j < in_vector_length - i - 3; j += 4) {
55 sum += (in_vector[j + 0] * in_vector[i + j + 0]) >> scaling;
56 sum += (in_vector[j + 1] * in_vector[i + j + 1]) >> scaling;
57 sum += (in_vector[j + 2] * in_vector[i + j + 2]) >> scaling;
58 sum += (in_vector[j + 3] * in_vector[i + j + 3]) >> scaling;
59 }
60 for (; j < in_vector_length - i; j++) {
61 sum += (in_vector[j] * in_vector[i + j]) >> scaling;
62 }
63 *result++ = sum;
64 }
65
66 *scale = scaling;
67 return order + 1;
68 }
69