• 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