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