Home | History | Annotate | Download | only in neteq
      1 /*
      2  *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
      3  *
      4  *  Use of this source code is governed by a BSD-style license
      5  *  that can be found in the LICENSE file in the root of the source
      6  *  tree. An additional intellectual property rights grant can be found
      7  *  in the file PATENTS.  All contributing project authors may
      8  *  be found in the AUTHORS file in the root of the source tree.
      9  */
     10 
     11 #include "webrtc/modules/audio_coding/neteq/dsp_helper.h"
     12 
     13 #include <assert.h>
     14 #include <string.h>  // Access to memset.
     15 
     16 #include <algorithm>  // Access to min, max.
     17 
     18 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
     19 
     20 namespace webrtc {
     21 
     22 // Table of constants used in method DspHelper::ParabolicFit().
     23 const int16_t DspHelper::kParabolaCoefficients[17][3] = {
     24     { 120, 32, 64 },
     25     { 140, 44, 75 },
     26     { 150, 50, 80 },
     27     { 160, 57, 85 },
     28     { 180, 72, 96 },
     29     { 200, 89, 107 },
     30     { 210, 98, 112 },
     31     { 220, 108, 117 },
     32     { 240, 128, 128 },
     33     { 260, 150, 139 },
     34     { 270, 162, 144 },
     35     { 280, 174, 149 },
     36     { 300, 200, 160 },
     37     { 320, 228, 171 },
     38     { 330, 242, 176 },
     39     { 340, 257, 181 },
     40     { 360, 288, 192 } };
     41 
     42 // Filter coefficients used when downsampling from the indicated sample rates
     43 // (8, 16, 32, 48 kHz) to 4 kHz. Coefficients are in Q12. The corresponding Q0
     44 // values are provided in the comments before each array.
     45 
     46 // Q0 values: {0.3, 0.4, 0.3}.
     47 const int16_t DspHelper::kDownsample8kHzTbl[3] = { 1229, 1638, 1229 };
     48 
     49 // Q0 values: {0.15, 0.2, 0.3, 0.2, 0.15}.
     50 const int16_t DspHelper::kDownsample16kHzTbl[5] = { 614, 819, 1229, 819, 614 };
     51 
     52 // Q0 values: {0.1425, 0.1251, 0.1525, 0.1628, 0.1525, 0.1251, 0.1425}.
     53 const int16_t DspHelper::kDownsample32kHzTbl[7] = {
     54     584, 512, 625, 667, 625, 512, 584 };
     55 
     56 // Q0 values: {0.2487, 0.0952, 0.1042, 0.1074, 0.1042, 0.0952, 0.2487}.
     57 const int16_t DspHelper::kDownsample48kHzTbl[7] = {
     58     1019, 390, 427, 440, 427, 390, 1019 };
     59 
     60 int DspHelper::RampSignal(const int16_t* input,
     61                           size_t length,
     62                           int factor,
     63                           int increment,
     64                           int16_t* output) {
     65   int factor_q20 = (factor << 6) + 32;
     66   // TODO(hlundin): Add 32 to factor_q20 when converting back to Q14?
     67   for (size_t i = 0; i < length; ++i) {
     68     output[i] = (factor * input[i] + 8192) >> 14;
     69     factor_q20 += increment;
     70     factor_q20 = std::max(factor_q20, 0);  // Never go negative.
     71     factor = std::min(factor_q20 >> 6, 16384);
     72   }
     73   return factor;
     74 }
     75 
     76 int DspHelper::RampSignal(int16_t* signal,
     77                           size_t length,
     78                           int factor,
     79                           int increment) {
     80   return RampSignal(signal, length, factor, increment, signal);
     81 }
     82 
     83 int DspHelper::RampSignal(AudioMultiVector* signal,
     84                           size_t start_index,
     85                           size_t length,
     86                           int factor,
     87                           int increment) {
     88   assert(start_index + length <= signal->Size());
     89   if (start_index + length > signal->Size()) {
     90     // Wrong parameters. Do nothing and return the scale factor unaltered.
     91     return factor;
     92   }
     93   int end_factor = 0;
     94   // Loop over the channels, starting at the same |factor| each time.
     95   for (size_t channel = 0; channel < signal->Channels(); ++channel) {
     96     end_factor =
     97         RampSignal(&(*signal)[channel][start_index], length, factor, increment);
     98   }
     99   return end_factor;
    100 }
    101 
    102 void DspHelper::PeakDetection(int16_t* data, int data_length,
    103                               int num_peaks, int fs_mult,
    104                               int* peak_index, int16_t* peak_value) {
    105   int16_t min_index = 0;
    106   int16_t max_index = 0;
    107 
    108   for (int i = 0; i <= num_peaks - 1; i++) {
    109     if (num_peaks == 1) {
    110       // Single peak.  The parabola fit assumes that an extra point is
    111       // available; worst case it gets a zero on the high end of the signal.
    112       // TODO(hlundin): This can potentially get much worse. It breaks the
    113       // API contract, that the length of |data| is |data_length|.
    114       data_length++;
    115     }
    116 
    117     peak_index[i] = WebRtcSpl_MaxIndexW16(data, data_length - 1);
    118 
    119     if (i != num_peaks - 1) {
    120       min_index = std::max(0, peak_index[i] - 2);
    121       max_index = std::min(data_length - 1, peak_index[i] + 2);
    122     }
    123 
    124     if ((peak_index[i] != 0) && (peak_index[i] != (data_length - 2))) {
    125       ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
    126                    &peak_value[i]);
    127     } else {
    128       if (peak_index[i] == data_length - 2) {
    129         if (data[peak_index[i]] > data[peak_index[i] + 1]) {
    130           ParabolicFit(&data[peak_index[i] - 1], fs_mult, &peak_index[i],
    131                        &peak_value[i]);
    132         } else if (data[peak_index[i]] <= data[peak_index[i] + 1]) {
    133           // Linear approximation.
    134           peak_value[i] = (data[peak_index[i]] + data[peak_index[i] + 1]) >> 1;
    135           peak_index[i] = (peak_index[i] * 2 + 1) * fs_mult;
    136         }
    137       } else {
    138         peak_value[i] = data[peak_index[i]];
    139         peak_index[i] = peak_index[i] * 2 * fs_mult;
    140       }
    141     }
    142 
    143     if (i != num_peaks - 1) {
    144       memset(&data[min_index], 0,
    145              sizeof(data[0]) * (max_index - min_index + 1));
    146     }
    147   }
    148 }
    149 
    150 void DspHelper::ParabolicFit(int16_t* signal_points, int fs_mult,
    151                              int* peak_index, int16_t* peak_value) {
    152   uint16_t fit_index[13];
    153   if (fs_mult == 1) {
    154     fit_index[0] = 0;
    155     fit_index[1] = 8;
    156     fit_index[2] = 16;
    157   } else if (fs_mult == 2) {
    158     fit_index[0] = 0;
    159     fit_index[1] = 4;
    160     fit_index[2] = 8;
    161     fit_index[3] = 12;
    162     fit_index[4] = 16;
    163   } else if (fs_mult == 4) {
    164     fit_index[0] = 0;
    165     fit_index[1] = 2;
    166     fit_index[2] = 4;
    167     fit_index[3] = 6;
    168     fit_index[4] = 8;
    169     fit_index[5] = 10;
    170     fit_index[6] = 12;
    171     fit_index[7] = 14;
    172     fit_index[8] = 16;
    173   } else {
    174     fit_index[0] = 0;
    175     fit_index[1] = 1;
    176     fit_index[2] = 3;
    177     fit_index[3] = 4;
    178     fit_index[4] = 5;
    179     fit_index[5] = 7;
    180     fit_index[6] = 8;
    181     fit_index[7] = 9;
    182     fit_index[8] = 11;
    183     fit_index[9] = 12;
    184     fit_index[10] = 13;
    185     fit_index[11] = 15;
    186     fit_index[12] = 16;
    187   }
    188 
    189   //  num = -3 * signal_points[0] + 4 * signal_points[1] - signal_points[2];
    190   //  den =      signal_points[0] - 2 * signal_points[1] + signal_points[2];
    191   int32_t num = (signal_points[0] * -3) + (signal_points[1] * 4)
    192       - signal_points[2];
    193   int32_t den = signal_points[0] + (signal_points[1] * -2) + signal_points[2];
    194   int32_t temp = num * 120;
    195   int flag = 1;
    196   int16_t stp = kParabolaCoefficients[fit_index[fs_mult]][0]
    197       - kParabolaCoefficients[fit_index[fs_mult - 1]][0];
    198   int16_t strt = (kParabolaCoefficients[fit_index[fs_mult]][0]
    199       + kParabolaCoefficients[fit_index[fs_mult - 1]][0]) / 2;
    200   int16_t lmt;
    201   if (temp < -den * strt) {
    202     lmt = strt - stp;
    203     while (flag) {
    204       if ((flag == fs_mult) || (temp > -den * lmt)) {
    205         *peak_value = (den * kParabolaCoefficients[fit_index[fs_mult - flag]][1]
    206             + num * kParabolaCoefficients[fit_index[fs_mult - flag]][2]
    207             + signal_points[0] * 256) / 256;
    208         *peak_index = *peak_index * 2 * fs_mult - flag;
    209         flag = 0;
    210       } else {
    211         flag++;
    212         lmt -= stp;
    213       }
    214     }
    215   } else if (temp > -den * (strt + stp)) {
    216     lmt = strt + 2 * stp;
    217     while (flag) {
    218       if ((flag == fs_mult) || (temp < -den * lmt)) {
    219         int32_t temp_term_1 =
    220             den * kParabolaCoefficients[fit_index[fs_mult+flag]][1];
    221         int32_t temp_term_2 =
    222             num * kParabolaCoefficients[fit_index[fs_mult+flag]][2];
    223         int32_t temp_term_3 = signal_points[0] * 256;
    224         *peak_value = (temp_term_1 + temp_term_2 + temp_term_3) / 256;
    225         *peak_index = *peak_index * 2 * fs_mult + flag;
    226         flag = 0;
    227       } else {
    228         flag++;
    229         lmt += stp;
    230       }
    231     }
    232   } else {
    233     *peak_value = signal_points[1];
    234     *peak_index = *peak_index * 2 * fs_mult;
    235   }
    236 }
    237 
    238 int DspHelper::MinDistortion(const int16_t* signal, int min_lag,
    239                              int max_lag, int length,
    240                              int32_t* distortion_value) {
    241   int best_index = -1;
    242   int32_t min_distortion = WEBRTC_SPL_WORD32_MAX;
    243   for (int i = min_lag; i <= max_lag; i++) {
    244     int32_t sum_diff = 0;
    245     const int16_t* data1 = signal;
    246     const int16_t* data2 = signal - i;
    247     for (int j = 0; j < length; j++) {
    248       sum_diff += WEBRTC_SPL_ABS_W32(data1[j] - data2[j]);
    249     }
    250     // Compare with previous minimum.
    251     if (sum_diff < min_distortion) {
    252       min_distortion = sum_diff;
    253       best_index = i;
    254     }
    255   }
    256   *distortion_value = min_distortion;
    257   return best_index;
    258 }
    259 
    260 void DspHelper::CrossFade(const int16_t* input1, const int16_t* input2,
    261                           size_t length, int16_t* mix_factor,
    262                           int16_t factor_decrement, int16_t* output) {
    263   int16_t factor = *mix_factor;
    264   int16_t complement_factor = 16384 - factor;
    265   for (size_t i = 0; i < length; i++) {
    266     output[i] =
    267         (factor * input1[i] + complement_factor * input2[i] + 8192) >> 14;
    268     factor -= factor_decrement;
    269     complement_factor += factor_decrement;
    270   }
    271   *mix_factor = factor;
    272 }
    273 
    274 void DspHelper::UnmuteSignal(const int16_t* input, size_t length,
    275                              int16_t* factor, int16_t increment,
    276                              int16_t* output) {
    277   uint16_t factor_16b = *factor;
    278   int32_t factor_32b = (static_cast<int32_t>(factor_16b) << 6) + 32;
    279   for (size_t i = 0; i < length; i++) {
    280     output[i] = (factor_16b * input[i] + 8192) >> 14;
    281     factor_32b = std::max(factor_32b + increment, 0);
    282     factor_16b = std::min(16384, factor_32b >> 6);
    283   }
    284   *factor = factor_16b;
    285 }
    286 
    287 void DspHelper::MuteSignal(int16_t* signal, int16_t mute_slope, size_t length) {
    288   int32_t factor = (16384 << 6) + 32;
    289   for (size_t i = 0; i < length; i++) {
    290     signal[i] = ((factor >> 6) * signal[i] + 8192) >> 14;
    291     factor -= mute_slope;
    292   }
    293 }
    294 
    295 int DspHelper::DownsampleTo4kHz(const int16_t* input, size_t input_length,
    296                                 int output_length, int input_rate_hz,
    297                                 bool compensate_delay, int16_t* output) {
    298   // Set filter parameters depending on input frequency.
    299   // NOTE: The phase delay values are wrong compared to the true phase delay
    300   // of the filters. However, the error is preserved (through the +1 term) for
    301   // consistency.
    302   const int16_t* filter_coefficients;  // Filter coefficients.
    303   int16_t filter_length;  // Number of coefficients.
    304   int16_t filter_delay;  // Phase delay in samples.
    305   int16_t factor;  // Conversion rate (inFsHz / 8000).
    306   switch (input_rate_hz) {
    307     case 8000: {
    308       filter_length = 3;
    309       factor = 2;
    310       filter_coefficients = kDownsample8kHzTbl;
    311       filter_delay = 1 + 1;
    312       break;
    313     }
    314     case 16000: {
    315       filter_length = 5;
    316       factor = 4;
    317       filter_coefficients = kDownsample16kHzTbl;
    318       filter_delay = 2 + 1;
    319       break;
    320     }
    321     case 32000: {
    322       filter_length = 7;
    323       factor = 8;
    324       filter_coefficients = kDownsample32kHzTbl;
    325       filter_delay = 3 + 1;
    326       break;
    327     }
    328     case 48000: {
    329       filter_length = 7;
    330       factor = 12;
    331       filter_coefficients = kDownsample48kHzTbl;
    332       filter_delay = 3 + 1;
    333       break;
    334     }
    335     default: {
    336       assert(false);
    337       return -1;
    338     }
    339   }
    340 
    341   if (!compensate_delay) {
    342     // Disregard delay compensation.
    343     filter_delay = 0;
    344   }
    345 
    346   // Returns -1 if input signal is too short; 0 otherwise.
    347   return WebRtcSpl_DownsampleFast(
    348       &input[filter_length - 1], static_cast<int>(input_length) -
    349       (filter_length - 1), output, output_length, filter_coefficients,
    350       filter_length, factor, filter_delay);
    351 }
    352 
    353 }  // namespace webrtc
    354