Line data    Source code

       1             : /*
       2             :  *  Copyright (c) 2011 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/common_audio/vad/vad_gmm.h"
      12             : 
      13             : #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
      14             : #include "webrtc/typedefs.h"
      15             : 
      16             : static const int32_t kCompVar = 22005;
      17             : static const int16_t kLog2Exp = 5909;  // log2(exp(1)) in Q12.
      18             : 
      19             : // For a normal distribution, the probability of |input| is calculated and
      20             : // returned (in Q20). The formula for normal distributed probability is
      21             : //
      22             : // 1 / s * exp(-(x - m)^2 / (2 * s^2))
      23             : //
      24             : // where the parameters are given in the following Q domains:
      25             : // m = |mean| (Q7)
      26             : // s = |std| (Q7)
      27             : // x = |input| (Q4)
      28             : // in addition to the probability we output |delta| (in Q11) used when updating
      29             : // the noise/speech model.
      30           0 : int32_t WebRtcVad_GaussianProbability(int16_t input,
      31             :                                       int16_t mean,
      32             :                                       int16_t std,
      33             :                                       int16_t* delta) {
      34           0 :   int16_t tmp16, inv_std, inv_std2, exp_value = 0;
      35             :   int32_t tmp32;
      36             : 
      37             :   // Calculate |inv_std| = 1 / s, in Q10.
      38             :   // 131072 = 1 in Q17, and (|std| >> 1) is for rounding instead of truncation.
      39             :   // Q-domain: Q17 / Q7 = Q10.
      40           0 :   tmp32 = (int32_t) 131072 + (int32_t) (std >> 1);
      41           0 :   inv_std = (int16_t) WebRtcSpl_DivW32W16(tmp32, std);
      42             : 
      43             :   // Calculate |inv_std2| = 1 / s^2, in Q14.
      44           0 :   tmp16 = (inv_std >> 2);  // Q10 -> Q8.
      45             :   // Q-domain: (Q8 * Q8) >> 2 = Q14.
      46           0 :   inv_std2 = (int16_t)((tmp16 * tmp16) >> 2);
      47             :   // TODO(bjornv): Investigate if changing to
      48             :   // inv_std2 = (int16_t)((inv_std * inv_std) >> 6);
      49             :   // gives better accuracy.
      50             : 
      51           0 :   tmp16 = (input << 3);  // Q4 -> Q7
      52           0 :   tmp16 = tmp16 - mean;  // Q7 - Q7 = Q7
      53             : 
      54             :   // To be used later, when updating noise/speech model.
      55             :   // |delta| = (x - m) / s^2, in Q11.
      56             :   // Q-domain: (Q14 * Q7) >> 10 = Q11.
      57           0 :   *delta = (int16_t)((inv_std2 * tmp16) >> 10);
      58             : 
      59             :   // Calculate the exponent |tmp32| = (x - m)^2 / (2 * s^2), in Q10. Replacing
      60             :   // division by two with one shift.
      61             :   // Q-domain: (Q11 * Q7) >> 8 = Q10.
      62           0 :   tmp32 = (*delta * tmp16) >> 9;
      63             : 
      64             :   // If the exponent is small enough to give a non-zero probability we calculate
      65             :   // |exp_value| ~= exp(-(x - m)^2 / (2 * s^2))
      66             :   //             ~= exp2(-log2(exp(1)) * |tmp32|).
      67           0 :   if (tmp32 < kCompVar) {
      68             :     // Calculate |tmp16| = log2(exp(1)) * |tmp32|, in Q10.
      69             :     // Q-domain: (Q12 * Q10) >> 12 = Q10.
      70           0 :     tmp16 = (int16_t)((kLog2Exp * tmp32) >> 12);
      71           0 :     tmp16 = -tmp16;
      72           0 :     exp_value = (0x0400 | (tmp16 & 0x03FF));
      73           0 :     tmp16 ^= 0xFFFF;
      74           0 :     tmp16 >>= 10;
      75           0 :     tmp16 += 1;
      76             :     // Get |exp_value| = exp(-|tmp32|) in Q10.
      77           0 :     exp_value >>= tmp16;
      78             :   }
      79             : 
      80             :   // Calculate and return (1 / s) * exp(-(x - m)^2 / (2 * s^2)), in Q20.
      81             :   // Q-domain: Q10 * Q10 = Q20.
      82           0 :   return inv_std * exp_value;
      83             : }