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/modules/audio_coding/neteq/dsp_helper.h"
12
13#include <assert.h>
14#include <string.h>  // Access to memset.
15
16#include <algorithm>  // Access to min, max.
17
18#include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
19
20namespace webrtc {
21
22// Table of constants used in method DspHelper::ParabolicFit().
23const int16_t DspHelper::kParabolaCoefficients[17][3] = {
24    { 120, 32, 64 },
25    { 140, 44, 75 },
26    { 150, 50, 80 },
27    { 160, 57, 85 },
28    { 180, 72, 96 },
29    { 200, 89, 107 },
30    { 210, 98, 112 },
31    { 220, 108, 117 },
32    { 240, 128, 128 },
33    { 260, 150, 139 },
34    { 270, 162, 144 },
35    { 280, 174, 149 },
36    { 300, 200, 160 },
37    { 320, 228, 171 },
38    { 330, 242, 176 },
39    { 340, 257, 181 },
40    { 360, 288, 192 } };
41
42// Filter coefficients used when downsampling from the indicated sample rates
43// (8, 16, 32, 48 kHz) to 4 kHz. Coefficients are in Q12. The corresponding Q0
44// values are provided in the comments before each array.
45
46// Q0 values: {0.3, 0.4, 0.3}.
47const int16_t DspHelper::kDownsample8kHzTbl[3] = { 1229, 1638, 1229 };
48
49// Q0 values: {0.15, 0.2, 0.3, 0.2, 0.15}.
50const int16_t DspHelper::kDownsample16kHzTbl[5] = { 614, 819, 1229, 819, 614 };
51
52// Q0 values: {0.1425, 0.1251, 0.1525, 0.1628, 0.1525, 0.1251, 0.1425}.
53const int16_t DspHelper::kDownsample32kHzTbl[7] = {
54    584, 512, 625, 667, 625, 512, 584 };
55
56// Q0 values: {0.2487, 0.0952, 0.1042, 0.1074, 0.1042, 0.0952, 0.2487}.
57const int16_t DspHelper::kDownsample48kHzTbl[7] = {
58    1019, 390, 427, 440, 427, 390, 1019 };
59
60int DspHelper::RampSignal(const int16_t* input,
61                          size_t length,
62                          int factor,
63                          int increment,
64                          int16_t* output) {
65  int factor_q20 = (factor << 6) + 32;
66  // TODO(hlundin): Add 32 to factor_q20 when converting back to Q14?
67  for (size_t i = 0; i < length; ++i) {
68    output[i] = (factor * input[i] + 8192) >> 14;
69    factor_q20 += increment;
70    factor_q20 = std::max(factor_q20, 0);  // Never go negative.
71    factor = std::min(factor_q20 >> 6, 16384);
72  }
73  return factor;
74}
75
76int DspHelper::RampSignal(int16_t* signal,
77                          size_t length,
78                          int factor,
79                          int increment) {
80  return RampSignal(signal, length, factor, increment, signal);
81}
82
83int DspHelper::RampSignal(AudioMultiVector* signal,
84                          size_t start_index,
85                          size_t length,
86                          int factor,
87                          int increment) {
88  assert(start_index + length <= signal->Size());
89  if (start_index + length > signal->Size()) {
90    // Wrong parameters. Do nothing and return the scale factor unaltered.
91    return factor;
92  }
93  int end_factor = 0;
94  // Loop over the channels, starting at the same |factor| each time.
95  for (size_t channel = 0; channel < signal->Channels(); ++channel) {
96    end_factor =
97        RampSignal(&(*signal)[channel][start_index], length, factor, increment);
98  }
99  return end_factor;
100}
101
102void DspHelper::PeakDetection(int16_t* data, int data_length,
103                              int num_peaks, int fs_mult,
104                              int* peak_index, int16_t* peak_value) {
105  int16_t min_index = 0;
106  int16_t max_index = 0;
107
108  for (int i = 0; i <= num_peaks - 1; i++) {
109    if (num_peaks == 1) {
110      // Single peak.  The parabola fit assumes that an extra point is
111      // available; worst case it gets a zero on the high end of the signal.
112      // TODO(hlundin): This can potentially get much worse. It breaks the
113      // API contract, that the length of |data| is |data_length|.
114      data_length++;
115    }
116
117    peak_index[i] = WebRtcSpl_MaxIndexW16(data, data_length - 1);
118
119    if (i != num_peaks - 1) {
120      min_index = std::max(0, peak_index[i] - 2);
121      max_index = std::min(data_length - 1, peak_index[i] + 2);
122    }
123
124    if ((peak_index[i] != 0) && (peak_index[i] != (data_length - 2))) {
125      ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
126                   &peak_value[i]);
127    } else {
128      if (peak_index[i] == data_length - 2) {
129        if (data[peak_index[i]] > data[peak_index[i] + 1]) {
130          ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
131                       &peak_value[i]);
132        } else if (data[peak_index[i]] <= data[peak_index[i] + 1]) {
133          // Linear approximation.
134          peak_value[i] = (data[peak_index[i]] + data[peak_index[i] + 1]) >> 1;
135          peak_index[i] = (peak_index[i] * 2 + 1) * fs_mult;
136        }
137      } else {
138        peak_value[i] = data[peak_index[i]];
139        peak_index[i] = peak_index[i] * 2 * fs_mult;
140      }
141    }
142
143    if (i != num_peaks - 1) {
144      memset(&data[min_index], 0,
145             sizeof(data[0]) * (max_index - min_index + 1));
146    }
147  }
148}
149
150void DspHelper::ParabolicFit(int16_t* signal_points, int fs_mult,
151                             int* peak_index, int16_t* peak_value) {
152  uint16_t fit_index[13];
153  if (fs_mult == 1) {
154    fit_index[0] = 0;
155    fit_index[1] = 8;
156    fit_index[2] = 16;
157  } else if (fs_mult == 2) {
158    fit_index[0] = 0;
159    fit_index[1] = 4;
160    fit_index[2] = 8;
161    fit_index[3] = 12;
162    fit_index[4] = 16;
163  } else if (fs_mult == 4) {
164    fit_index[0] = 0;
165    fit_index[1] = 2;
166    fit_index[2] = 4;
167    fit_index[3] = 6;
168    fit_index[4] = 8;
169    fit_index[5] = 10;
170    fit_index[6] = 12;
171    fit_index[7] = 14;
172    fit_index[8] = 16;
173  } else {
174    fit_index[0] = 0;
175    fit_index[1] = 1;
176    fit_index[2] = 3;
177    fit_index[3] = 4;
178    fit_index[4] = 5;
179    fit_index[5] = 7;
180    fit_index[6] = 8;
181    fit_index[7] = 9;
182    fit_index[8] = 11;
183    fit_index[9] = 12;
184    fit_index[10] = 13;
185    fit_index[11] = 15;
186    fit_index[12] = 16;
187  }
188
189  //  num = -3 * signal_points[0] + 4 * signal_points[1] - signal_points[2];
190  //  den =      signal_points[0] - 2 * signal_points[1] + signal_points[2];
191  int32_t num = (signal_points[0] * -3) + (signal_points[1] * 4)
192      - signal_points[2];
193  int32_t den = signal_points[0] + (signal_points[1] * -2) + signal_points[2];
194  int32_t temp = num * 120;
195  int flag = 1;
196  int16_t stp = kParabolaCoefficients[fit_index[fs_mult]][0]
197      - kParabolaCoefficients[fit_index[fs_mult - 1]][0];
198  int16_t strt = (kParabolaCoefficients[fit_index[fs_mult]][0]
199      + kParabolaCoefficients[fit_index[fs_mult - 1]][0]) / 2;
200  int16_t lmt;
201  if (temp < -den * strt) {
202    lmt = strt - stp;
203    while (flag) {
204      if ((flag == fs_mult) || (temp > -den * lmt)) {
205        *peak_value = (den * kParabolaCoefficients[fit_index[fs_mult - flag]][1]
206            + num * kParabolaCoefficients[fit_index[fs_mult - flag]][2]
207            + signal_points[0] * 256) / 256;
208        *peak_index = *peak_index * 2 * fs_mult - flag;
209        flag = 0;
210      } else {
211        flag++;
212        lmt -= stp;
213      }
214    }
215  } else if (temp > -den * (strt + stp)) {
216    lmt = strt + 2 * stp;
217    while (flag) {
218      if ((flag == fs_mult) || (temp < -den * lmt)) {
219        int32_t temp_term_1 =
220            den * kParabolaCoefficients[fit_index[fs_mult+flag]][1];
221        int32_t temp_term_2 =
222            num * kParabolaCoefficients[fit_index[fs_mult+flag]][2];
223        int32_t temp_term_3 = signal_points[0] * 256;
224        *peak_value = (temp_term_1 + temp_term_2 + temp_term_3) / 256;
225        *peak_index = *peak_index * 2 * fs_mult + flag;
226        flag = 0;
227      } else {
228        flag++;
229        lmt += stp;
230      }
231    }
232  } else {
233    *peak_value = signal_points[1];
234    *peak_index = *peak_index * 2 * fs_mult;
235  }
236}
237
238int DspHelper::MinDistortion(const int16_t* signal, int min_lag,
239                             int max_lag, int length,
240                             int32_t* distortion_value) {
241  int best_index = -1;
242  int32_t min_distortion = WEBRTC_SPL_WORD32_MAX;
243  for (int i = min_lag; i <= max_lag; i++) {
244    int32_t sum_diff = 0;
245    const int16_t* data1 = signal;
246    const int16_t* data2 = signal - i;
247    for (int j = 0; j < length; j++) {
248      sum_diff += WEBRTC_SPL_ABS_W32(data1[j] - data2[j]);
249    }
250    // Compare with previous minimum.
251    if (sum_diff < min_distortion) {
252      min_distortion = sum_diff;
253      best_index = i;
254    }
255  }
256  *distortion_value = min_distortion;
257  return best_index;
258}
259
260void DspHelper::CrossFade(const int16_t* input1, const int16_t* input2,
261                          size_t length, int16_t* mix_factor,
262                          int16_t factor_decrement, int16_t* output) {
263  int16_t factor = *mix_factor;
264  int16_t complement_factor = 16384 - factor;
265  for (size_t i = 0; i < length; i++) {
266    output[i] =
267        (factor * input1[i] + complement_factor * input2[i] + 8192) >> 14;
268    factor -= factor_decrement;
269    complement_factor += factor_decrement;
270  }
271  *mix_factor = factor;
272}
273
274void DspHelper::UnmuteSignal(const int16_t* input, size_t length,
275                             int16_t* factor, int16_t increment,
276                             int16_t* output) {
277  uint16_t factor_16b = *factor;
278  int32_t factor_32b = (static_cast<int32_t>(factor_16b) << 6) + 32;
279  for (size_t i = 0; i < length; i++) {
280    output[i] = (factor_16b * input[i] + 8192) >> 14;
281    factor_32b = std::max(factor_32b + increment, 0);
282    factor_16b = std::min(16384, factor_32b >> 6);
283  }
284  *factor = factor_16b;
285}
286
287void DspHelper::MuteSignal(int16_t* signal, int16_t mute_slope, size_t length) {
288  int32_t factor = (16384 << 6) + 32;
289  for (size_t i = 0; i < length; i++) {
290    signal[i] = ((factor >> 6) * signal[i] + 8192) >> 14;
291    factor -= mute_slope;
292  }
293}
294
295int DspHelper::DownsampleTo4kHz(const int16_t* input, size_t input_length,
296                                int output_length, int input_rate_hz,
297                                bool compensate_delay, int16_t* output) {
298  // Set filter parameters depending on input frequency.
299  // NOTE: The phase delay values are wrong compared to the true phase delay
300  // of the filters. However, the error is preserved (through the +1 term) for
301  // consistency.
302  const int16_t* filter_coefficients;  // Filter coefficients.
303  int16_t filter_length;  // Number of coefficients.
304  int16_t filter_delay;  // Phase delay in samples.
305  int16_t factor;  // Conversion rate (inFsHz / 8000).
306  switch (input_rate_hz) {
307    case 8000: {
308      filter_length = 3;
309      factor = 2;
310      filter_coefficients = kDownsample8kHzTbl;
311      filter_delay = 1 + 1;
312      break;
313    }
314    case 16000: {
315      filter_length = 5;
316      factor = 4;
317      filter_coefficients = kDownsample16kHzTbl;
318      filter_delay = 2 + 1;
319      break;
320    }
321    case 32000: {
322      filter_length = 7;
323      factor = 8;
324      filter_coefficients = kDownsample32kHzTbl;
325      filter_delay = 3 + 1;
326      break;
327    }
328    case 48000: {
329      filter_length = 7;
330      factor = 12;
331      filter_coefficients = kDownsample48kHzTbl;
332      filter_delay = 3 + 1;
333      break;
334    }
335    default: {
336      assert(false);
337      return -1;
338    }
339  }
340
341  if (!compensate_delay) {
342    // Disregard delay compensation.
343    filter_delay = 0;
344  }
345
346  // Returns -1 if input signal is too short; 0 otherwise.
347  return WebRtcSpl_DownsampleFast(
348      &input[filter_length - 1], static_cast<int>(input_length) -
349      (filter_length - 1), output, output_length, filter_coefficients,
350      filter_length, factor, filter_delay);
351}
352
353}  // namespace webrtc
354