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