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
WebRtcSpl_AutoCorrelation(const int16_t * in_vector,size_t in_vector_length,size_t order,int32_t * result,int * scale)15 size_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