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
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(AudioMultiVector * signal,size_t start_index,size_t length,int factor,int increment)83 int 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
PeakDetection(int16_t * data,size_t data_length,size_t num_peaks,int fs_mult,size_t * peak_index,int16_t * peak_value)102 void DspHelper::PeakDetection(int16_t* data, size_t data_length,
103 size_t num_peaks, int fs_mult,
104 size_t* peak_index, int16_t* peak_value) {
105 size_t min_index = 0;
106 size_t max_index = 0;
107
108 for (size_t 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 = (peak_index[i] > 2) ? (peak_index[i] - 2) : 0;
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
ParabolicFit(int16_t * signal_points,int fs_mult,size_t * peak_index,int16_t * peak_value)150 void DspHelper::ParabolicFit(int16_t* signal_points, int fs_mult,
151 size_t* 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
MinDistortion(const int16_t * signal,size_t min_lag,size_t max_lag,size_t length,int32_t * distortion_value)238 size_t DspHelper::MinDistortion(const int16_t* signal, size_t min_lag,
239 size_t max_lag, size_t length,
240 int32_t* distortion_value) {
241 size_t best_index = 0;
242 int32_t min_distortion = WEBRTC_SPL_WORD32_MAX;
243 for (size_t 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 (size_t 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
CrossFade(const int16_t * input1,const int16_t * input2,size_t length,int16_t * mix_factor,int16_t factor_decrement,int16_t * output)260 void 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
UnmuteSignal(const int16_t * input,size_t length,int16_t * factor,int increment,int16_t * output)274 void DspHelper::UnmuteSignal(const int16_t* input, size_t length,
275 int16_t* factor, int 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
MuteSignal(int16_t * signal,int mute_slope,size_t length)287 void DspHelper::MuteSignal(int16_t* signal, int 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
DownsampleTo4kHz(const int16_t * input,size_t input_length,size_t output_length,int input_rate_hz,bool compensate_delay,int16_t * output)295 int DspHelper::DownsampleTo4kHz(const int16_t* input, size_t input_length,
296 size_t 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 size_t filter_length; // Number of coefficients.
304 size_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], input_length - filter_length + 1, output,
349 output_length, filter_coefficients, filter_length, factor, filter_delay);
350 }
351
352 } // namespace webrtc
353