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/common_audio/vad/vad_filterbank.h" 12 13 #include <assert.h> 14 15 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h" 16 #include "webrtc/typedefs.h" 17 18 // Constants used in LogOfEnergy(). 19 static const int16_t kLogConst = 24660; // 160*log10(2) in Q9. 20 static const int16_t kLogEnergyIntPart = 14336; // 14 in Q10 21 22 // Coefficients used by HighPassFilter, Q14. 23 static const int16_t kHpZeroCoefs[3] = { 6631, -13262, 6631 }; 24 static const int16_t kHpPoleCoefs[3] = { 16384, -7756, 5620 }; 25 26 // Allpass filter coefficients, upper and lower, in Q15. 27 // Upper: 0.64, Lower: 0.17 28 static const int16_t kAllPassCoefsQ15[2] = { 20972, 5571 }; 29 30 // Adjustment for division with two in SplitFilter. 31 static const int16_t kOffsetVector[6] = { 368, 368, 272, 176, 176, 176 }; 32 33 // High pass filtering, with a cut-off frequency at 80 Hz, if the |data_in| is 34 // sampled at 500 Hz. 35 // 36 // - data_in [i] : Input audio data sampled at 500 Hz. 37 // - data_length [i] : Length of input and output data. 38 // - filter_state [i/o] : State of the filter. 39 // - data_out [o] : Output audio data in the frequency interval 40 // 80 - 250 Hz. 41 static void HighPassFilter(const int16_t* data_in, size_t data_length, 42 int16_t* filter_state, int16_t* data_out) { 43 size_t i; 44 const int16_t* in_ptr = data_in; 45 int16_t* out_ptr = data_out; 46 int32_t tmp32 = 0; 47 48 49 // The sum of the absolute values of the impulse response: 50 // The zero/pole-filter has a max amplification of a single sample of: 1.4546 51 // Impulse response: 0.4047 -0.6179 -0.0266 0.1993 0.1035 -0.0194 52 // The all-zero section has a max amplification of a single sample of: 1.6189 53 // Impulse response: 0.4047 -0.8094 0.4047 0 0 0 54 // The all-pole section has a max amplification of a single sample of: 1.9931 55 // Impulse response: 1.0000 0.4734 -0.1189 -0.2187 -0.0627 0.04532 56 57 for (i = 0; i < data_length; i++) { 58 // All-zero section (filter coefficients in Q14). 59 tmp32 = kHpZeroCoefs[0] * *in_ptr; 60 tmp32 += kHpZeroCoefs[1] * filter_state[0]; 61 tmp32 += kHpZeroCoefs[2] * filter_state[1]; 62 filter_state[1] = filter_state[0]; 63 filter_state[0] = *in_ptr++; 64 65 // All-pole section (filter coefficients in Q14). 66 tmp32 -= kHpPoleCoefs[1] * filter_state[2]; 67 tmp32 -= kHpPoleCoefs[2] * filter_state[3]; 68 filter_state[3] = filter_state[2]; 69 filter_state[2] = (int16_t) (tmp32 >> 14); 70 *out_ptr++ = filter_state[2]; 71 } 72 } 73 74 // All pass filtering of |data_in|, used before splitting the signal into two 75 // frequency bands (low pass vs high pass). 76 // Note that |data_in| and |data_out| can NOT correspond to the same address. 77 // 78 // - data_in [i] : Input audio signal given in Q0. 79 // - data_length [i] : Length of input and output data. 80 // - filter_coefficient [i] : Given in Q15. 81 // - filter_state [i/o] : State of the filter given in Q(-1). 82 // - data_out [o] : Output audio signal given in Q(-1). 83 static void AllPassFilter(const int16_t* data_in, size_t data_length, 84 int16_t filter_coefficient, int16_t* filter_state, 85 int16_t* data_out) { 86 // The filter can only cause overflow (in the w16 output variable) 87 // if more than 4 consecutive input numbers are of maximum value and 88 // has the the same sign as the impulse responses first taps. 89 // First 6 taps of the impulse response: 90 // 0.6399 0.5905 -0.3779 0.2418 -0.1547 0.0990 91 92 size_t i; 93 int16_t tmp16 = 0; 94 int32_t tmp32 = 0; 95 int32_t state32 = ((int32_t) (*filter_state) << 16); // Q15 96 97 for (i = 0; i < data_length; i++) { 98 tmp32 = state32 + filter_coefficient * *data_in; 99 tmp16 = (int16_t) (tmp32 >> 16); // Q(-1) 100 *data_out++ = tmp16; 101 state32 = (*data_in << 14) - filter_coefficient * tmp16; // Q14 102 state32 <<= 1; // Q15. 103 data_in += 2; 104 } 105 106 *filter_state = (int16_t) (state32 >> 16); // Q(-1) 107 } 108 109 // Splits |data_in| into |hp_data_out| and |lp_data_out| corresponding to 110 // an upper (high pass) part and a lower (low pass) part respectively. 111 // 112 // - data_in [i] : Input audio data to be split into two frequency bands. 113 // - data_length [i] : Length of |data_in|. 114 // - upper_state [i/o] : State of the upper filter, given in Q(-1). 115 // - lower_state [i/o] : State of the lower filter, given in Q(-1). 116 // - hp_data_out [o] : Output audio data of the upper half of the spectrum. 117 // The length is |data_length| / 2. 118 // - lp_data_out [o] : Output audio data of the lower half of the spectrum. 119 // The length is |data_length| / 2. 120 static void SplitFilter(const int16_t* data_in, size_t data_length, 121 int16_t* upper_state, int16_t* lower_state, 122 int16_t* hp_data_out, int16_t* lp_data_out) { 123 size_t i; 124 size_t half_length = data_length >> 1; // Downsampling by 2. 125 int16_t tmp_out; 126 127 // All-pass filtering upper branch. 128 AllPassFilter(&data_in[0], half_length, kAllPassCoefsQ15[0], upper_state, 129 hp_data_out); 130 131 // All-pass filtering lower branch. 132 AllPassFilter(&data_in[1], half_length, kAllPassCoefsQ15[1], lower_state, 133 lp_data_out); 134 135 // Make LP and HP signals. 136 for (i = 0; i < half_length; i++) { 137 tmp_out = *hp_data_out; 138 *hp_data_out++ -= *lp_data_out; 139 *lp_data_out++ += tmp_out; 140 } 141 } 142 143 // Calculates the energy of |data_in| in dB, and also updates an overall 144 // |total_energy| if necessary. 145 // 146 // - data_in [i] : Input audio data for energy calculation. 147 // - data_length [i] : Length of input data. 148 // - offset [i] : Offset value added to |log_energy|. 149 // - total_energy [i/o] : An external energy updated with the energy of 150 // |data_in|. 151 // NOTE: |total_energy| is only updated if 152 // |total_energy| <= |kMinEnergy|. 153 // - log_energy [o] : 10 * log10("energy of |data_in|") given in Q4. 154 static void LogOfEnergy(const int16_t* data_in, size_t data_length, 155 int16_t offset, int16_t* total_energy, 156 int16_t* log_energy) { 157 // |tot_rshifts| accumulates the number of right shifts performed on |energy|. 158 int tot_rshifts = 0; 159 // The |energy| will be normalized to 15 bits. We use unsigned integer because 160 // we eventually will mask out the fractional part. 161 uint32_t energy = 0; 162 163 assert(data_in != NULL); 164 assert(data_length > 0); 165 166 energy = (uint32_t) WebRtcSpl_Energy((int16_t*) data_in, data_length, 167 &tot_rshifts); 168 169 if (energy != 0) { 170 // By construction, normalizing to 15 bits is equivalent with 17 leading 171 // zeros of an unsigned 32 bit value. 172 int normalizing_rshifts = 17 - WebRtcSpl_NormU32(energy); 173 // In a 15 bit representation the leading bit is 2^14. log2(2^14) in Q10 is 174 // (14 << 10), which is what we initialize |log2_energy| with. For a more 175 // detailed derivations, see below. 176 int16_t log2_energy = kLogEnergyIntPart; 177 178 tot_rshifts += normalizing_rshifts; 179 // Normalize |energy| to 15 bits. 180 // |tot_rshifts| is now the total number of right shifts performed on 181 // |energy| after normalization. This means that |energy| is in 182 // Q(-tot_rshifts). 183 if (normalizing_rshifts < 0) { 184 energy <<= -normalizing_rshifts; 185 } else { 186 energy >>= normalizing_rshifts; 187 } 188 189 // Calculate the energy of |data_in| in dB, in Q4. 190 // 191 // 10 * log10("true energy") in Q4 = 2^4 * 10 * log10("true energy") = 192 // 160 * log10(|energy| * 2^|tot_rshifts|) = 193 // 160 * log10(2) * log2(|energy| * 2^|tot_rshifts|) = 194 // 160 * log10(2) * (log2(|energy|) + log2(2^|tot_rshifts|)) = 195 // (160 * log10(2)) * (log2(|energy|) + |tot_rshifts|) = 196 // |kLogConst| * (|log2_energy| + |tot_rshifts|) 197 // 198 // We know by construction that |energy| is normalized to 15 bits. Hence, 199 // |energy| = 2^14 + frac_Q15, where frac_Q15 is a fractional part in Q15. 200 // Further, we'd like |log2_energy| in Q10 201 // log2(|energy|) in Q10 = 2^10 * log2(2^14 + frac_Q15) = 202 // 2^10 * log2(2^14 * (1 + frac_Q15 * 2^-14)) = 203 // 2^10 * (14 + log2(1 + frac_Q15 * 2^-14)) ~= 204 // (14 << 10) + 2^10 * (frac_Q15 * 2^-14) = 205 // (14 << 10) + (frac_Q15 * 2^-4) = (14 << 10) + (frac_Q15 >> 4) 206 // 207 // Note that frac_Q15 = (|energy| & 0x00003FFF) 208 209 // Calculate and add the fractional part to |log2_energy|. 210 log2_energy += (int16_t) ((energy & 0x00003FFF) >> 4); 211 212 // |kLogConst| is in Q9, |log2_energy| in Q10 and |tot_rshifts| in Q0. 213 // Note that we in our derivation above have accounted for an output in Q4. 214 *log_energy = (int16_t)(((kLogConst * log2_energy) >> 19) + 215 ((tot_rshifts * kLogConst) >> 9)); 216 217 if (*log_energy < 0) { 218 *log_energy = 0; 219 } 220 } else { 221 *log_energy = offset; 222 return; 223 } 224 225 *log_energy += offset; 226 227 // Update the approximate |total_energy| with the energy of |data_in|, if 228 // |total_energy| has not exceeded |kMinEnergy|. |total_energy| is used as an 229 // energy indicator in WebRtcVad_GmmProbability() in vad_core.c. 230 if (*total_energy <= kMinEnergy) { 231 if (tot_rshifts >= 0) { 232 // We know by construction that the |energy| > |kMinEnergy| in Q0, so add 233 // an arbitrary value such that |total_energy| exceeds |kMinEnergy|. 234 *total_energy += kMinEnergy + 1; 235 } else { 236 // By construction |energy| is represented by 15 bits, hence any number of 237 // right shifted |energy| will fit in an int16_t. In addition, adding the 238 // value to |total_energy| is wrap around safe as long as 239 // |kMinEnergy| < 8192. 240 *total_energy += (int16_t) (energy >> -tot_rshifts); // Q0. 241 } 242 } 243 } 244 245 int16_t WebRtcVad_CalculateFeatures(VadInstT* self, const int16_t* data_in, 246 size_t data_length, int16_t* features) { 247 int16_t total_energy = 0; 248 // We expect |data_length| to be 80, 160 or 240 samples, which corresponds to 249 // 10, 20 or 30 ms in 8 kHz. Therefore, the intermediate downsampled data will 250 // have at most 120 samples after the first split and at most 60 samples after 251 // the second split. 252 int16_t hp_120[120], lp_120[120]; 253 int16_t hp_60[60], lp_60[60]; 254 const size_t half_data_length = data_length >> 1; 255 size_t length = half_data_length; // |data_length| / 2, corresponds to 256 // bandwidth = 2000 Hz after downsampling. 257 258 // Initialize variables for the first SplitFilter(). 259 int frequency_band = 0; 260 const int16_t* in_ptr = data_in; // [0 - 4000] Hz. 261 int16_t* hp_out_ptr = hp_120; // [2000 - 4000] Hz. 262 int16_t* lp_out_ptr = lp_120; // [0 - 2000] Hz. 263 264 assert(data_length <= 240); 265 assert(4 < kNumChannels - 1); // Checking maximum |frequency_band|. 266 267 // Split at 2000 Hz and downsample. 268 SplitFilter(in_ptr, data_length, &self->upper_state[frequency_band], 269 &self->lower_state[frequency_band], hp_out_ptr, lp_out_ptr); 270 271 // For the upper band (2000 Hz - 4000 Hz) split at 3000 Hz and downsample. 272 frequency_band = 1; 273 in_ptr = hp_120; // [2000 - 4000] Hz. 274 hp_out_ptr = hp_60; // [3000 - 4000] Hz. 275 lp_out_ptr = lp_60; // [2000 - 3000] Hz. 276 SplitFilter(in_ptr, length, &self->upper_state[frequency_band], 277 &self->lower_state[frequency_band], hp_out_ptr, lp_out_ptr); 278 279 // Energy in 3000 Hz - 4000 Hz. 280 length >>= 1; // |data_length| / 4 <=> bandwidth = 1000 Hz. 281 282 LogOfEnergy(hp_60, length, kOffsetVector[5], &total_energy, &features[5]); 283 284 // Energy in 2000 Hz - 3000 Hz. 285 LogOfEnergy(lp_60, length, kOffsetVector[4], &total_energy, &features[4]); 286 287 // For the lower band (0 Hz - 2000 Hz) split at 1000 Hz and downsample. 288 frequency_band = 2; 289 in_ptr = lp_120; // [0 - 2000] Hz. 290 hp_out_ptr = hp_60; // [1000 - 2000] Hz. 291 lp_out_ptr = lp_60; // [0 - 1000] Hz. 292 length = half_data_length; // |data_length| / 2 <=> bandwidth = 2000 Hz. 293 SplitFilter(in_ptr, length, &self->upper_state[frequency_band], 294 &self->lower_state[frequency_band], hp_out_ptr, lp_out_ptr); 295 296 // Energy in 1000 Hz - 2000 Hz. 297 length >>= 1; // |data_length| / 4 <=> bandwidth = 1000 Hz. 298 LogOfEnergy(hp_60, length, kOffsetVector[3], &total_energy, &features[3]); 299 300 // For the lower band (0 Hz - 1000 Hz) split at 500 Hz and downsample. 301 frequency_band = 3; 302 in_ptr = lp_60; // [0 - 1000] Hz. 303 hp_out_ptr = hp_120; // [500 - 1000] Hz. 304 lp_out_ptr = lp_120; // [0 - 500] Hz. 305 SplitFilter(in_ptr, length, &self->upper_state[frequency_band], 306 &self->lower_state[frequency_band], hp_out_ptr, lp_out_ptr); 307 308 // Energy in 500 Hz - 1000 Hz. 309 length >>= 1; // |data_length| / 8 <=> bandwidth = 500 Hz. 310 LogOfEnergy(hp_120, length, kOffsetVector[2], &total_energy, &features[2]); 311 312 // For the lower band (0 Hz - 500 Hz) split at 250 Hz and downsample. 313 frequency_band = 4; 314 in_ptr = lp_120; // [0 - 500] Hz. 315 hp_out_ptr = hp_60; // [250 - 500] Hz. 316 lp_out_ptr = lp_60; // [0 - 250] Hz. 317 SplitFilter(in_ptr, length, &self->upper_state[frequency_band], 318 &self->lower_state[frequency_band], hp_out_ptr, lp_out_ptr); 319 320 // Energy in 250 Hz - 500 Hz. 321 length >>= 1; // |data_length| / 16 <=> bandwidth = 250 Hz. 322 LogOfEnergy(hp_60, length, kOffsetVector[1], &total_energy, &features[1]); 323 324 // Remove 0 Hz - 80 Hz, by high pass filtering the lower band. 325 HighPassFilter(lp_60, length, self->hp_filter_state, hp_120); 326 327 // Energy in 80 Hz - 250 Hz. 328 LogOfEnergy(hp_120, length, kOffsetVector[0], &total_energy, &features[0]); 329 330 return total_energy; 331 } 332
