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