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