michael@0: /***********************************************************************
michael@0: Copyright (c) 2006-2011, Skype Limited. All rights reserved.
michael@0: Redistribution and use in source and binary forms, with or without
michael@0: modification, are permitted provided that the following conditions
michael@0: are met:
michael@0: - Redistributions of source code must retain the above copyright notice,
michael@0: this list of conditions and the following disclaimer.
michael@0: - Redistributions in binary form must reproduce the above copyright
michael@0: notice, this list of conditions and the following disclaimer in the
michael@0: documentation and/or other materials provided with the distribution.
michael@0: - Neither the name of Internet Society, IETF or IETF Trust, nor the
michael@0: names of specific contributors, may be used to endorse or promote
michael@0: products derived from this software without specific prior written
michael@0: permission.
michael@0: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
michael@0: AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
michael@0: IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
michael@0: ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
michael@0: LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
michael@0: CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
michael@0: SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
michael@0: INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
michael@0: CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
michael@0: ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
michael@0: POSSIBILITY OF SUCH DAMAGE.
michael@0: ***********************************************************************/
michael@0: 
michael@0: #ifdef HAVE_CONFIG_H
michael@0: #include "config.h"
michael@0: #endif
michael@0: 
michael@0: #include "main.h"
michael@0: #include "stack_alloc.h"
michael@0: 
michael@0: /* Silk VAD noise level estimation */
michael@0: static OPUS_INLINE void silk_VAD_GetNoiseLevels(
michael@0:     const opus_int32             pX[ VAD_N_BANDS ], /* I    subband energies                            */
michael@0:     silk_VAD_state              *psSilk_VAD         /* I/O  Pointer to Silk VAD state                   */
michael@0: );
michael@0: 
michael@0: /**********************************/
michael@0: /* Initialization of the Silk VAD */
michael@0: /**********************************/
michael@0: opus_int silk_VAD_Init(                                         /* O    Return value, 0 if success                  */
michael@0:     silk_VAD_state              *psSilk_VAD                     /* I/O  Pointer to Silk VAD state                   */
michael@0: )
michael@0: {
michael@0:     opus_int b, ret = 0;
michael@0: 
michael@0:     /* reset state memory */
michael@0:     silk_memset( psSilk_VAD, 0, sizeof( silk_VAD_state ) );
michael@0: 
michael@0:     /* init noise levels */
michael@0:     /* Initialize array with approx pink noise levels (psd proportional to inverse of frequency) */
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         psSilk_VAD->NoiseLevelBias[ b ] = silk_max_32( silk_DIV32_16( VAD_NOISE_LEVELS_BIAS, b + 1 ), 1 );
michael@0:     }
michael@0: 
michael@0:     /* Initialize state */
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         psSilk_VAD->NL[ b ]     = silk_MUL( 100, psSilk_VAD->NoiseLevelBias[ b ] );
michael@0:         psSilk_VAD->inv_NL[ b ] = silk_DIV32( silk_int32_MAX, psSilk_VAD->NL[ b ] );
michael@0:     }
michael@0:     psSilk_VAD->counter = 15;
michael@0: 
michael@0:     /* init smoothed energy-to-noise ratio*/
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         psSilk_VAD->NrgRatioSmth_Q8[ b ] = 100 * 256;       /* 100 * 256 --> 20 dB SNR */
michael@0:     }
michael@0: 
michael@0:     return( ret );
michael@0: }
michael@0: 
michael@0: /* Weighting factors for tilt measure */
michael@0: static const opus_int32 tiltWeights[ VAD_N_BANDS ] = { 30000, 6000, -12000, -12000 };
michael@0: 
michael@0: /***************************************/
michael@0: /* Get the speech activity level in Q8 */
michael@0: /***************************************/
michael@0: opus_int silk_VAD_GetSA_Q8(                                     /* O    Return value, 0 if success                  */
michael@0:     silk_encoder_state          *psEncC,                        /* I/O  Encoder state                               */
michael@0:     const opus_int16            pIn[]                           /* I    PCM input                                   */
michael@0: )
michael@0: {
michael@0:     opus_int   SA_Q15, pSNR_dB_Q7, input_tilt;
michael@0:     opus_int   decimated_framelength1, decimated_framelength2;
michael@0:     opus_int   decimated_framelength;
michael@0:     opus_int   dec_subframe_length, dec_subframe_offset, SNR_Q7, i, b, s;
michael@0:     opus_int32 sumSquared, smooth_coef_Q16;
michael@0:     opus_int16 HPstateTmp;
michael@0:     VARDECL( opus_int16, X );
michael@0:     opus_int32 Xnrg[ VAD_N_BANDS ];
michael@0:     opus_int32 NrgToNoiseRatio_Q8[ VAD_N_BANDS ];
michael@0:     opus_int32 speech_nrg, x_tmp;
michael@0:     opus_int   X_offset[ VAD_N_BANDS ];
michael@0:     opus_int   ret = 0;
michael@0:     silk_VAD_state *psSilk_VAD = &psEncC->sVAD;
michael@0:     SAVE_STACK;
michael@0: 
michael@0:     /* Safety checks */
michael@0:     silk_assert( VAD_N_BANDS == 4 );
michael@0:     silk_assert( MAX_FRAME_LENGTH >= psEncC->frame_length );
michael@0:     silk_assert( psEncC->frame_length <= 512 );
michael@0:     silk_assert( psEncC->frame_length == 8 * silk_RSHIFT( psEncC->frame_length, 3 ) );
michael@0: 
michael@0:     /***********************/
michael@0:     /* Filter and Decimate */
michael@0:     /***********************/
michael@0:     decimated_framelength1 = silk_RSHIFT( psEncC->frame_length, 1 );
michael@0:     decimated_framelength2 = silk_RSHIFT( psEncC->frame_length, 2 );
michael@0:     decimated_framelength = silk_RSHIFT( psEncC->frame_length, 3 );
michael@0:     /* Decimate into 4 bands:
michael@0:        0       L      3L       L              3L                             5L
michael@0:                -      --       -              --                             --
michael@0:                8       8       2               4                              4
michael@0: 
michael@0:        [0-1 kHz| temp. |1-2 kHz|    2-4 kHz    |            4-8 kHz           |
michael@0: 
michael@0:        They're arranged to allow the minimal ( frame_length / 4 ) extra
michael@0:        scratch space during the downsampling process */
michael@0:     X_offset[ 0 ] = 0;
michael@0:     X_offset[ 1 ] = decimated_framelength + decimated_framelength2;
michael@0:     X_offset[ 2 ] = X_offset[ 1 ] + decimated_framelength;
michael@0:     X_offset[ 3 ] = X_offset[ 2 ] + decimated_framelength2;
michael@0:     ALLOC( X, X_offset[ 3 ] + decimated_framelength1, opus_int16 );
michael@0: 
michael@0:     /* 0-8 kHz to 0-4 kHz and 4-8 kHz */
michael@0:     silk_ana_filt_bank_1( pIn, &psSilk_VAD->AnaState[  0 ],
michael@0:         X, &X[ X_offset[ 3 ] ], psEncC->frame_length );
michael@0: 
michael@0:     /* 0-4 kHz to 0-2 kHz and 2-4 kHz */
michael@0:     silk_ana_filt_bank_1( X, &psSilk_VAD->AnaState1[ 0 ],
michael@0:         X, &X[ X_offset[ 2 ] ], decimated_framelength1 );
michael@0: 
michael@0:     /* 0-2 kHz to 0-1 kHz and 1-2 kHz */
michael@0:     silk_ana_filt_bank_1( X, &psSilk_VAD->AnaState2[ 0 ],
michael@0:         X, &X[ X_offset[ 1 ] ], decimated_framelength2 );
michael@0: 
michael@0:     /*********************************************/
michael@0:     /* HP filter on lowest band (differentiator) */
michael@0:     /*********************************************/
michael@0:     X[ decimated_framelength - 1 ] = silk_RSHIFT( X[ decimated_framelength - 1 ], 1 );
michael@0:     HPstateTmp = X[ decimated_framelength - 1 ];
michael@0:     for( i = decimated_framelength - 1; i > 0; i-- ) {
michael@0:         X[ i - 1 ]  = silk_RSHIFT( X[ i - 1 ], 1 );
michael@0:         X[ i ]     -= X[ i - 1 ];
michael@0:     }
michael@0:     X[ 0 ] -= psSilk_VAD->HPstate;
michael@0:     psSilk_VAD->HPstate = HPstateTmp;
michael@0: 
michael@0:     /*************************************/
michael@0:     /* Calculate the energy in each band */
michael@0:     /*************************************/
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         /* Find the decimated framelength in the non-uniformly divided bands */
michael@0:         decimated_framelength = silk_RSHIFT( psEncC->frame_length, silk_min_int( VAD_N_BANDS - b, VAD_N_BANDS - 1 ) );
michael@0: 
michael@0:         /* Split length into subframe lengths */
michael@0:         dec_subframe_length = silk_RSHIFT( decimated_framelength, VAD_INTERNAL_SUBFRAMES_LOG2 );
michael@0:         dec_subframe_offset = 0;
michael@0: 
michael@0:         /* Compute energy per sub-frame */
michael@0:         /* initialize with summed energy of last subframe */
michael@0:         Xnrg[ b ] = psSilk_VAD->XnrgSubfr[ b ];
michael@0:         for( s = 0; s < VAD_INTERNAL_SUBFRAMES; s++ ) {
michael@0:             sumSquared = 0;
michael@0:             for( i = 0; i < dec_subframe_length; i++ ) {
michael@0:                 /* The energy will be less than dec_subframe_length * ( silk_int16_MIN / 8 ) ^ 2.            */
michael@0:                 /* Therefore we can accumulate with no risk of overflow (unless dec_subframe_length > 128)  */
michael@0:                 x_tmp = silk_RSHIFT(
michael@0:                     X[ X_offset[ b ] + i + dec_subframe_offset ], 3 );
michael@0:                 sumSquared = silk_SMLABB( sumSquared, x_tmp, x_tmp );
michael@0: 
michael@0:                 /* Safety check */
michael@0:                 silk_assert( sumSquared >= 0 );
michael@0:             }
michael@0: 
michael@0:             /* Add/saturate summed energy of current subframe */
michael@0:             if( s < VAD_INTERNAL_SUBFRAMES - 1 ) {
michael@0:                 Xnrg[ b ] = silk_ADD_POS_SAT32( Xnrg[ b ], sumSquared );
michael@0:             } else {
michael@0:                 /* Look-ahead subframe */
michael@0:                 Xnrg[ b ] = silk_ADD_POS_SAT32( Xnrg[ b ], silk_RSHIFT( sumSquared, 1 ) );
michael@0:             }
michael@0: 
michael@0:             dec_subframe_offset += dec_subframe_length;
michael@0:         }
michael@0:         psSilk_VAD->XnrgSubfr[ b ] = sumSquared;
michael@0:     }
michael@0: 
michael@0:     /********************/
michael@0:     /* Noise estimation */
michael@0:     /********************/
michael@0:     silk_VAD_GetNoiseLevels( &Xnrg[ 0 ], psSilk_VAD );
michael@0: 
michael@0:     /***********************************************/
michael@0:     /* Signal-plus-noise to noise ratio estimation */
michael@0:     /***********************************************/
michael@0:     sumSquared = 0;
michael@0:     input_tilt = 0;
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         speech_nrg = Xnrg[ b ] - psSilk_VAD->NL[ b ];
michael@0:         if( speech_nrg > 0 ) {
michael@0:             /* Divide, with sufficient resolution */
michael@0:             if( ( Xnrg[ b ] & 0xFF800000 ) == 0 ) {
michael@0:                 NrgToNoiseRatio_Q8[ b ] = silk_DIV32( silk_LSHIFT( Xnrg[ b ], 8 ), psSilk_VAD->NL[ b ] + 1 );
michael@0:             } else {
michael@0:                 NrgToNoiseRatio_Q8[ b ] = silk_DIV32( Xnrg[ b ], silk_RSHIFT( psSilk_VAD->NL[ b ], 8 ) + 1 );
michael@0:             }
michael@0: 
michael@0:             /* Convert to log domain */
michael@0:             SNR_Q7 = silk_lin2log( NrgToNoiseRatio_Q8[ b ] ) - 8 * 128;
michael@0: 
michael@0:             /* Sum-of-squares */
michael@0:             sumSquared = silk_SMLABB( sumSquared, SNR_Q7, SNR_Q7 );          /* Q14 */
michael@0: 
michael@0:             /* Tilt measure */
michael@0:             if( speech_nrg < ( (opus_int32)1 << 20 ) ) {
michael@0:                 /* Scale down SNR value for small subband speech energies */
michael@0:                 SNR_Q7 = silk_SMULWB( silk_LSHIFT( silk_SQRT_APPROX( speech_nrg ), 6 ), SNR_Q7 );
michael@0:             }
michael@0:             input_tilt = silk_SMLAWB( input_tilt, tiltWeights[ b ], SNR_Q7 );
michael@0:         } else {
michael@0:             NrgToNoiseRatio_Q8[ b ] = 256;
michael@0:         }
michael@0:     }
michael@0: 
michael@0:     /* Mean-of-squares */
michael@0:     sumSquared = silk_DIV32_16( sumSquared, VAD_N_BANDS ); /* Q14 */
michael@0: 
michael@0:     /* Root-mean-square approximation, scale to dBs, and write to output pointer */
michael@0:     pSNR_dB_Q7 = (opus_int16)( 3 * silk_SQRT_APPROX( sumSquared ) ); /* Q7 */
michael@0: 
michael@0:     /*********************************/
michael@0:     /* Speech Probability Estimation */
michael@0:     /*********************************/
michael@0:     SA_Q15 = silk_sigm_Q15( silk_SMULWB( VAD_SNR_FACTOR_Q16, pSNR_dB_Q7 ) - VAD_NEGATIVE_OFFSET_Q5 );
michael@0: 
michael@0:     /**************************/
michael@0:     /* Frequency Tilt Measure */
michael@0:     /**************************/
michael@0:     psEncC->input_tilt_Q15 = silk_LSHIFT( silk_sigm_Q15( input_tilt ) - 16384, 1 );
michael@0: 
michael@0:     /**************************************************/
michael@0:     /* Scale the sigmoid output based on power levels */
michael@0:     /**************************************************/
michael@0:     speech_nrg = 0;
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         /* Accumulate signal-without-noise energies, higher frequency bands have more weight */
michael@0:         speech_nrg += ( b + 1 ) * silk_RSHIFT( Xnrg[ b ] - psSilk_VAD->NL[ b ], 4 );
michael@0:     }
michael@0: 
michael@0:     /* Power scaling */
michael@0:     if( speech_nrg <= 0 ) {
michael@0:         SA_Q15 = silk_RSHIFT( SA_Q15, 1 );
michael@0:     } else if( speech_nrg < 32768 ) {
michael@0:         if( psEncC->frame_length == 10 * psEncC->fs_kHz ) {
michael@0:             speech_nrg = silk_LSHIFT_SAT32( speech_nrg, 16 );
michael@0:         } else {
michael@0:             speech_nrg = silk_LSHIFT_SAT32( speech_nrg, 15 );
michael@0:         }
michael@0: 
michael@0:         /* square-root */
michael@0:         speech_nrg = silk_SQRT_APPROX( speech_nrg );
michael@0:         SA_Q15 = silk_SMULWB( 32768 + speech_nrg, SA_Q15 );
michael@0:     }
michael@0: 
michael@0:     /* Copy the resulting speech activity in Q8 */
michael@0:     psEncC->speech_activity_Q8 = silk_min_int( silk_RSHIFT( SA_Q15, 7 ), silk_uint8_MAX );
michael@0: 
michael@0:     /***********************************/
michael@0:     /* Energy Level and SNR estimation */
michael@0:     /***********************************/
michael@0:     /* Smoothing coefficient */
michael@0:     smooth_coef_Q16 = silk_SMULWB( VAD_SNR_SMOOTH_COEF_Q18, silk_SMULWB( (opus_int32)SA_Q15, SA_Q15 ) );
michael@0: 
michael@0:     if( psEncC->frame_length == 10 * psEncC->fs_kHz ) {
michael@0:         smooth_coef_Q16 >>= 1;
michael@0:     }
michael@0: 
michael@0:     for( b = 0; b < VAD_N_BANDS; b++ ) {
michael@0:         /* compute smoothed energy-to-noise ratio per band */
michael@0:         psSilk_VAD->NrgRatioSmth_Q8[ b ] = silk_SMLAWB( psSilk_VAD->NrgRatioSmth_Q8[ b ],
michael@0:             NrgToNoiseRatio_Q8[ b ] - psSilk_VAD->NrgRatioSmth_Q8[ b ], smooth_coef_Q16 );
michael@0: 
michael@0:         /* signal to noise ratio in dB per band */
michael@0:         SNR_Q7 = 3 * ( silk_lin2log( psSilk_VAD->NrgRatioSmth_Q8[b] ) - 8 * 128 );
michael@0:         /* quality = sigmoid( 0.25 * ( SNR_dB - 16 ) ); */
michael@0:         psEncC->input_quality_bands_Q15[ b ] = silk_sigm_Q15( silk_RSHIFT( SNR_Q7 - 16 * 128, 4 ) );
michael@0:     }
michael@0: 
michael@0:     RESTORE_STACK;
michael@0:     return( ret );
michael@0: }
michael@0: 
michael@0: /**************************/
michael@0: /* Noise level estimation */
michael@0: /**************************/
michael@0: static OPUS_INLINE void silk_VAD_GetNoiseLevels(
michael@0:     const opus_int32            pX[ VAD_N_BANDS ],  /* I    subband energies                            */
michael@0:     silk_VAD_state              *psSilk_VAD         /* I/O  Pointer to Silk VAD state                   */
michael@0: )
michael@0: {
michael@0:     opus_int   k;
michael@0:     opus_int32 nl, nrg, inv_nrg;
michael@0:     opus_int   coef, min_coef;
michael@0: 
michael@0:     /* Initially faster smoothing */
michael@0:     if( psSilk_VAD->counter < 1000 ) { /* 1000 = 20 sec */
michael@0:         min_coef = silk_DIV32_16( silk_int16_MAX, silk_RSHIFT( psSilk_VAD->counter, 4 ) + 1 );
michael@0:     } else {
michael@0:         min_coef = 0;
michael@0:     }
michael@0: 
michael@0:     for( k = 0; k < VAD_N_BANDS; k++ ) {
michael@0:         /* Get old noise level estimate for current band */
michael@0:         nl = psSilk_VAD->NL[ k ];
michael@0:         silk_assert( nl >= 0 );
michael@0: 
michael@0:         /* Add bias */
michael@0:         nrg = silk_ADD_POS_SAT32( pX[ k ], psSilk_VAD->NoiseLevelBias[ k ] );
michael@0:         silk_assert( nrg > 0 );
michael@0: 
michael@0:         /* Invert energies */
michael@0:         inv_nrg = silk_DIV32( silk_int32_MAX, nrg );
michael@0:         silk_assert( inv_nrg >= 0 );
michael@0: 
michael@0:         /* Less update when subband energy is high */
michael@0:         if( nrg > silk_LSHIFT( nl, 3 ) ) {
michael@0:             coef = VAD_NOISE_LEVEL_SMOOTH_COEF_Q16 >> 3;
michael@0:         } else if( nrg < nl ) {
michael@0:             coef = VAD_NOISE_LEVEL_SMOOTH_COEF_Q16;
michael@0:         } else {
michael@0:             coef = silk_SMULWB( silk_SMULWW( inv_nrg, nl ), VAD_NOISE_LEVEL_SMOOTH_COEF_Q16 << 1 );
michael@0:         }
michael@0: 
michael@0:         /* Initially faster smoothing */
michael@0:         coef = silk_max_int( coef, min_coef );
michael@0: 
michael@0:         /* Smooth inverse energies */
michael@0:         psSilk_VAD->inv_NL[ k ] = silk_SMLAWB( psSilk_VAD->inv_NL[ k ], inv_nrg - psSilk_VAD->inv_NL[ k ], coef );
michael@0:         silk_assert( psSilk_VAD->inv_NL[ k ] >= 0 );
michael@0: 
michael@0:         /* Compute noise level by inverting again */
michael@0:         nl = silk_DIV32( silk_int32_MAX, psSilk_VAD->inv_NL[ k ] );
michael@0:         silk_assert( nl >= 0 );
michael@0: 
michael@0:         /* Limit noise levels (guarantee 7 bits of head room) */
michael@0:         nl = silk_min( nl, 0x00FFFFFF );
michael@0: 
michael@0:         /* Store as part of state */
michael@0:         psSilk_VAD->NL[ k ] = nl;
michael@0:     }
michael@0: 
michael@0:     /* Increment frame counter */
michael@0:     psSilk_VAD->counter++;
michael@0: }