expand.cc (39446B)
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/expand.h" 12 13 #include <algorithm> // min, max 14 #include <cstdint> 15 #include <cstring> // memset 16 #include <limits> // numeric_limits<T> 17 #include <memory> 18 19 #include "common_audio/signal_processing/dot_product_with_scale.h" 20 #include "common_audio/signal_processing/include/signal_processing_library.h" 21 #include "common_audio/signal_processing/include/spl_inl.h" 22 #include "modules/audio_coding/neteq/audio_multi_vector.h" 23 #include "modules/audio_coding/neteq/background_noise.h" 24 #include "modules/audio_coding/neteq/cross_correlation.h" 25 #include "modules/audio_coding/neteq/dsp_helper.h" 26 #include "modules/audio_coding/neteq/random_vector.h" 27 #include "modules/audio_coding/neteq/statistics_calculator.h" 28 #include "modules/audio_coding/neteq/sync_buffer.h" 29 #include "rtc_base/checks.h" 30 #include "rtc_base/numerics/safe_conversions.h" 31 32 namespace webrtc { 33 34 Expand::Expand(BackgroundNoise* background_noise, 35 SyncBuffer* sync_buffer, 36 RandomVector* random_vector, 37 StatisticsCalculator* statistics, 38 int fs, 39 size_t num_channels) 40 : random_vector_(random_vector), 41 sync_buffer_(sync_buffer), 42 first_expand_(true), 43 fs_hz_(fs), 44 num_channels_(num_channels), 45 consecutive_expands_(0), 46 background_noise_(background_noise), 47 statistics_(statistics), 48 overlap_length_(5 * fs / 8000), 49 lag_index_direction_(0), 50 current_lag_index_(0), 51 stop_muting_(false), 52 expand_duration_samples_(0), 53 channel_parameters_(new ChannelParameters[num_channels_]) { 54 RTC_DCHECK(fs == 8000 || fs == 16000 || fs == 32000 || fs == 48000); 55 RTC_DCHECK_LE(fs, 56 static_cast<int>(kMaxSampleRate)); // Should not be possible. 57 RTC_DCHECK_GT(num_channels_, 0); 58 memset(expand_lags_, 0, sizeof(expand_lags_)); 59 Reset(); 60 } 61 62 Expand::~Expand() = default; 63 64 void Expand::Reset() { 65 first_expand_ = true; 66 consecutive_expands_ = 0; 67 max_lag_ = 0; 68 for (size_t ix = 0; ix < num_channels_; ++ix) { 69 channel_parameters_[ix].expand_vector0.Clear(); 70 channel_parameters_[ix].expand_vector1.Clear(); 71 } 72 } 73 74 int Expand::Process(AudioMultiVector* output) { 75 int16_t random_vector[kMaxSampleRate / 8000 * 120 + 30]; 76 int16_t scaled_random_vector[kMaxSampleRate / 8000 * 125]; 77 static const int kTempDataSize = 3600; 78 int16_t temp_data[kTempDataSize]; // TODO(hlundin) Remove this. 79 int16_t* voiced_vector_storage = temp_data; 80 int16_t* voiced_vector = &voiced_vector_storage[overlap_length_]; 81 static const size_t kNoiseLpcOrder = BackgroundNoise::kMaxLpcOrder; 82 int16_t unvoiced_array_memory[kNoiseLpcOrder + kMaxSampleRate / 8000 * 125]; 83 int16_t* unvoiced_vector = unvoiced_array_memory + kUnvoicedLpcOrder; 84 int16_t* noise_vector = unvoiced_array_memory + kNoiseLpcOrder; 85 86 int fs_mult = fs_hz_ / 8000; 87 88 if (first_expand_) { 89 // Perform initial setup if this is the first expansion since last reset. 90 AnalyzeSignal(random_vector); 91 first_expand_ = false; 92 expand_duration_samples_ = 0; 93 } else { 94 // This is not the first expansion, parameters are already estimated. 95 // Extract a noise segment. 96 size_t rand_length = max_lag_; 97 // This only applies to SWB where length could be larger than 256. 98 RTC_DCHECK_LE(rand_length, kMaxSampleRate / 8000 * 120 + 30); 99 GenerateRandomVector(2, rand_length, random_vector); 100 } 101 102 // Generate signal. 103 UpdateLagIndex(); 104 105 // Voiced part. 106 // Generate a weighted vector with the current lag. 107 const size_t expansion_vector_length = max_lag_ + overlap_length_; 108 const size_t current_lag = expand_lags_[current_lag_index_]; 109 // Copy lag+overlap data. 110 const size_t expansion_vector_position = 111 expansion_vector_length - current_lag - overlap_length_; 112 const size_t expansion_temp_length = current_lag + overlap_length_; 113 for (size_t channel_ix = 0; channel_ix < num_channels_; ++channel_ix) { 114 ChannelParameters& parameters = channel_parameters_[channel_ix]; 115 if (current_lag_index_ == 0) { 116 // Use only expand_vector0. 117 RTC_DCHECK_LE(expansion_vector_position + expansion_temp_length, 118 parameters.expand_vector0.Size()); 119 parameters.expand_vector0.CopyTo(expansion_temp_length, 120 expansion_vector_position, 121 voiced_vector_storage); 122 } else if (current_lag_index_ == 1) { 123 std::unique_ptr<int16_t[]> temp_0(new int16_t[expansion_temp_length]); 124 parameters.expand_vector0.CopyTo(expansion_temp_length, 125 expansion_vector_position, temp_0.get()); 126 std::unique_ptr<int16_t[]> temp_1(new int16_t[expansion_temp_length]); 127 parameters.expand_vector1.CopyTo(expansion_temp_length, 128 expansion_vector_position, temp_1.get()); 129 // Mix 3/4 of expand_vector0 with 1/4 of expand_vector1. 130 WebRtcSpl_ScaleAndAddVectorsWithRound(temp_0.get(), 3, temp_1.get(), 1, 2, 131 voiced_vector_storage, 132 expansion_temp_length); 133 } else if (current_lag_index_ == 2) { 134 // Mix 1/2 of expand_vector0 with 1/2 of expand_vector1. 135 RTC_DCHECK_LE(expansion_vector_position + expansion_temp_length, 136 parameters.expand_vector0.Size()); 137 RTC_DCHECK_LE(expansion_vector_position + expansion_temp_length, 138 parameters.expand_vector1.Size()); 139 140 std::unique_ptr<int16_t[]> temp_0(new int16_t[expansion_temp_length]); 141 parameters.expand_vector0.CopyTo(expansion_temp_length, 142 expansion_vector_position, temp_0.get()); 143 std::unique_ptr<int16_t[]> temp_1(new int16_t[expansion_temp_length]); 144 parameters.expand_vector1.CopyTo(expansion_temp_length, 145 expansion_vector_position, temp_1.get()); 146 WebRtcSpl_ScaleAndAddVectorsWithRound(temp_0.get(), 1, temp_1.get(), 1, 1, 147 voiced_vector_storage, 148 expansion_temp_length); 149 } 150 151 // Get tapering window parameters. Values are in Q15. 152 int16_t muting_window, muting_window_increment; 153 int16_t unmuting_window, unmuting_window_increment; 154 if (fs_hz_ == 8000) { 155 muting_window = DspHelper::kMuteFactorStart8kHz; 156 muting_window_increment = DspHelper::kMuteFactorIncrement8kHz; 157 unmuting_window = DspHelper::kUnmuteFactorStart8kHz; 158 unmuting_window_increment = DspHelper::kUnmuteFactorIncrement8kHz; 159 } else if (fs_hz_ == 16000) { 160 muting_window = DspHelper::kMuteFactorStart16kHz; 161 muting_window_increment = DspHelper::kMuteFactorIncrement16kHz; 162 unmuting_window = DspHelper::kUnmuteFactorStart16kHz; 163 unmuting_window_increment = DspHelper::kUnmuteFactorIncrement16kHz; 164 } else if (fs_hz_ == 32000) { 165 muting_window = DspHelper::kMuteFactorStart32kHz; 166 muting_window_increment = DspHelper::kMuteFactorIncrement32kHz; 167 unmuting_window = DspHelper::kUnmuteFactorStart32kHz; 168 unmuting_window_increment = DspHelper::kUnmuteFactorIncrement32kHz; 169 } else { // fs_ == 48000 170 muting_window = DspHelper::kMuteFactorStart48kHz; 171 muting_window_increment = DspHelper::kMuteFactorIncrement48kHz; 172 unmuting_window = DspHelper::kUnmuteFactorStart48kHz; 173 unmuting_window_increment = DspHelper::kUnmuteFactorIncrement48kHz; 174 } 175 176 // Smooth the expanded if it has not been muted to a low amplitude and 177 // `current_voice_mix_factor` is larger than 0.5. 178 if ((parameters.mute_factor > 819) && 179 (parameters.current_voice_mix_factor > 8192)) { 180 size_t start_ix = sync_buffer_->Size() - overlap_length_; 181 for (size_t i = 0; i < overlap_length_; i++) { 182 // Do overlap add between new vector and overlap. 183 (*sync_buffer_)[channel_ix][start_ix + i] = 184 (((*sync_buffer_)[channel_ix][start_ix + i] * muting_window) + 185 (((parameters.mute_factor * voiced_vector_storage[i]) >> 14) * 186 unmuting_window) + 187 16384) >> 188 15; 189 muting_window += muting_window_increment; 190 unmuting_window += unmuting_window_increment; 191 } 192 } else if (parameters.mute_factor == 0) { 193 // The expanded signal will consist of only comfort noise if 194 // mute_factor = 0. Set the output length to 15 ms for best noise 195 // production. 196 // TODO(hlundin): This has been disabled since the length of 197 // parameters.expand_vector0 and parameters.expand_vector1 no longer 198 // match with expand_lags_, causing invalid reads and writes. Is it a good 199 // idea to enable this again, and solve the vector size problem? 200 // max_lag_ = fs_mult * 120; 201 // expand_lags_[0] = fs_mult * 120; 202 // expand_lags_[1] = fs_mult * 120; 203 // expand_lags_[2] = fs_mult * 120; 204 } 205 206 // Unvoiced part. 207 // Filter `scaled_random_vector` through `ar_filter_`. 208 memcpy(unvoiced_vector - kUnvoicedLpcOrder, parameters.ar_filter_state, 209 sizeof(int16_t) * kUnvoicedLpcOrder); 210 int32_t add_constant = 0; 211 if (parameters.ar_gain_scale > 0) { 212 add_constant = 1 << (parameters.ar_gain_scale - 1); 213 } 214 WebRtcSpl_AffineTransformVector(scaled_random_vector, random_vector, 215 parameters.ar_gain, add_constant, 216 parameters.ar_gain_scale, current_lag); 217 WebRtcSpl_FilterARFastQ12(scaled_random_vector, unvoiced_vector, 218 parameters.ar_filter, kUnvoicedLpcOrder + 1, 219 current_lag); 220 memcpy(parameters.ar_filter_state, 221 &(unvoiced_vector[current_lag - kUnvoicedLpcOrder]), 222 sizeof(int16_t) * kUnvoicedLpcOrder); 223 224 // Combine voiced and unvoiced contributions. 225 226 // Set a suitable cross-fading slope. 227 // For lag = 228 // <= 31 * fs_mult => go from 1 to 0 in about 8 ms; 229 // (>= 31 .. <= 63) * fs_mult => go from 1 to 0 in about 16 ms; 230 // >= 64 * fs_mult => go from 1 to 0 in about 32 ms. 231 // temp_shift = getbits(max_lag_) - 5. 232 int temp_shift = 233 (31 - WebRtcSpl_NormW32(dchecked_cast<int32_t>(max_lag_))) - 5; 234 int16_t mix_factor_increment = 256 >> temp_shift; 235 if (stop_muting_) { 236 mix_factor_increment = 0; 237 } 238 239 // Create combined signal by shifting in more and more of unvoiced part. 240 temp_shift = 8 - temp_shift; // = getbits(mix_factor_increment). 241 size_t temp_length = 242 (parameters.current_voice_mix_factor - parameters.voice_mix_factor) >> 243 temp_shift; 244 temp_length = std::min(temp_length, current_lag); 245 DspHelper::CrossFade(voiced_vector, unvoiced_vector, temp_length, 246 ¶meters.current_voice_mix_factor, 247 mix_factor_increment, temp_data); 248 249 // End of cross-fading period was reached before end of expanded signal 250 // path. Mix the rest with a fixed mixing factor. 251 if (temp_length < current_lag) { 252 if (mix_factor_increment != 0) { 253 parameters.current_voice_mix_factor = parameters.voice_mix_factor; 254 } 255 int16_t temp_scale = 16384 - parameters.current_voice_mix_factor; 256 WebRtcSpl_ScaleAndAddVectorsWithRound( 257 voiced_vector + temp_length, parameters.current_voice_mix_factor, 258 unvoiced_vector + temp_length, temp_scale, 14, 259 temp_data + temp_length, current_lag - temp_length); 260 } 261 262 // Select muting slope depending on how many consecutive expands we have 263 // done. 264 if (consecutive_expands_ == 3) { 265 // Let the mute factor decrease from 1.0 to 0.95 in 6.25 ms. 266 // mute_slope = 0.0010 / fs_mult in Q20. 267 parameters.mute_slope = std::max(parameters.mute_slope, 1049 / fs_mult); 268 } 269 if (consecutive_expands_ == 7) { 270 // Let the mute factor decrease from 1.0 to 0.90 in 6.25 ms. 271 // mute_slope = 0.0020 / fs_mult in Q20. 272 parameters.mute_slope = std::max(parameters.mute_slope, 2097 / fs_mult); 273 } 274 275 // Mute segment according to slope value. 276 if ((consecutive_expands_ != 0) || !parameters.onset) { 277 // Mute to the previous level, then continue with the muting. 278 WebRtcSpl_AffineTransformVector( 279 temp_data, temp_data, parameters.mute_factor, 8192, 14, current_lag); 280 281 if (!stop_muting_) { 282 DspHelper::MuteSignal(temp_data, parameters.mute_slope, current_lag); 283 284 // Shift by 6 to go from Q20 to Q14. 285 // TODO(hlundin): Adding 8192 before shifting 6 steps seems wrong. 286 // Legacy. 287 int16_t gain = static_cast<int16_t>( 288 16384 - (((current_lag * parameters.mute_slope) + 8192) >> 6)); 289 gain = ((gain * parameters.mute_factor) + 8192) >> 14; 290 291 // Guard against getting stuck with very small (but sometimes audible) 292 // gain. 293 if ((consecutive_expands_ > 3) && (gain >= parameters.mute_factor)) { 294 parameters.mute_factor = 0; 295 } else { 296 parameters.mute_factor = gain; 297 } 298 } 299 } 300 301 // Background noise part. 302 background_noise_->GenerateBackgroundNoise( 303 random_vector, channel_ix, channel_parameters_[channel_ix].mute_slope, 304 TooManyExpands(), current_lag, unvoiced_array_memory); 305 306 // Add background noise to the combined voiced-unvoiced signal. 307 for (size_t i = 0; i < current_lag; i++) { 308 temp_data[i] = temp_data[i] + noise_vector[i]; 309 } 310 if (channel_ix == 0) { 311 output->AssertSize(current_lag); 312 } else { 313 RTC_DCHECK_EQ(output->Size(), current_lag); 314 } 315 (*output)[channel_ix].OverwriteAt(temp_data, current_lag, 0); 316 } 317 318 // Increase call number and cap it. 319 consecutive_expands_ = consecutive_expands_ >= kMaxConsecutiveExpands 320 ? kMaxConsecutiveExpands 321 : consecutive_expands_ + 1; 322 expand_duration_samples_ += output->Size(); 323 // Clamp the duration counter at 2 seconds. 324 expand_duration_samples_ = 325 std::min(expand_duration_samples_, dchecked_cast<size_t>(fs_hz_ * 2)); 326 return 0; 327 } 328 329 void Expand::SetParametersForNormalAfterExpand() { 330 current_lag_index_ = 0; 331 lag_index_direction_ = 0; 332 stop_muting_ = true; // Do not mute signal any more. 333 statistics_->LogDelayedPacketOutageEvent(expand_duration_samples_, fs_hz_); 334 statistics_->EndExpandEvent(fs_hz_); 335 } 336 337 void Expand::SetParametersForMergeAfterExpand() { 338 current_lag_index_ = -1; /* out of the 3 possible ones */ 339 lag_index_direction_ = 1; /* make sure we get the "optimal" lag */ 340 stop_muting_ = true; 341 statistics_->EndExpandEvent(fs_hz_); 342 } 343 344 bool Expand::Muted() const { 345 if (first_expand_ || stop_muting_) 346 return false; 347 RTC_DCHECK(channel_parameters_); 348 for (size_t ch = 0; ch < num_channels_; ++ch) { 349 if (channel_parameters_[ch].mute_factor != 0) 350 return false; 351 } 352 return true; 353 } 354 355 size_t Expand::overlap_length() const { 356 return overlap_length_; 357 } 358 359 void Expand::InitializeForAnExpandPeriod() { 360 lag_index_direction_ = 1; 361 current_lag_index_ = -1; 362 stop_muting_ = false; 363 random_vector_->set_seed_increment(1); 364 consecutive_expands_ = 0; 365 for (size_t ix = 0; ix < num_channels_; ++ix) { 366 channel_parameters_[ix].current_voice_mix_factor = 16384; // 1.0 in Q14. 367 channel_parameters_[ix].mute_factor = 16384; // 1.0 in Q14. 368 // Start with 0 gain for background noise. 369 background_noise_->SetMuteFactor(ix, 0); 370 } 371 } 372 373 bool Expand::TooManyExpands() { 374 return consecutive_expands_ >= kMaxConsecutiveExpands; 375 } 376 377 void Expand::AnalyzeSignal(int16_t* random_vector) { 378 int32_t auto_correlation[kUnvoicedLpcOrder + 1]; 379 int16_t reflection_coeff[kUnvoicedLpcOrder]; 380 int16_t correlation_vector[kMaxSampleRate / 8000 * 102]; 381 size_t best_correlation_index[kNumCorrelationCandidates]; 382 int16_t best_correlation[kNumCorrelationCandidates]; 383 size_t best_distortion_index[kNumCorrelationCandidates]; 384 int16_t best_distortion[kNumCorrelationCandidates]; 385 int32_t correlation_vector2[(99 * kMaxSampleRate / 8000) + 1]; 386 int32_t best_distortion_w32[kNumCorrelationCandidates]; 387 static const size_t kNoiseLpcOrder = BackgroundNoise::kMaxLpcOrder; 388 int16_t unvoiced_array_memory[kNoiseLpcOrder + kMaxSampleRate / 8000 * 125]; 389 int16_t* unvoiced_vector = unvoiced_array_memory + kUnvoicedLpcOrder; 390 391 int fs_mult = fs_hz_ / 8000; 392 393 // Pre-calculate common multiplications with fs_mult. 394 size_t fs_mult_4 = static_cast<size_t>(fs_mult * 4); 395 size_t fs_mult_20 = static_cast<size_t>(fs_mult * 20); 396 size_t fs_mult_120 = static_cast<size_t>(fs_mult * 120); 397 size_t fs_mult_dist_len = fs_mult * kDistortionLength; 398 size_t fs_mult_lpc_analysis_len = fs_mult * kLpcAnalysisLength; 399 400 const size_t signal_length = static_cast<size_t>(256 * fs_mult); 401 402 const size_t audio_history_position = sync_buffer_->Size() - signal_length; 403 std::unique_ptr<int16_t[]> audio_history(new int16_t[signal_length]); 404 (*sync_buffer_)[0].CopyTo(signal_length, audio_history_position, 405 audio_history.get()); 406 407 // Initialize. 408 InitializeForAnExpandPeriod(); 409 410 // Calculate correlation in downsampled domain (4 kHz sample rate). 411 size_t correlation_length = 51; // TODO(hlundin): Legacy bit-exactness. 412 // If it is decided to break bit-exactness `correlation_length` should be 413 // initialized to the return value of Correlation(). 414 Correlation(audio_history.get(), signal_length, correlation_vector); 415 416 // Find peaks in correlation vector. 417 DspHelper::PeakDetection(correlation_vector, correlation_length, 418 kNumCorrelationCandidates, fs_mult, 419 best_correlation_index, best_correlation); 420 421 // Adjust peak locations; cross-correlation lags start at 2.5 ms 422 // (20 * fs_mult samples). 423 best_correlation_index[0] += fs_mult_20; 424 best_correlation_index[1] += fs_mult_20; 425 best_correlation_index[2] += fs_mult_20; 426 427 // Calculate distortion around the `kNumCorrelationCandidates` best lags. 428 int distortion_scale = 0; 429 for (size_t i = 0; i < kNumCorrelationCandidates; i++) { 430 size_t min_index = 431 std::max(fs_mult_20, best_correlation_index[i] - fs_mult_4); 432 size_t max_index = 433 std::min(fs_mult_120 - 1, best_correlation_index[i] + fs_mult_4); 434 best_distortion_index[i] = DspHelper::MinDistortion( 435 &(audio_history[signal_length - fs_mult_dist_len]), min_index, 436 max_index, fs_mult_dist_len, &best_distortion_w32[i]); 437 distortion_scale = std::max(16 - WebRtcSpl_NormW32(best_distortion_w32[i]), 438 distortion_scale); 439 } 440 // Shift the distortion values to fit in 16 bits. 441 WebRtcSpl_VectorBitShiftW32ToW16(best_distortion, kNumCorrelationCandidates, 442 best_distortion_w32, distortion_scale); 443 444 // Find the maximizing index `i` of the cost function 445 // f[i] = best_correlation[i] / best_distortion[i]. 446 int32_t best_ratio = std::numeric_limits<int32_t>::min(); 447 size_t best_index = std::numeric_limits<size_t>::max(); 448 for (size_t i = 0; i < kNumCorrelationCandidates; ++i) { 449 int32_t ratio; 450 if (best_distortion[i] > 0) { 451 ratio = (best_correlation[i] * (1 << 16)) / best_distortion[i]; 452 } else if (best_correlation[i] == 0) { 453 ratio = 0; // No correlation set result to zero. 454 } else { 455 ratio = std::numeric_limits<int32_t>::max(); // Denominator is zero. 456 } 457 if (ratio > best_ratio) { 458 best_index = i; 459 best_ratio = ratio; 460 } 461 } 462 463 size_t distortion_lag = best_distortion_index[best_index]; 464 size_t correlation_lag = best_correlation_index[best_index]; 465 max_lag_ = std::max(distortion_lag, correlation_lag); 466 467 // Calculate the exact best correlation in the range between 468 // `correlation_lag` and `distortion_lag`. 469 correlation_length = std::max(std::min(distortion_lag + 10, fs_mult_120), 470 static_cast<size_t>(60 * fs_mult)); 471 472 size_t start_index = std::min(distortion_lag, correlation_lag); 473 size_t correlation_lags = static_cast<size_t>( 474 WEBRTC_SPL_ABS_W16((distortion_lag - correlation_lag)) + 1); 475 RTC_DCHECK_LE(correlation_lags, static_cast<size_t>(99 * fs_mult + 1)); 476 477 for (size_t channel_ix = 0; channel_ix < num_channels_; ++channel_ix) { 478 ChannelParameters& parameters = channel_parameters_[channel_ix]; 479 if (channel_ix > 0) { 480 // When channel_ix == 0, audio_history contains the correct audio. For the 481 // other cases, we will have to copy the correct channel into 482 // audio_history. 483 (*sync_buffer_)[channel_ix].CopyTo(signal_length, audio_history_position, 484 audio_history.get()); 485 } 486 487 // Calculate suitable scaling. 488 int16_t signal_max = WebRtcSpl_MaxAbsValueW16( 489 &audio_history[signal_length - correlation_length - start_index - 490 correlation_lags], 491 correlation_length + start_index + correlation_lags - 1); 492 int correlation_scale = 493 (31 - WebRtcSpl_NormW32(signal_max * signal_max)) + 494 (31 - WebRtcSpl_NormW32(static_cast<int32_t>(correlation_length))) - 31; 495 correlation_scale = std::max(0, correlation_scale); 496 497 // Calculate the correlation, store in `correlation_vector2`. 498 WebRtcSpl_CrossCorrelation( 499 correlation_vector2, 500 &(audio_history[signal_length - correlation_length]), 501 &(audio_history[signal_length - correlation_length - start_index]), 502 correlation_length, correlation_lags, correlation_scale, -1); 503 504 // Find maximizing index. 505 best_index = WebRtcSpl_MaxIndexW32(correlation_vector2, correlation_lags); 506 int32_t max_correlation = correlation_vector2[best_index]; 507 // Compensate index with start offset. 508 best_index = best_index + start_index; 509 510 // Calculate energies. 511 int32_t energy1 = WebRtcSpl_DotProductWithScale( 512 &(audio_history[signal_length - correlation_length]), 513 &(audio_history[signal_length - correlation_length]), 514 correlation_length, correlation_scale); 515 int32_t energy2 = WebRtcSpl_DotProductWithScale( 516 &(audio_history[signal_length - correlation_length - best_index]), 517 &(audio_history[signal_length - correlation_length - best_index]), 518 correlation_length, correlation_scale); 519 520 // Calculate the correlation coefficient between the two portions of the 521 // signal. 522 int32_t corr_coefficient; 523 if ((energy1 > 0) && (energy2 > 0)) { 524 int energy1_scale = std::max(16 - WebRtcSpl_NormW32(energy1), 0); 525 int energy2_scale = std::max(16 - WebRtcSpl_NormW32(energy2), 0); 526 // Make sure total scaling is even (to simplify scale factor after sqrt). 527 if ((energy1_scale + energy2_scale) & 1) { 528 // If sum is odd, add 1 to make it even. 529 energy1_scale += 1; 530 } 531 int32_t scaled_energy1 = energy1 >> energy1_scale; 532 int32_t scaled_energy2 = energy2 >> energy2_scale; 533 int16_t sqrt_energy_product = static_cast<int16_t>( 534 WebRtcSpl_SqrtFloor(scaled_energy1 * scaled_energy2)); 535 // Calculate max_correlation / sqrt(energy1 * energy2) in Q14. 536 int cc_shift = 14 - (energy1_scale + energy2_scale) / 2; 537 max_correlation = WEBRTC_SPL_SHIFT_W32(max_correlation, cc_shift); 538 corr_coefficient = 539 WebRtcSpl_DivW32W16(max_correlation, sqrt_energy_product); 540 // Cap at 1.0 in Q14. 541 corr_coefficient = std::min(16384, corr_coefficient); 542 } else { 543 corr_coefficient = 0; 544 } 545 546 // Extract the two vectors expand_vector0 and expand_vector1 from 547 // `audio_history`. 548 size_t expansion_length = max_lag_ + overlap_length_; 549 const int16_t* vector1 = &(audio_history[signal_length - expansion_length]); 550 const int16_t* vector2 = vector1 - distortion_lag; 551 // Normalize the second vector to the same energy as the first. 552 energy1 = WebRtcSpl_DotProductWithScale(vector1, vector1, expansion_length, 553 correlation_scale); 554 energy2 = WebRtcSpl_DotProductWithScale(vector2, vector2, expansion_length, 555 correlation_scale); 556 // Confirm that amplitude ratio sqrt(energy1 / energy2) is within 0.5 - 2.0, 557 // i.e., energy1 / energy2 is within 0.25 - 4. 558 int16_t amplitude_ratio; 559 if ((energy1 / 4 < energy2) && (energy1 > energy2 / 4)) { 560 // Energy constraint fulfilled. Use both vectors and scale them 561 // accordingly. 562 int32_t scaled_energy2 = std::max(16 - WebRtcSpl_NormW32(energy2), 0); 563 int32_t scaled_energy1 = scaled_energy2 - 13; 564 // Calculate scaled_energy1 / scaled_energy2 in Q13. 565 int32_t energy_ratio = 566 WebRtcSpl_DivW32W16(WEBRTC_SPL_SHIFT_W32(energy1, -scaled_energy1), 567 static_cast<int16_t>(energy2 >> scaled_energy2)); 568 // Calculate sqrt ratio in Q13 (sqrt of en1/en2 in Q26). 569 amplitude_ratio = 570 static_cast<int16_t>(WebRtcSpl_SqrtFloor(energy_ratio << 13)); 571 // Copy the two vectors and give them the same energy. 572 parameters.expand_vector0.Clear(); 573 parameters.expand_vector0.PushBack(vector1, expansion_length); 574 parameters.expand_vector1.Clear(); 575 if (parameters.expand_vector1.Size() < expansion_length) { 576 parameters.expand_vector1.Extend(expansion_length - 577 parameters.expand_vector1.Size()); 578 } 579 std::unique_ptr<int16_t[]> temp_1(new int16_t[expansion_length]); 580 WebRtcSpl_AffineTransformVector( 581 temp_1.get(), const_cast<int16_t*>(vector2), amplitude_ratio, 4096, 582 13, expansion_length); 583 parameters.expand_vector1.OverwriteAt(temp_1.get(), expansion_length, 0); 584 } else { 585 // Energy change constraint not fulfilled. Only use last vector. 586 parameters.expand_vector0.Clear(); 587 parameters.expand_vector0.PushBack(vector1, expansion_length); 588 // Copy from expand_vector0 to expand_vector1. 589 parameters.expand_vector0.CopyTo(¶meters.expand_vector1); 590 // Set the energy_ratio since it is used by muting slope. 591 if ((energy1 / 4 < energy2) || (energy2 == 0)) { 592 amplitude_ratio = 4096; // 0.5 in Q13. 593 } else { 594 amplitude_ratio = 16384; // 2.0 in Q13. 595 } 596 } 597 598 // Set the 3 lag values. 599 if (distortion_lag == correlation_lag) { 600 expand_lags_[0] = distortion_lag; 601 expand_lags_[1] = distortion_lag; 602 expand_lags_[2] = distortion_lag; 603 } else { 604 // `distortion_lag` and `correlation_lag` are not equal; use different 605 // combinations of the two. 606 // First lag is `distortion_lag` only. 607 expand_lags_[0] = distortion_lag; 608 // Second lag is the average of the two. 609 expand_lags_[1] = (distortion_lag + correlation_lag) / 2; 610 // Third lag is the average again, but rounding towards `correlation_lag`. 611 if (distortion_lag > correlation_lag) { 612 expand_lags_[2] = (distortion_lag + correlation_lag - 1) / 2; 613 } else { 614 expand_lags_[2] = (distortion_lag + correlation_lag + 1) / 2; 615 } 616 } 617 618 // Calculate the LPC and the gain of the filters. 619 620 // Calculate kUnvoicedLpcOrder + 1 lags of the auto-correlation function. 621 size_t temp_index = 622 signal_length - fs_mult_lpc_analysis_len - kUnvoicedLpcOrder; 623 // Copy signal to temporary vector to be able to pad with leading zeros. 624 int16_t* temp_signal = 625 new int16_t[fs_mult_lpc_analysis_len + kUnvoicedLpcOrder]; 626 memset(temp_signal, 0, 627 sizeof(int16_t) * (fs_mult_lpc_analysis_len + kUnvoicedLpcOrder)); 628 memcpy(&temp_signal[kUnvoicedLpcOrder], 629 &audio_history[temp_index + kUnvoicedLpcOrder], 630 sizeof(int16_t) * fs_mult_lpc_analysis_len); 631 CrossCorrelationWithAutoShift( 632 &temp_signal[kUnvoicedLpcOrder], &temp_signal[kUnvoicedLpcOrder], 633 fs_mult_lpc_analysis_len, kUnvoicedLpcOrder + 1, -1, auto_correlation); 634 delete[] temp_signal; 635 636 // Verify that variance is positive. 637 if (auto_correlation[0] > 0) { 638 // Estimate AR filter parameters using Levinson-Durbin algorithm; 639 // kUnvoicedLpcOrder + 1 filter coefficients. 640 int16_t stability = 641 WebRtcSpl_LevinsonDurbin(auto_correlation, parameters.ar_filter, 642 reflection_coeff, kUnvoicedLpcOrder); 643 644 // Keep filter parameters only if filter is stable. 645 if (stability != 1) { 646 // Set first coefficient to 4096 (1.0 in Q12). 647 parameters.ar_filter[0] = 4096; 648 // Set remaining `kUnvoicedLpcOrder` coefficients to zero. 649 WebRtcSpl_MemSetW16(parameters.ar_filter + 1, 0, kUnvoicedLpcOrder); 650 } 651 } 652 653 if (channel_ix == 0) { 654 // Extract a noise segment. 655 size_t noise_length; 656 if (distortion_lag < 40) { 657 noise_length = 2 * distortion_lag + 30; 658 } else { 659 noise_length = distortion_lag + 30; 660 } 661 if (noise_length <= RandomVector::kRandomTableSize) { 662 memcpy(random_vector, RandomVector::kRandomTable, 663 sizeof(int16_t) * noise_length); 664 } else { 665 // Only applies to SWB where length could be larger than 666 // `kRandomTableSize`. 667 memcpy(random_vector, RandomVector::kRandomTable, 668 sizeof(int16_t) * RandomVector::kRandomTableSize); 669 RTC_DCHECK_LE(noise_length, kMaxSampleRate / 8000 * 120 + 30); 670 random_vector_->IncreaseSeedIncrement(2); 671 random_vector_->Generate( 672 noise_length - RandomVector::kRandomTableSize, 673 &random_vector[RandomVector::kRandomTableSize]); 674 } 675 } 676 677 // Set up state vector and calculate scale factor for unvoiced filtering. 678 memcpy(parameters.ar_filter_state, 679 &(audio_history[signal_length - kUnvoicedLpcOrder]), 680 sizeof(int16_t) * kUnvoicedLpcOrder); 681 memcpy(unvoiced_vector - kUnvoicedLpcOrder, 682 &(audio_history[signal_length - 128 - kUnvoicedLpcOrder]), 683 sizeof(int16_t) * kUnvoicedLpcOrder); 684 WebRtcSpl_FilterMAFastQ12(&audio_history[signal_length - 128], 685 unvoiced_vector, parameters.ar_filter, 686 kUnvoicedLpcOrder + 1, 128); 687 const int unvoiced_max_abs = [&] { 688 const int16_t max_abs = WebRtcSpl_MaxAbsValueW16(unvoiced_vector, 128); 689 // Since WebRtcSpl_MaxAbsValueW16 returns 2^15 - 1 when the input contains 690 // -2^15, we have to conservatively bump the return value by 1 691 // if it is 2^15 - 1. 692 return max_abs == WEBRTC_SPL_WORD16_MAX ? max_abs + 1 : max_abs; 693 }(); 694 // Pick the smallest n such that 2^n > unvoiced_max_abs; then the maximum 695 // value of the dot product is less than 2^7 * 2^(2*n) = 2^(2*n + 7), so to 696 // prevent overflows we want 2n + 7 <= 31, which means we should shift by 697 // 2n + 7 - 31 bits, if this value is greater than zero. 698 int unvoiced_prescale = 699 std::max(0, 2 * WebRtcSpl_GetSizeInBits(unvoiced_max_abs) - 24); 700 701 int32_t unvoiced_energy = WebRtcSpl_DotProductWithScale( 702 unvoiced_vector, unvoiced_vector, 128, unvoiced_prescale); 703 704 // Normalize `unvoiced_energy` to 28 or 29 bits to preserve sqrt() accuracy. 705 int16_t unvoiced_scale = WebRtcSpl_NormW32(unvoiced_energy) - 3; 706 // Make sure we do an odd number of shifts since we already have 7 shifts 707 // from dividing with 128 earlier. This will make the total scale factor 708 // even, which is suitable for the sqrt. 709 unvoiced_scale += ((unvoiced_scale & 0x1) ^ 0x1); 710 unvoiced_energy = WEBRTC_SPL_SHIFT_W32(unvoiced_energy, unvoiced_scale); 711 int16_t unvoiced_gain = 712 static_cast<int16_t>(WebRtcSpl_SqrtFloor(unvoiced_energy)); 713 parameters.ar_gain_scale = 714 13 + (unvoiced_scale + 7 - unvoiced_prescale) / 2; 715 parameters.ar_gain = unvoiced_gain; 716 717 // Calculate voice_mix_factor from corr_coefficient. 718 // Let x = corr_coefficient. Then, we compute: 719 // if (x > 0.48) 720 // voice_mix_factor = (-5179 + 19931x - 16422x^2 + 5776x^3) / 4096; 721 // else 722 // voice_mix_factor = 0; 723 if (corr_coefficient > 7875) { 724 int16_t x1, x2, x3; 725 // `corr_coefficient` is in Q14. 726 x1 = static_cast<int16_t>(corr_coefficient); 727 x2 = (x1 * x1) >> 14; // Shift 14 to keep result in Q14. 728 x3 = (x1 * x2) >> 14; 729 static const int kCoefficients[4] = {-5179, 19931, -16422, 5776}; 730 int32_t temp_sum = kCoefficients[0] * 16384; 731 temp_sum += kCoefficients[1] * x1; 732 temp_sum += kCoefficients[2] * x2; 733 temp_sum += kCoefficients[3] * x3; 734 parameters.voice_mix_factor = 735 static_cast<int16_t>(std::min(temp_sum / 4096, 16384)); 736 parameters.voice_mix_factor = 737 std::max(parameters.voice_mix_factor, static_cast<int16_t>(0)); 738 } else { 739 parameters.voice_mix_factor = 0; 740 } 741 742 // Calculate muting slope. Reuse value from earlier scaling of 743 // `expand_vector0` and `expand_vector1`. 744 int16_t slope = amplitude_ratio; 745 if (slope > 12288) { 746 // slope > 1.5. 747 // Calculate (1 - (1 / slope)) / distortion_lag = 748 // (slope - 1) / (distortion_lag * slope). 749 // `slope` is in Q13, so 1 corresponds to 8192. Shift up to Q25 before 750 // the division. 751 // Shift the denominator from Q13 to Q5 before the division. The result of 752 // the division will then be in Q20. 753 int16_t denom = saturated_cast<int16_t>((distortion_lag * slope) >> 8); 754 int temp_ratio = WebRtcSpl_DivW32W16((slope - 8192) << 12, denom); 755 if (slope > 14746) { 756 // slope > 1.8. 757 // Divide by 2, with proper rounding. 758 parameters.mute_slope = (temp_ratio + 1) / 2; 759 } else { 760 // Divide by 8, with proper rounding. 761 parameters.mute_slope = (temp_ratio + 4) / 8; 762 } 763 parameters.onset = true; 764 } else { 765 // Calculate (1 - slope) / distortion_lag. 766 // Shift `slope` by 7 to Q20 before the division. The result is in Q20. 767 parameters.mute_slope = WebRtcSpl_DivW32W16( 768 (8192 - slope) * 128, static_cast<int16_t>(distortion_lag)); 769 if (parameters.voice_mix_factor <= 13107) { 770 // Make sure the mute factor decreases from 1.0 to 0.9 in no more than 771 // 6.25 ms. 772 // mute_slope >= 0.005 / fs_mult in Q20. 773 parameters.mute_slope = std::max(5243 / fs_mult, parameters.mute_slope); 774 } else if (slope > 8028) { 775 parameters.mute_slope = 0; 776 } 777 parameters.onset = false; 778 } 779 } 780 } 781 782 Expand::ChannelParameters::ChannelParameters() 783 : mute_factor(16384), 784 ar_gain(0), 785 ar_gain_scale(0), 786 voice_mix_factor(0), 787 current_voice_mix_factor(0), 788 onset(false), 789 mute_slope(0) { 790 memset(ar_filter, 0, sizeof(ar_filter)); 791 memset(ar_filter_state, 0, sizeof(ar_filter_state)); 792 } 793 794 void Expand::Correlation(const int16_t* input, 795 size_t input_length, 796 int16_t* output) const { 797 // Set parameters depending on sample rate. 798 const int16_t* filter_coefficients; 799 size_t num_coefficients; 800 int16_t downsampling_factor; 801 if (fs_hz_ == 8000) { 802 num_coefficients = 3; 803 downsampling_factor = 2; 804 filter_coefficients = DspHelper::kDownsample8kHzTbl; 805 } else if (fs_hz_ == 16000) { 806 num_coefficients = 5; 807 downsampling_factor = 4; 808 filter_coefficients = DspHelper::kDownsample16kHzTbl; 809 } else if (fs_hz_ == 32000) { 810 num_coefficients = 7; 811 downsampling_factor = 8; 812 filter_coefficients = DspHelper::kDownsample32kHzTbl; 813 } else { // fs_hz_ == 48000. 814 num_coefficients = 7; 815 downsampling_factor = 12; 816 filter_coefficients = DspHelper::kDownsample48kHzTbl; 817 } 818 819 // Correlate from lag 10 to lag 60 in downsampled domain. 820 // (Corresponds to 20-120 for narrow-band, 40-240 for wide-band, and so on.) 821 static const size_t kCorrelationStartLag = 10; 822 static const size_t kNumCorrelationLags = 54; 823 static const size_t kCorrelationLength = 60; 824 // Downsample to 4 kHz sample rate. 825 static const size_t kDownsampledLength = 826 kCorrelationStartLag + kNumCorrelationLags + kCorrelationLength; 827 int16_t downsampled_input[kDownsampledLength]; 828 static const size_t kFilterDelay = 0; 829 WebRtcSpl_DownsampleFast( 830 input + input_length - kDownsampledLength * downsampling_factor, 831 kDownsampledLength * downsampling_factor, downsampled_input, 832 kDownsampledLength, filter_coefficients, num_coefficients, 833 downsampling_factor, kFilterDelay); 834 835 // Normalize `downsampled_input` to using all 16 bits. 836 int16_t max_value = 837 WebRtcSpl_MaxAbsValueW16(downsampled_input, kDownsampledLength); 838 int16_t norm_shift = 16 - WebRtcSpl_NormW32(max_value); 839 WebRtcSpl_VectorBitShiftW16(downsampled_input, kDownsampledLength, 840 downsampled_input, norm_shift); 841 842 int32_t correlation[kNumCorrelationLags]; 843 CrossCorrelationWithAutoShift( 844 &downsampled_input[kDownsampledLength - kCorrelationLength], 845 &downsampled_input[kDownsampledLength - kCorrelationLength - 846 kCorrelationStartLag], 847 kCorrelationLength, kNumCorrelationLags, -1, correlation); 848 849 // Normalize and move data from 32-bit to 16-bit vector. 850 int32_t max_correlation = 851 WebRtcSpl_MaxAbsValueW32(correlation, kNumCorrelationLags); 852 int16_t norm_shift2 = static_cast<int16_t>( 853 std::max(18 - WebRtcSpl_NormW32(max_correlation), 0)); 854 WebRtcSpl_VectorBitShiftW32ToW16(output, kNumCorrelationLags, correlation, 855 norm_shift2); 856 } 857 858 void Expand::UpdateLagIndex() { 859 current_lag_index_ = current_lag_index_ + lag_index_direction_; 860 // Change direction if needed. 861 if (current_lag_index_ <= 0) { 862 lag_index_direction_ = 1; 863 } 864 if (current_lag_index_ >= kNumLags - 1) { 865 lag_index_direction_ = -1; 866 } 867 } 868 869 Expand* ExpandFactory::Create(BackgroundNoise* background_noise, 870 SyncBuffer* sync_buffer, 871 RandomVector* random_vector, 872 StatisticsCalculator* statistics, 873 int fs, 874 size_t num_channels) const { 875 return new Expand(background_noise, sync_buffer, random_vector, statistics, 876 fs, num_channels); 877 } 878 879 void Expand::GenerateRandomVector(int16_t seed_increment, 880 size_t length, 881 int16_t* random_vector) { 882 // TODO(turajs): According to hlundin The loop should not be needed. Should be 883 // just as good to generate all of the vector in one call. 884 size_t samples_generated = 0; 885 const size_t kMaxRandSamples = RandomVector::kRandomTableSize; 886 while (samples_generated < length) { 887 size_t rand_length = std::min(length - samples_generated, kMaxRandSamples); 888 random_vector_->IncreaseSeedIncrement(seed_increment); 889 random_vector_->Generate(rand_length, &random_vector[samples_generated]); 890 samples_generated += rand_length; 891 } 892 } 893 894 } // namespace webrtc