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