1/* Copyright (c) 2014, Cisco Systems, INC 2 Written by XiangMingZhu WeiZhou MinPeng YanWang 3 4 Redistribution and use in source and binary forms, with or without 5 modification, are permitted provided that the following conditions 6 are met: 7 8 - Redistributions of source code must retain the above copyright 9 notice, this list of conditions and the following disclaimer. 10 11 - Redistributions in binary form must reproduce the above copyright 12 notice, this list of conditions and the following disclaimer in the 13 documentation and/or other materials provided with the distribution. 14 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 16 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 17 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 18 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 19 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 20 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 21 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 22 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 23 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 24 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26*/ 27 28#ifdef HAVE_CONFIG_H 29#include "config.h" 30#endif 31 32#include <xmmintrin.h> 33#include <emmintrin.h> 34#include <smmintrin.h> 35 36#include "main.h" 37#include "stack_alloc.h" 38 39/* Weighting factors for tilt measure */ 40static const opus_int32 tiltWeights[ VAD_N_BANDS ] = { 30000, 6000, -12000, -12000 }; 41 42/***************************************/ 43/* Get the speech activity level in Q8 */ 44/***************************************/ 45opus_int silk_VAD_GetSA_Q8_sse4_1( /* O Return value, 0 if success */ 46 silk_encoder_state *psEncC, /* I/O Encoder state */ 47 const opus_int16 pIn[] /* I PCM input */ 48) 49{ 50 opus_int SA_Q15, pSNR_dB_Q7, input_tilt; 51 opus_int decimated_framelength1, decimated_framelength2; 52 opus_int decimated_framelength; 53 opus_int dec_subframe_length, dec_subframe_offset, SNR_Q7, i, b, s; 54 opus_int32 sumSquared, smooth_coef_Q16; 55 opus_int16 HPstateTmp; 56 VARDECL( opus_int16, X ); 57 opus_int32 Xnrg[ VAD_N_BANDS ]; 58 opus_int32 NrgToNoiseRatio_Q8[ VAD_N_BANDS ]; 59 opus_int32 speech_nrg, x_tmp; 60 opus_int X_offset[ VAD_N_BANDS ]; 61 opus_int ret = 0; 62 silk_VAD_state *psSilk_VAD = &psEncC->sVAD; 63 64 SAVE_STACK; 65 66 /* Safety checks */ 67 silk_assert( VAD_N_BANDS == 4 ); 68 silk_assert( MAX_FRAME_LENGTH >= psEncC->frame_length ); 69 silk_assert( psEncC->frame_length <= 512 ); 70 silk_assert( psEncC->frame_length == 8 * silk_RSHIFT( psEncC->frame_length, 3 ) ); 71 72 /***********************/ 73 /* Filter and Decimate */ 74 /***********************/ 75 decimated_framelength1 = silk_RSHIFT( psEncC->frame_length, 1 ); 76 decimated_framelength2 = silk_RSHIFT( psEncC->frame_length, 2 ); 77 decimated_framelength = silk_RSHIFT( psEncC->frame_length, 3 ); 78 /* Decimate into 4 bands: 79 0 L 3L L 3L 5L 80 - -- - -- -- 81 8 8 2 4 4 82 83 [0-1 kHz| temp. |1-2 kHz| 2-4 kHz | 4-8 kHz | 84 85 They're arranged to allow the minimal ( frame_length / 4 ) extra 86 scratch space during the downsampling process */ 87 X_offset[ 0 ] = 0; 88 X_offset[ 1 ] = decimated_framelength + decimated_framelength2; 89 X_offset[ 2 ] = X_offset[ 1 ] + decimated_framelength; 90 X_offset[ 3 ] = X_offset[ 2 ] + decimated_framelength2; 91 ALLOC( X, X_offset[ 3 ] + decimated_framelength1, opus_int16 ); 92 93 /* 0-8 kHz to 0-4 kHz and 4-8 kHz */ 94 silk_ana_filt_bank_1( pIn, &psSilk_VAD->AnaState[ 0 ], 95 X, &X[ X_offset[ 3 ] ], psEncC->frame_length ); 96 97 /* 0-4 kHz to 0-2 kHz and 2-4 kHz */ 98 silk_ana_filt_bank_1( X, &psSilk_VAD->AnaState1[ 0 ], 99 X, &X[ X_offset[ 2 ] ], decimated_framelength1 ); 100 101 /* 0-2 kHz to 0-1 kHz and 1-2 kHz */ 102 silk_ana_filt_bank_1( X, &psSilk_VAD->AnaState2[ 0 ], 103 X, &X[ X_offset[ 1 ] ], decimated_framelength2 ); 104 105 /*********************************************/ 106 /* HP filter on lowest band (differentiator) */ 107 /*********************************************/ 108 X[ decimated_framelength - 1 ] = silk_RSHIFT( X[ decimated_framelength - 1 ], 1 ); 109 HPstateTmp = X[ decimated_framelength - 1 ]; 110 for( i = decimated_framelength - 1; i > 0; i-- ) { 111 X[ i - 1 ] = silk_RSHIFT( X[ i - 1 ], 1 ); 112 X[ i ] -= X[ i - 1 ]; 113 } 114 X[ 0 ] -= psSilk_VAD->HPstate; 115 psSilk_VAD->HPstate = HPstateTmp; 116 117 /*************************************/ 118 /* Calculate the energy in each band */ 119 /*************************************/ 120 for( b = 0; b < VAD_N_BANDS; b++ ) { 121 /* Find the decimated framelength in the non-uniformly divided bands */ 122 decimated_framelength = silk_RSHIFT( psEncC->frame_length, silk_min_int( VAD_N_BANDS - b, VAD_N_BANDS - 1 ) ); 123 124 /* Split length into subframe lengths */ 125 dec_subframe_length = silk_RSHIFT( decimated_framelength, VAD_INTERNAL_SUBFRAMES_LOG2 ); 126 dec_subframe_offset = 0; 127 128 /* Compute energy per sub-frame */ 129 /* initialize with summed energy of last subframe */ 130 Xnrg[ b ] = psSilk_VAD->XnrgSubfr[ b ]; 131 for( s = 0; s < VAD_INTERNAL_SUBFRAMES; s++ ) { 132 __m128i xmm_X, xmm_acc; 133 sumSquared = 0; 134 135 xmm_acc = _mm_setzero_si128(); 136 137 for( i = 0; i < dec_subframe_length - 7; i += 8 ) 138 { 139 xmm_X = _mm_loadu_si128( (__m128i *)&(X[ X_offset[ b ] + i + dec_subframe_offset ] ) ); 140 xmm_X = _mm_srai_epi16( xmm_X, 3 ); 141 xmm_X = _mm_madd_epi16( xmm_X, xmm_X ); 142 xmm_acc = _mm_add_epi32( xmm_acc, xmm_X ); 143 } 144 145 xmm_acc = _mm_add_epi32( xmm_acc, _mm_unpackhi_epi64( xmm_acc, xmm_acc ) ); 146 xmm_acc = _mm_add_epi32( xmm_acc, _mm_shufflelo_epi16( xmm_acc, 0x0E ) ); 147 148 sumSquared += _mm_cvtsi128_si32( xmm_acc ); 149 150 for( ; i < dec_subframe_length; i++ ) { 151 /* The energy will be less than dec_subframe_length * ( silk_int16_MIN / 8 ) ^ 2. */ 152 /* Therefore we can accumulate with no risk of overflow (unless dec_subframe_length > 128) */ 153 x_tmp = silk_RSHIFT( 154 X[ X_offset[ b ] + i + dec_subframe_offset ], 3 ); 155 sumSquared = silk_SMLABB( sumSquared, x_tmp, x_tmp ); 156 157 /* Safety check */ 158 silk_assert( sumSquared >= 0 ); 159 } 160 161 /* Add/saturate summed energy of current subframe */ 162 if( s < VAD_INTERNAL_SUBFRAMES - 1 ) { 163 Xnrg[ b ] = silk_ADD_POS_SAT32( Xnrg[ b ], sumSquared ); 164 } else { 165 /* Look-ahead subframe */ 166 Xnrg[ b ] = silk_ADD_POS_SAT32( Xnrg[ b ], silk_RSHIFT( sumSquared, 1 ) ); 167 } 168 169 dec_subframe_offset += dec_subframe_length; 170 } 171 psSilk_VAD->XnrgSubfr[ b ] = sumSquared; 172 } 173 174 /********************/ 175 /* Noise estimation */ 176 /********************/ 177 silk_VAD_GetNoiseLevels( &Xnrg[ 0 ], psSilk_VAD ); 178 179 /***********************************************/ 180 /* Signal-plus-noise to noise ratio estimation */ 181 /***********************************************/ 182 sumSquared = 0; 183 input_tilt = 0; 184 for( b = 0; b < VAD_N_BANDS; b++ ) { 185 speech_nrg = Xnrg[ b ] - psSilk_VAD->NL[ b ]; 186 if( speech_nrg > 0 ) { 187 /* Divide, with sufficient resolution */ 188 if( ( Xnrg[ b ] & 0xFF800000 ) == 0 ) { 189 NrgToNoiseRatio_Q8[ b ] = silk_DIV32( silk_LSHIFT( Xnrg[ b ], 8 ), psSilk_VAD->NL[ b ] + 1 ); 190 } else { 191 NrgToNoiseRatio_Q8[ b ] = silk_DIV32( Xnrg[ b ], silk_RSHIFT( psSilk_VAD->NL[ b ], 8 ) + 1 ); 192 } 193 194 /* Convert to log domain */ 195 SNR_Q7 = silk_lin2log( NrgToNoiseRatio_Q8[ b ] ) - 8 * 128; 196 197 /* Sum-of-squares */ 198 sumSquared = silk_SMLABB( sumSquared, SNR_Q7, SNR_Q7 ); /* Q14 */ 199 200 /* Tilt measure */ 201 if( speech_nrg < ( (opus_int32)1 << 20 ) ) { 202 /* Scale down SNR value for small subband speech energies */ 203 SNR_Q7 = silk_SMULWB( silk_LSHIFT( silk_SQRT_APPROX( speech_nrg ), 6 ), SNR_Q7 ); 204 } 205 input_tilt = silk_SMLAWB( input_tilt, tiltWeights[ b ], SNR_Q7 ); 206 } else { 207 NrgToNoiseRatio_Q8[ b ] = 256; 208 } 209 } 210 211 /* Mean-of-squares */ 212 sumSquared = silk_DIV32_16( sumSquared, VAD_N_BANDS ); /* Q14 */ 213 214 /* Root-mean-square approximation, scale to dBs, and write to output pointer */ 215 pSNR_dB_Q7 = (opus_int16)( 3 * silk_SQRT_APPROX( sumSquared ) ); /* Q7 */ 216 217 /*********************************/ 218 /* Speech Probability Estimation */ 219 /*********************************/ 220 SA_Q15 = silk_sigm_Q15( silk_SMULWB( VAD_SNR_FACTOR_Q16, pSNR_dB_Q7 ) - VAD_NEGATIVE_OFFSET_Q5 ); 221 222 /**************************/ 223 /* Frequency Tilt Measure */ 224 /**************************/ 225 psEncC->input_tilt_Q15 = silk_LSHIFT( silk_sigm_Q15( input_tilt ) - 16384, 1 ); 226 227 /**************************************************/ 228 /* Scale the sigmoid output based on power levels */ 229 /**************************************************/ 230 speech_nrg = 0; 231 for( b = 0; b < VAD_N_BANDS; b++ ) { 232 /* Accumulate signal-without-noise energies, higher frequency bands have more weight */ 233 speech_nrg += ( b + 1 ) * silk_RSHIFT( Xnrg[ b ] - psSilk_VAD->NL[ b ], 4 ); 234 } 235 236 /* Power scaling */ 237 if( speech_nrg <= 0 ) { 238 SA_Q15 = silk_RSHIFT( SA_Q15, 1 ); 239 } else if( speech_nrg < 32768 ) { 240 if( psEncC->frame_length == 10 * psEncC->fs_kHz ) { 241 speech_nrg = silk_LSHIFT_SAT32( speech_nrg, 16 ); 242 } else { 243 speech_nrg = silk_LSHIFT_SAT32( speech_nrg, 15 ); 244 } 245 246 /* square-root */ 247 speech_nrg = silk_SQRT_APPROX( speech_nrg ); 248 SA_Q15 = silk_SMULWB( 32768 + speech_nrg, SA_Q15 ); 249 } 250 251 /* Copy the resulting speech activity in Q8 */ 252 psEncC->speech_activity_Q8 = silk_min_int( silk_RSHIFT( SA_Q15, 7 ), silk_uint8_MAX ); 253 254 /***********************************/ 255 /* Energy Level and SNR estimation */ 256 /***********************************/ 257 /* Smoothing coefficient */ 258 smooth_coef_Q16 = silk_SMULWB( VAD_SNR_SMOOTH_COEF_Q18, silk_SMULWB( (opus_int32)SA_Q15, SA_Q15 ) ); 259 260 if( psEncC->frame_length == 10 * psEncC->fs_kHz ) { 261 smooth_coef_Q16 >>= 1; 262 } 263 264 for( b = 0; b < VAD_N_BANDS; b++ ) { 265 /* compute smoothed energy-to-noise ratio per band */ 266 psSilk_VAD->NrgRatioSmth_Q8[ b ] = silk_SMLAWB( psSilk_VAD->NrgRatioSmth_Q8[ b ], 267 NrgToNoiseRatio_Q8[ b ] - psSilk_VAD->NrgRatioSmth_Q8[ b ], smooth_coef_Q16 ); 268 269 /* signal to noise ratio in dB per band */ 270 SNR_Q7 = 3 * ( silk_lin2log( psSilk_VAD->NrgRatioSmth_Q8[b] ) - 8 * 128 ); 271 /* quality = sigmoid( 0.25 * ( SNR_dB - 16 ) ); */ 272 psEncC->input_quality_bands_Q15[ b ] = silk_sigm_Q15( silk_RSHIFT( SNR_Q7 - 16 * 128, 4 ) ); 273 } 274 275 RESTORE_STACK; 276 return( ret ); 277} 278