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