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