modules/audio_coding/neteq/dsp_helper.cc

/*
 *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
 *
 *  Use of this source code is governed by a BSD-style license
 *  that can be found in the LICENSE file in the root of the source
 *  tree. An additional intellectual property rights grant can be found
 *  in the file PATENTS.  All contributing project authors may
 *  be found in the AUTHORS file in the root of the source tree.
 */

#include "modules/audio_coding/neteq/dsp_helper.h"

#include <assert.h>
#include <string.h>  // Access to memset.

#include <algorithm>  // Access to min, max.

#include "common_audio/signal_processing/include/signal_processing_library.h"

namespace webrtc {

// Table of constants used in method DspHelper::ParabolicFit().
const int16_t DspHelper::kParabolaCoefficients[17][3] = {
    {120, 32, 64},   {140, 44, 75},   {150, 50, 80},   {160, 57, 85},
    {180, 72, 96},   {200, 89, 107},  {210, 98, 112},  {220, 108, 117},
    {240, 128, 128}, {260, 150, 139}, {270, 162, 144}, {280, 174, 149},
    {300, 200, 160}, {320, 228, 171}, {330, 242, 176}, {340, 257, 181},
    {360, 288, 192}};

// Filter coefficients used when downsampling from the indicated sample rates
// (8, 16, 32, 48 kHz) to 4 kHz. Coefficients are in Q12. The corresponding Q0
// values are provided in the comments before each array.

// Q0 values: {0.3, 0.4, 0.3}.
const int16_t DspHelper::kDownsample8kHzTbl[3] = {1229, 1638, 1229};

// Q0 values: {0.15, 0.2, 0.3, 0.2, 0.15}.
const int16_t DspHelper::kDownsample16kHzTbl[5] = {614, 819, 1229, 819, 614};

// Q0 values: {0.1425, 0.1251, 0.1525, 0.1628, 0.1525, 0.1251, 0.1425}.
const int16_t DspHelper::kDownsample32kHzTbl[7] = {584, 512, 625, 667,
                                                   625, 512, 584};

// Q0 values: {0.2487, 0.0952, 0.1042, 0.1074, 0.1042, 0.0952, 0.2487}.
const int16_t DspHelper::kDownsample48kHzTbl[7] = {1019, 390, 427, 440,
                                                   427,  390, 1019};

int DspHelper::RampSignal(const int16_t* input,
                          size_t length,
                          int factor,
                          int increment,
                          int16_t* output) {
  int factor_q20 = (factor << 6) + 32;
  // TODO(hlundin): Add 32 to factor_q20 when converting back to Q14?
  for (size_t i = 0; i < length; ++i) {
    output[i] = (factor * input[i] + 8192) >> 14;
    factor_q20 += increment;
    factor_q20 = std::max(factor_q20, 0);  // Never go negative.
    factor = std::min(factor_q20 >> 6, 16384);
  }
  return factor;
}

int DspHelper::RampSignal(int16_t* signal,
                          size_t length,
                          int factor,
                          int increment) {
  return RampSignal(signal, length, factor, increment, signal);
}

int DspHelper::RampSignal(AudioVector* signal,
                          size_t start_index,
                          size_t length,
                          int factor,
                          int increment) {
  int factor_q20 = (factor << 6) + 32;
  // TODO(hlundin): Add 32 to factor_q20 when converting back to Q14?
  for (size_t i = start_index; i < start_index + length; ++i) {
    (*signal)[i] = (factor * (*signal)[i] + 8192) >> 14;
    factor_q20 += increment;
    factor_q20 = std::max(factor_q20, 0);  // Never go negative.
    factor = std::min(factor_q20 >> 6, 16384);
  }
  return factor;
}

int DspHelper::RampSignal(AudioMultiVector* signal,
                          size_t start_index,
                          size_t length,
                          int factor,
                          int increment) {
  assert(start_index + length <= signal->Size());
  if (start_index + length > signal->Size()) {
    // Wrong parameters. Do nothing and return the scale factor unaltered.
    return factor;
  }
  int end_factor = 0;
  // Loop over the channels, starting at the same |factor| each time.
  for (size_t channel = 0; channel < signal->Channels(); ++channel) {
    end_factor =
        RampSignal(&(*signal)[channel], start_index, length, factor, increment);
  }
  return end_factor;
}

void DspHelper::PeakDetection(int16_t* data,
                              size_t data_length,
                              size_t num_peaks,
                              int fs_mult,
                              size_t* peak_index,
                              int16_t* peak_value) {
  size_t min_index = 0;
  size_t max_index = 0;

  for (size_t i = 0; i <= num_peaks - 1; i++) {
    if (num_peaks == 1) {
      // Single peak.  The parabola fit assumes that an extra point is
      // available; worst case it gets a zero on the high end of the signal.
      // TODO(hlundin): This can potentially get much worse. It breaks the
      // API contract, that the length of |data| is |data_length|.
      data_length++;
    }

    peak_index[i] = WebRtcSpl_MaxIndexW16(data, data_length - 1);

    if (i != num_peaks - 1) {
      min_index = (peak_index[i] > 2) ? (peak_index[i] - 2) : 0;
      max_index = std::min(data_length - 1, peak_index[i] + 2);
    }

    if ((peak_index[i] != 0) && (peak_index[i] != (data_length - 2))) {
      ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
                   &peak_value[i]);
    } else {
      if (peak_index[i] == data_length - 2) {
        if (data[peak_index[i]] > data[peak_index[i] + 1]) {
          ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
                       &peak_value[i]);
        } else if (data[peak_index[i]] <= data[peak_index[i] + 1]) {
          // Linear approximation.
          peak_value[i] = (data[peak_index[i]] + data[peak_index[i] + 1]) >> 1;
          peak_index[i] = (peak_index[i] * 2 + 1) * fs_mult;
        }
      } else {
        peak_value[i] = data[peak_index[i]];
        peak_index[i] = peak_index[i] * 2 * fs_mult;
      }
    }

    if (i != num_peaks - 1) {
      memset(&data[min_index], 0,
             sizeof(data[0]) * (max_index - min_index + 1));
    }
  }
}

void DspHelper::ParabolicFit(int16_t* signal_points,
                             int fs_mult,
                             size_t* peak_index,
                             int16_t* peak_value) {
  uint16_t fit_index[13];
  if (fs_mult == 1) {
    fit_index[0] = 0;
    fit_index[1] = 8;
    fit_index[2] = 16;
  } else if (fs_mult == 2) {
    fit_index[0] = 0;
    fit_index[1] = 4;
    fit_index[2] = 8;
    fit_index[3] = 12;
    fit_index[4] = 16;
  } else if (fs_mult == 4) {
    fit_index[0] = 0;
    fit_index[1] = 2;
    fit_index[2] = 4;
    fit_index[3] = 6;
    fit_index[4] = 8;
    fit_index[5] = 10;
    fit_index[6] = 12;
    fit_index[7] = 14;
    fit_index[8] = 16;
  } else {
    fit_index[0] = 0;
    fit_index[1] = 1;
    fit_index[2] = 3;
    fit_index[3] = 4;
    fit_index[4] = 5;
    fit_index[5] = 7;
    fit_index[6] = 8;
    fit_index[7] = 9;
    fit_index[8] = 11;
    fit_index[9] = 12;
    fit_index[10] = 13;
    fit_index[11] = 15;
    fit_index[12] = 16;
  }

  //  num = -3 * signal_points[0] + 4 * signal_points[1] - signal_points[2];
  //  den =      signal_points[0] - 2 * signal_points[1] + signal_points[2];
  int32_t num =
      (signal_points[0] * -3) + (signal_points[1] * 4) - signal_points[2];
  int32_t den = signal_points[0] + (signal_points[1] * -2) + signal_points[2];
  int32_t temp = num * 120;
  int flag = 1;
  int16_t stp = kParabolaCoefficients[fit_index[fs_mult]][0] -
                kParabolaCoefficients[fit_index[fs_mult - 1]][0];
  int16_t strt = (kParabolaCoefficients[fit_index[fs_mult]][0] +
                  kParabolaCoefficients[fit_index[fs_mult - 1]][0]) /
                 2;
  int16_t lmt;
  if (temp < -den * strt) {
    lmt = strt - stp;
    while (flag) {
      if ((flag == fs_mult) || (temp > -den * lmt)) {
        *peak_value =
            (den * kParabolaCoefficients[fit_index[fs_mult - flag]][1] +
             num * kParabolaCoefficients[fit_index[fs_mult - flag]][2] +
             signal_points[0] * 256) /
            256;
        *peak_index = *peak_index * 2 * fs_mult - flag;
        flag = 0;
      } else {
        flag++;
        lmt -= stp;
      }
    }
  } else if (temp > -den * (strt + stp)) {
    lmt = strt + 2 * stp;
    while (flag) {
      if ((flag == fs_mult) || (temp < -den * lmt)) {
        int32_t temp_term_1 =
            den * kParabolaCoefficients[fit_index[fs_mult + flag]][1];
        int32_t temp_term_2 =
            num * kParabolaCoefficients[fit_index[fs_mult + flag]][2];
        int32_t temp_term_3 = signal_points[0] * 256;
        *peak_value = (temp_term_1 + temp_term_2 + temp_term_3) / 256;
        *peak_index = *peak_index * 2 * fs_mult + flag;
        flag = 0;
      } else {
        flag++;
        lmt += stp;
      }
    }
  } else {
    *peak_value = signal_points[1];
    *peak_index = *peak_index * 2 * fs_mult;
  }
}

size_t DspHelper::MinDistortion(const int16_t* signal,
                                size_t min_lag,
                                size_t max_lag,
                                size_t length,
                                int32_t* distortion_value) {
  size_t best_index = 0;
  int32_t min_distortion = WEBRTC_SPL_WORD32_MAX;
  for (size_t i = min_lag; i <= max_lag; i++) {
    int32_t sum_diff = 0;
    const int16_t* data1 = signal;
    const int16_t* data2 = signal - i;
    for (size_t j = 0; j < length; j++) {
      sum_diff += WEBRTC_SPL_ABS_W32(data1[j] - data2[j]);
    }
    // Compare with previous minimum.
    if (sum_diff < min_distortion) {
      min_distortion = sum_diff;
      best_index = i;
    }
  }
  *distortion_value = min_distortion;
  return best_index;
}

void DspHelper::CrossFade(const int16_t* input1,
                          const int16_t* input2,
                          size_t length,
                          int16_t* mix_factor,
                          int16_t factor_decrement,
                          int16_t* output) {
  int16_t factor = *mix_factor;
  int16_t complement_factor = 16384 - factor;
  for (size_t i = 0; i < length; i++) {
    output[i] =
        (factor * input1[i] + complement_factor * input2[i] + 8192) >> 14;
    factor -= factor_decrement;
    complement_factor += factor_decrement;
  }
  *mix_factor = factor;
}

void DspHelper::UnmuteSignal(const int16_t* input,
                             size_t length,
                             int16_t* factor,
                             int increment,
                             int16_t* output) {
  uint16_t factor_16b = *factor;
  int32_t factor_32b = (static_cast<int32_t>(factor_16b) << 6) + 32;
  for (size_t i = 0; i < length; i++) {
    output[i] = (factor_16b * input[i] + 8192) >> 14;
    factor_32b = std::max(factor_32b + increment, 0);
    factor_16b = std::min(16384, factor_32b >> 6);
  }
  *factor = factor_16b;
}

void DspHelper::MuteSignal(int16_t* signal, int mute_slope, size_t length) {
  int32_t factor = (16384 << 6) + 32;
  for (size_t i = 0; i < length; i++) {
    signal[i] = ((factor >> 6) * signal[i] + 8192) >> 14;
    factor -= mute_slope;
  }
}

int DspHelper::DownsampleTo4kHz(const int16_t* input,
                                size_t input_length,
                                size_t output_length,
                                int input_rate_hz,
                                bool compensate_delay,
                                int16_t* output) {
  // Set filter parameters depending on input frequency.
  // NOTE: The phase delay values are wrong compared to the true phase delay
  // of the filters. However, the error is preserved (through the +1 term) for
  // consistency.
  const int16_t* filter_coefficients;  // Filter coefficients.
  size_t filter_length;                // Number of coefficients.
  size_t filter_delay;                 // Phase delay in samples.
  int16_t factor;                      // Conversion rate (inFsHz / 8000).
  switch (input_rate_hz) {
    case 8000: {
      filter_length = 3;
      factor = 2;
      filter_coefficients = kDownsample8kHzTbl;
      filter_delay = 1 + 1;
      break;
    }
    case 16000: {
      filter_length = 5;
      factor = 4;
      filter_coefficients = kDownsample16kHzTbl;
      filter_delay = 2 + 1;
      break;
    }
    case 32000: {
      filter_length = 7;
      factor = 8;
      filter_coefficients = kDownsample32kHzTbl;
      filter_delay = 3 + 1;
      break;
    }
    case 48000: {
      filter_length = 7;
      factor = 12;
      filter_coefficients = kDownsample48kHzTbl;
      filter_delay = 3 + 1;
      break;
    }
    default: {
      assert(false);
      return -1;
    }
  }

  if (!compensate_delay) {
    // Disregard delay compensation.
    filter_delay = 0;
  }

  // Returns -1 if input signal is too short; 0 otherwise.
  return WebRtcSpl_DownsampleFast(
      &input[filter_length - 1], input_length - filter_length + 1, output,
      output_length, filter_coefficients, filter_length, factor, filter_delay);
}

}  // namespace webrtc