Line data Source code
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 :
60 0 : int DspHelper::RampSignal(const int16_t* input,
61 : size_t length,
62 : int factor,
63 : int increment,
64 : int16_t* output) {
65 0 : int factor_q20 = (factor << 6) + 32;
66 : // TODO(hlundin): Add 32 to factor_q20 when converting back to Q14?
67 0 : for (size_t i = 0; i < length; ++i) {
68 0 : output[i] = (factor * input[i] + 8192) >> 14;
69 0 : factor_q20 += increment;
70 0 : factor_q20 = std::max(factor_q20, 0); // Never go negative.
71 0 : factor = std::min(factor_q20 >> 6, 16384);
72 : }
73 0 : return factor;
74 : }
75 :
76 0 : int DspHelper::RampSignal(int16_t* signal,
77 : size_t length,
78 : int factor,
79 : int increment) {
80 0 : return RampSignal(signal, length, factor, increment, signal);
81 : }
82 :
83 0 : int DspHelper::RampSignal(AudioVector* signal,
84 : size_t start_index,
85 : size_t length,
86 : int factor,
87 : int increment) {
88 0 : int factor_q20 = (factor << 6) + 32;
89 : // TODO(hlundin): Add 32 to factor_q20 when converting back to Q14?
90 0 : for (size_t i = start_index; i < start_index + length; ++i) {
91 0 : (*signal)[i] = (factor * (*signal)[i] + 8192) >> 14;
92 0 : factor_q20 += increment;
93 0 : factor_q20 = std::max(factor_q20, 0); // Never go negative.
94 0 : factor = std::min(factor_q20 >> 6, 16384);
95 : }
96 0 : return factor;
97 : }
98 :
99 0 : int DspHelper::RampSignal(AudioMultiVector* signal,
100 : size_t start_index,
101 : size_t length,
102 : int factor,
103 : int increment) {
104 0 : assert(start_index + length <= signal->Size());
105 0 : if (start_index + length > signal->Size()) {
106 : // Wrong parameters. Do nothing and return the scale factor unaltered.
107 0 : return factor;
108 : }
109 0 : int end_factor = 0;
110 : // Loop over the channels, starting at the same |factor| each time.
111 0 : for (size_t channel = 0; channel < signal->Channels(); ++channel) {
112 : end_factor =
113 0 : RampSignal(&(*signal)[channel], start_index, length, factor, increment);
114 : }
115 0 : return end_factor;
116 : }
117 :
118 0 : 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 0 : size_t min_index = 0;
122 0 : size_t max_index = 0;
123 :
124 0 : for (size_t i = 0; i <= num_peaks - 1; i++) {
125 0 : 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 0 : data_length++;
131 : }
132 :
133 0 : peak_index[i] = WebRtcSpl_MaxIndexW16(data, data_length - 1);
134 :
135 0 : if (i != num_peaks - 1) {
136 0 : min_index = (peak_index[i] > 2) ? (peak_index[i] - 2) : 0;
137 0 : max_index = std::min(data_length - 1, peak_index[i] + 2);
138 : }
139 :
140 0 : if ((peak_index[i] != 0) && (peak_index[i] != (data_length - 2))) {
141 0 : ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
142 0 : &peak_value[i]);
143 : } else {
144 0 : if (peak_index[i] == data_length - 2) {
145 0 : if (data[peak_index[i]] > data[peak_index[i] + 1]) {
146 0 : ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
147 0 : &peak_value[i]);
148 0 : } else if (data[peak_index[i]] <= data[peak_index[i] + 1]) {
149 : // Linear approximation.
150 0 : peak_value[i] = (data[peak_index[i]] + data[peak_index[i] + 1]) >> 1;
151 0 : peak_index[i] = (peak_index[i] * 2 + 1) * fs_mult;
152 : }
153 : } else {
154 0 : peak_value[i] = data[peak_index[i]];
155 0 : peak_index[i] = peak_index[i] * 2 * fs_mult;
156 : }
157 : }
158 :
159 0 : if (i != num_peaks - 1) {
160 0 : memset(&data[min_index], 0,
161 0 : sizeof(data[0]) * (max_index - min_index + 1));
162 : }
163 : }
164 0 : }
165 :
166 0 : 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 0 : if (fs_mult == 1) {
170 0 : fit_index[0] = 0;
171 0 : fit_index[1] = 8;
172 0 : fit_index[2] = 16;
173 0 : } else if (fs_mult == 2) {
174 0 : fit_index[0] = 0;
175 0 : fit_index[1] = 4;
176 0 : fit_index[2] = 8;
177 0 : fit_index[3] = 12;
178 0 : fit_index[4] = 16;
179 0 : } else if (fs_mult == 4) {
180 0 : fit_index[0] = 0;
181 0 : fit_index[1] = 2;
182 0 : fit_index[2] = 4;
183 0 : fit_index[3] = 6;
184 0 : fit_index[4] = 8;
185 0 : fit_index[5] = 10;
186 0 : fit_index[6] = 12;
187 0 : fit_index[7] = 14;
188 0 : fit_index[8] = 16;
189 : } else {
190 0 : fit_index[0] = 0;
191 0 : fit_index[1] = 1;
192 0 : fit_index[2] = 3;
193 0 : fit_index[3] = 4;
194 0 : fit_index[4] = 5;
195 0 : fit_index[5] = 7;
196 0 : fit_index[6] = 8;
197 0 : fit_index[7] = 9;
198 0 : fit_index[8] = 11;
199 0 : fit_index[9] = 12;
200 0 : fit_index[10] = 13;
201 0 : fit_index[11] = 15;
202 0 : 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 0 : int32_t num = (signal_points[0] * -3) + (signal_points[1] * 4)
208 0 : - signal_points[2];
209 0 : int32_t den = signal_points[0] + (signal_points[1] * -2) + signal_points[2];
210 0 : int32_t temp = num * 120;
211 0 : int flag = 1;
212 0 : int16_t stp = kParabolaCoefficients[fit_index[fs_mult]][0]
213 0 : - kParabolaCoefficients[fit_index[fs_mult - 1]][0];
214 0 : int16_t strt = (kParabolaCoefficients[fit_index[fs_mult]][0]
215 0 : + kParabolaCoefficients[fit_index[fs_mult - 1]][0]) / 2;
216 : int16_t lmt;
217 0 : if (temp < -den * strt) {
218 0 : lmt = strt - stp;
219 0 : while (flag) {
220 0 : if ((flag == fs_mult) || (temp > -den * lmt)) {
221 0 : *peak_value = (den * kParabolaCoefficients[fit_index[fs_mult - flag]][1]
222 0 : + num * kParabolaCoefficients[fit_index[fs_mult - flag]][2]
223 0 : + signal_points[0] * 256) / 256;
224 0 : *peak_index = *peak_index * 2 * fs_mult - flag;
225 0 : flag = 0;
226 : } else {
227 0 : flag++;
228 0 : lmt -= stp;
229 : }
230 : }
231 0 : } else if (temp > -den * (strt + stp)) {
232 0 : lmt = strt + 2 * stp;
233 0 : while (flag) {
234 0 : if ((flag == fs_mult) || (temp < -den * lmt)) {
235 : int32_t temp_term_1 =
236 0 : den * kParabolaCoefficients[fit_index[fs_mult+flag]][1];
237 : int32_t temp_term_2 =
238 0 : num * kParabolaCoefficients[fit_index[fs_mult+flag]][2];
239 0 : int32_t temp_term_3 = signal_points[0] * 256;
240 0 : *peak_value = (temp_term_1 + temp_term_2 + temp_term_3) / 256;
241 0 : *peak_index = *peak_index * 2 * fs_mult + flag;
242 0 : flag = 0;
243 : } else {
244 0 : flag++;
245 0 : lmt += stp;
246 : }
247 : }
248 : } else {
249 0 : *peak_value = signal_points[1];
250 0 : *peak_index = *peak_index * 2 * fs_mult;
251 : }
252 0 : }
253 :
254 0 : 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 0 : size_t best_index = 0;
258 0 : int32_t min_distortion = WEBRTC_SPL_WORD32_MAX;
259 0 : for (size_t i = min_lag; i <= max_lag; i++) {
260 0 : int32_t sum_diff = 0;
261 0 : const int16_t* data1 = signal;
262 0 : const int16_t* data2 = signal - i;
263 0 : for (size_t j = 0; j < length; j++) {
264 0 : sum_diff += WEBRTC_SPL_ABS_W32(data1[j] - data2[j]);
265 : }
266 : // Compare with previous minimum.
267 0 : if (sum_diff < min_distortion) {
268 0 : min_distortion = sum_diff;
269 0 : best_index = i;
270 : }
271 : }
272 0 : *distortion_value = min_distortion;
273 0 : return best_index;
274 : }
275 :
276 0 : 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 0 : int16_t factor = *mix_factor;
280 0 : int16_t complement_factor = 16384 - factor;
281 0 : for (size_t i = 0; i < length; i++) {
282 0 : output[i] =
283 0 : (factor * input1[i] + complement_factor * input2[i] + 8192) >> 14;
284 0 : factor -= factor_decrement;
285 0 : complement_factor += factor_decrement;
286 : }
287 0 : *mix_factor = factor;
288 0 : }
289 :
290 0 : void DspHelper::UnmuteSignal(const int16_t* input, size_t length,
291 : int16_t* factor, int increment,
292 : int16_t* output) {
293 0 : uint16_t factor_16b = *factor;
294 0 : int32_t factor_32b = (static_cast<int32_t>(factor_16b) << 6) + 32;
295 0 : for (size_t i = 0; i < length; i++) {
296 0 : output[i] = (factor_16b * input[i] + 8192) >> 14;
297 0 : factor_32b = std::max(factor_32b + increment, 0);
298 0 : factor_16b = std::min(16384, factor_32b >> 6);
299 : }
300 0 : *factor = factor_16b;
301 0 : }
302 :
303 0 : void DspHelper::MuteSignal(int16_t* signal, int mute_slope, size_t length) {
304 0 : int32_t factor = (16384 << 6) + 32;
305 0 : for (size_t i = 0; i < length; i++) {
306 0 : signal[i] = ((factor >> 6) * signal[i] + 8192) >> 14;
307 0 : factor -= mute_slope;
308 : }
309 0 : }
310 :
311 0 : 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 0 : switch (input_rate_hz) {
323 : case 8000: {
324 0 : filter_length = 3;
325 0 : factor = 2;
326 0 : filter_coefficients = kDownsample8kHzTbl;
327 0 : filter_delay = 1 + 1;
328 0 : break;
329 : }
330 : case 16000: {
331 0 : filter_length = 5;
332 0 : factor = 4;
333 0 : filter_coefficients = kDownsample16kHzTbl;
334 0 : filter_delay = 2 + 1;
335 0 : break;
336 : }
337 : case 32000: {
338 0 : filter_length = 7;
339 0 : factor = 8;
340 0 : filter_coefficients = kDownsample32kHzTbl;
341 0 : filter_delay = 3 + 1;
342 0 : break;
343 : }
344 : case 48000: {
345 0 : filter_length = 7;
346 0 : factor = 12;
347 0 : filter_coefficients = kDownsample48kHzTbl;
348 0 : filter_delay = 3 + 1;
349 0 : break;
350 : }
351 : default: {
352 0 : assert(false);
353 : return -1;
354 : }
355 : }
356 :
357 0 : if (!compensate_delay) {
358 : // Disregard delay compensation.
359 0 : filter_delay = 0;
360 : }
361 :
362 : // Returns -1 if input signal is too short; 0 otherwise.
363 0 : return WebRtcSpl_DownsampleFast(
364 0 : &input[filter_length - 1], input_length - filter_length + 1, output,
365 0 : output_length, filter_coefficients, filter_length, factor, filter_delay);
366 : }
367 :
368 : } // namespace webrtc
|
