Initial Skype commit taken from FreeSwitch, which got it from the IETF draft.
[opus.git] / src / SKP_Silk_VAD.c
1 /***********************************************************************\r
2 Copyright (c) 2006-2010, Skype Limited. All rights reserved. \r
3 Redistribution and use in source and binary forms, with or without \r
4 modification, (subject to the limitations in the disclaimer below) \r
5 are permitted provided that the following conditions are met:\r
6 - Redistributions of source code must retain the above copyright notice,\r
7 this list of conditions and the following disclaimer.\r
8 - Redistributions in binary form must reproduce the above copyright \r
9 notice, this list of conditions and the following disclaimer in the \r
10 documentation and/or other materials provided with the distribution.\r
11 - Neither the name of Skype Limited, nor the names of specific \r
12 contributors, may be used to endorse or promote products derived from \r
13 this software without specific prior written permission.\r
14 NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED \r
15 BY THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND \r
16 CONTRIBUTORS ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,\r
17 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND \r
18 FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE \r
19 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, \r
20 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT\r
21 NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF \r
22 USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON \r
23 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT \r
24 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE \r
25 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
26 ***********************************************************************/\r
27 \r
28 /*\r
29  * File Name:   SKP_Silk_VAD.c\r
30  * Description: Silk VAD.\r
31  */\r
32 \r
33 #include <stdlib.h>\r
34 #include "SKP_Silk_main.h"\r
35 \r
36 /**********************************/\r
37 /* Initialization of the Silk VAD */\r
38 /**********************************/\r
39 SKP_int SKP_Silk_VAD_Init(                              /* O    Return value, 0 if success                  */ \r
40     SKP_Silk_VAD_state              *psSilk_VAD         /* I/O  Pointer to Silk VAD state                   */ \r
41 )\r
42 {\r
43     SKP_int b, ret = 0;\r
44 \r
45     /* reset state memory */\r
46     SKP_memset( psSilk_VAD, 0, sizeof( SKP_Silk_VAD_state ) );\r
47 \r
48     /* init noise levels */\r
49     /* Initialize array with approx pink noise levels (psd proportional to inverse of frequency) */\r
50     for( b = 0; b < VAD_N_BANDS; b++ ) {\r
51         psSilk_VAD->NoiseLevelBias[ b ] = SKP_max_32( SKP_DIV32_16( VAD_NOISE_LEVELS_BIAS, b + 1 ), 1 );\r
52     }\r
53 \r
54     /* Initialize state */\r
55     for( b = 0; b < VAD_N_BANDS; b++ ) {\r
56         psSilk_VAD->NL[ b ]     = SKP_MUL( 100, psSilk_VAD->NoiseLevelBias[ b ] );\r
57         psSilk_VAD->inv_NL[ b ] = SKP_DIV32( SKP_int32_MAX, psSilk_VAD->NL[ b ] );\r
58     }\r
59     psSilk_VAD->counter = 15;\r
60 \r
61     /* init smoothed energy-to-noise ratio*/\r
62     for( b = 0; b < VAD_N_BANDS; b++ ) {\r
63         psSilk_VAD->NrgRatioSmth_Q8[ b ] = 100 * 256;       /* 100 * 256 --> 20 dB SNR */\r
64     }\r
65 \r
66     return( ret );\r
67 }\r
68 \r
69 /* Weighting factors for tilt measure */\r
70 const static SKP_int32 tiltWeights[ VAD_N_BANDS ] = { 30000, 6000, -12000, -12000 };\r
71 \r
72 /***************************************/\r
73 /* Get the speech activity level in Q8 */\r
74 /***************************************/\r
75 SKP_int SKP_Silk_VAD_GetSA_Q8(                                      /* O    Return value, 0 if success      */\r
76     SKP_Silk_VAD_state              *psSilk_VAD,                    /* I/O  Silk VAD state                  */\r
77     SKP_int                         *pSA_Q8,                        /* O    Speech activity level in Q8     */\r
78     SKP_int                         *pSNR_dB_Q7,                    /* O    SNR for current frame in Q7     */\r
79     SKP_int                         pQuality_Q15[ VAD_N_BANDS ],    /* O    Smoothed SNR for each band      */\r
80     SKP_int                         *pTilt_Q15,                     /* O    current frame's frequency tilt  */\r
81     const SKP_int16                 pIn[],                          /* I    PCM input       [framelength]   */\r
82     const SKP_int                   framelength                     /* I    Input frame length              */\r
83 )\r
84 {\r
85     SKP_int   SA_Q15, input_tilt;\r
86     SKP_int32 scratch[ 3 * MAX_FRAME_LENGTH / 2 ];\r
87     SKP_int   decimated_framelength, dec_subframe_length, dec_subframe_offset, SNR_Q7, i, b, s;\r
88     SKP_int32 sumSquared, smooth_coef_Q16;\r
89     SKP_int16 HPstateTmp;\r
90 \r
91     SKP_int16 X[ VAD_N_BANDS ][ MAX_FRAME_LENGTH / 2 ];\r
92     SKP_int32 Xnrg[ VAD_N_BANDS ];\r
93     SKP_int32 NrgToNoiseRatio_Q8[ VAD_N_BANDS ];\r
94     SKP_int32 speech_nrg, x_tmp;\r
95     SKP_int   ret = 0;\r
96 \r
97     /* Safety checks */\r
98     SKP_assert( VAD_N_BANDS == 4 );\r
99     SKP_assert( MAX_FRAME_LENGTH >= framelength );\r
100     SKP_assert( framelength <= 512 );\r
101 \r
102     /***********************/\r
103     /* Filter and Decimate */\r
104     /***********************/\r
105     /* 0-8 kHz to 0-4 kHz and 4-8 kHz */\r
106     SKP_Silk_ana_filt_bank_1( pIn,          &psSilk_VAD->AnaState[  0 ], &X[ 0 ][ 0 ], &X[ 3 ][ 0 ], &scratch[ 0 ], framelength );        \r
107     \r
108     /* 0-4 kHz to 0-2 kHz and 2-4 kHz */\r
109     SKP_Silk_ana_filt_bank_1( &X[ 0 ][ 0 ], &psSilk_VAD->AnaState1[ 0 ], &X[ 0 ][ 0 ], &X[ 2 ][ 0 ], &scratch[ 0 ], SKP_RSHIFT( framelength, 1 ) );\r
110     \r
111     /* 0-2 kHz to 0-1 kHz and 1-2 kHz */\r
112     SKP_Silk_ana_filt_bank_1( &X[ 0 ][ 0 ], &psSilk_VAD->AnaState2[ 0 ], &X[ 0 ][ 0 ], &X[ 1 ][ 0 ], &scratch[ 0 ], SKP_RSHIFT( framelength, 2 ) );\r
113 \r
114     /*********************************************/\r
115     /* HP filter on lowest band (differentiator) */\r
116     /*********************************************/\r
117     decimated_framelength = SKP_RSHIFT( framelength, 3 );\r
118     X[ 0 ][ decimated_framelength - 1 ] = SKP_RSHIFT( X[ 0 ][ decimated_framelength - 1 ], 1 );\r
119     HPstateTmp = X[ 0 ][ decimated_framelength - 1 ];\r
120     for( i = decimated_framelength - 1; i > 0; i-- ) {\r
121         X[ 0 ][ i - 1 ]  = SKP_RSHIFT( X[ 0 ][ i - 1 ], 1 );\r
122         X[ 0 ][ i ]     -= X[ 0 ][ i - 1 ];\r
123     }\r
124     X[ 0 ][ 0 ] -= psSilk_VAD->HPstate;\r
125     psSilk_VAD->HPstate = HPstateTmp;\r
126 \r
127     /*************************************/\r
128     /* Calculate the energy in each band */\r
129     /*************************************/\r
130     for( b = 0; b < VAD_N_BANDS; b++ ) {        \r
131         /* Find the decimated framelength in the non-uniformly divided bands */\r
132         decimated_framelength = SKP_RSHIFT( framelength, SKP_min_int( VAD_N_BANDS - b, VAD_N_BANDS - 1 ) );\r
133 \r
134         /* Split length into subframe lengths */\r
135         dec_subframe_length = SKP_RSHIFT( decimated_framelength, VAD_INTERNAL_SUBFRAMES_LOG2 );\r
136         dec_subframe_offset = 0;\r
137 \r
138         /* Compute energy per sub-frame */\r
139         /* initialize with summed energy of last subframe */\r
140         Xnrg[ b ] = psSilk_VAD->XnrgSubfr[ b ];\r
141         for( s = 0; s < VAD_INTERNAL_SUBFRAMES; s++ ) {\r
142             sumSquared = 0;\r
143             for( i = 0; i < dec_subframe_length; i++ ) {\r
144                 /* The energy will be less than dec_subframe_length * ( SKP_int16_MIN / 8 )^2.              */\r
145                 /* Therefore we can accumulate with no risk of overflow (unless dec_subframe_length > 128)  */\r
146                 x_tmp = SKP_RSHIFT( X[ b ][ i + dec_subframe_offset ], 3 );\r
147                 sumSquared = SKP_SMLABB( sumSquared, x_tmp, x_tmp );\r
148 \r
149                 /* Safety check */\r
150                 SKP_assert( sumSquared >= 0 );\r
151             }\r
152 \r
153             /* add/saturate summed energy of current subframe */\r
154             if( s < VAD_INTERNAL_SUBFRAMES - 1 ) {\r
155                 Xnrg[ b ] = SKP_ADD_POS_SAT32( Xnrg[ b ], sumSquared );\r
156             } else {\r
157                 /* look-ahead subframe */\r
158                 Xnrg[ b ] = SKP_ADD_POS_SAT32( Xnrg[ b ], SKP_RSHIFT( sumSquared, 1 ) );\r
159             }\r
160 \r
161             dec_subframe_offset += dec_subframe_length;\r
162         }\r
163         psSilk_VAD->XnrgSubfr[ b ] = sumSquared; \r
164     }\r
165 \r
166     /********************/\r
167     /* Noise estimation */\r
168     /********************/\r
169     SKP_Silk_VAD_GetNoiseLevels( &Xnrg[ 0 ], psSilk_VAD );\r
170 \r
171     /***********************************************/\r
172     /* Signal-plus-noise to noise ratio estimation */\r
173     /***********************************************/\r
174     sumSquared = 0;\r
175     input_tilt = 0;\r
176     for( b = 0; b < VAD_N_BANDS; b++ ) {\r
177         speech_nrg = Xnrg[ b ] - psSilk_VAD->NL[ b ];\r
178         if( speech_nrg > 0 ) {\r
179             /* Divide, with sufficient resolution */\r
180             if( ( Xnrg[ b ] & 0xFF800000 ) == 0 ) {\r
181                 NrgToNoiseRatio_Q8[ b ] = SKP_DIV32( SKP_LSHIFT( Xnrg[ b ], 8 ), psSilk_VAD->NL[ b ] + 1 );\r
182             } else {\r
183                 NrgToNoiseRatio_Q8[ b ] = SKP_DIV32( Xnrg[ b ], SKP_RSHIFT( psSilk_VAD->NL[ b ], 8 ) + 1 );\r
184             }\r
185 \r
186             /* Convert to log domain */\r
187             SNR_Q7 = SKP_Silk_lin2log( NrgToNoiseRatio_Q8[ b ] ) - 8 * 128;\r
188 \r
189             /* Sum-of-squares */\r
190             sumSquared = SKP_SMLABB( sumSquared, SNR_Q7, SNR_Q7 );          /* Q14 */\r
191 \r
192             /* Tilt measure */\r
193             if( speech_nrg < ( 1 << 20 ) ) {\r
194                 /* Scale down SNR value for small subband speech energies */\r
195                 SNR_Q7 = SKP_SMULWB( SKP_LSHIFT( SKP_Silk_SQRT_APPROX( speech_nrg ), 6 ), SNR_Q7 );\r
196             }\r
197             input_tilt = SKP_SMLAWB( input_tilt, tiltWeights[ b ], SNR_Q7 );\r
198         } else {\r
199             NrgToNoiseRatio_Q8[ b ] = 256;\r
200         }\r
201     }\r
202 \r
203     /* Mean-of-squares */\r
204     sumSquared = SKP_DIV32_16( sumSquared, VAD_N_BANDS );           /* Q14 */\r
205 \r
206     /* Root-mean-square approximation, scale to dBs, and write to output pointer */\r
207     *pSNR_dB_Q7 = ( SKP_int16 )( 3 * SKP_Silk_SQRT_APPROX( sumSquared ) );  /* Q7 */\r
208 \r
209     /*********************************/\r
210     /* Speech Probability Estimation */\r
211     /*********************************/\r
212     SA_Q15 = SKP_Silk_sigm_Q15( SKP_SMULWB( VAD_SNR_FACTOR_Q16, *pSNR_dB_Q7 ) - VAD_NEGATIVE_OFFSET_Q5 );\r
213 \r
214     /**************************/\r
215     /* Frequency Tilt Measure */\r
216     /**************************/\r
217     *pTilt_Q15 = SKP_LSHIFT( SKP_Silk_sigm_Q15( input_tilt ) - 16384, 1 );\r
218 \r
219     /**************************************************/\r
220     /* Scale the sigmoid output based on power levels */\r
221     /**************************************************/\r
222     speech_nrg = 0;\r
223     for( b = 0; b < VAD_N_BANDS; b++ ) {\r
224         /* Accumulate signal-without-noise energies, higher frequency bands have more weight */\r
225         speech_nrg += ( b + 1 ) * SKP_RSHIFT( Xnrg[ b ] - psSilk_VAD->NL[ b ], 4 );\r
226     }\r
227 \r
228     /* Power scaling */\r
229     if( speech_nrg <= 0 ) {\r
230         SA_Q15 = SKP_RSHIFT( SA_Q15, 1 ); \r
231     } else if( speech_nrg < 32768 ) {\r
232         /* square-root */\r
233         speech_nrg = SKP_Silk_SQRT_APPROX( SKP_LSHIFT( speech_nrg, 15 ) );\r
234         SA_Q15 = SKP_SMULWB( 32768 + speech_nrg, SA_Q15 ); \r
235     }\r
236 \r
237     /* Copy the resulting speech activity in Q8 to *pSA_Q8 */\r
238     *pSA_Q8 = SKP_min_int( SKP_RSHIFT( SA_Q15, 7 ), SKP_uint8_MAX );\r
239 \r
240     /***********************************/\r
241     /* Energy Level and SNR estimation */\r
242     /***********************************/\r
243     /* smoothing coefficient */\r
244     smooth_coef_Q16 = SKP_SMULWB( VAD_SNR_SMOOTH_COEF_Q18, SKP_SMULWB( SA_Q15, SA_Q15 ) );\r
245     for( b = 0; b < VAD_N_BANDS; b++ ) {\r
246         /* compute smoothed energy-to-noise ratio per band */\r
247         psSilk_VAD->NrgRatioSmth_Q8[ b ] = SKP_SMLAWB( psSilk_VAD->NrgRatioSmth_Q8[ b ], \r
248             NrgToNoiseRatio_Q8[ b ] - psSilk_VAD->NrgRatioSmth_Q8[ b ], smooth_coef_Q16 );\r
249 \r
250         /* signal to noise ratio in dB per band */\r
251         SNR_Q7 = 3 * ( SKP_Silk_lin2log( psSilk_VAD->NrgRatioSmth_Q8[b] ) - 8 * 128 );\r
252         /* quality = sigmoid( 0.25 * ( SNR_dB - 16 ) ); */\r
253         pQuality_Q15[ b ] = SKP_Silk_sigm_Q15( SKP_RSHIFT( SNR_Q7 - 16 * 128, 4 ) );\r
254     }\r
255 \r
256     return( ret );\r
257 }\r
258 \r
259 /**************************/\r
260 /* Noise level estimation */\r
261 /**************************/\r
262 void SKP_Silk_VAD_GetNoiseLevels(\r
263     const SKP_int32                 pX[ VAD_N_BANDS ],  /* I    subband energies                            */\r
264     SKP_Silk_VAD_state              *psSilk_VAD         /* I/O  Pointer to Silk VAD state                   */ \r
265 )\r
266 {\r
267     SKP_int   k;\r
268     SKP_int32 nl, nrg, inv_nrg;\r
269     SKP_int   coef, min_coef;\r
270 \r
271     /* Initially faster smoothing */\r
272     if( psSilk_VAD->counter < 1000 ) { /* 1000 = 20 sec */\r
273         min_coef = SKP_DIV32_16( SKP_int16_MAX, SKP_RSHIFT( psSilk_VAD->counter, 4 ) + 1 );  \r
274     } else {\r
275         min_coef = 0;\r
276     }\r
277 \r
278     for( k = 0; k < VAD_N_BANDS; k++ ) {\r
279         /* Get old noise level estimate for current band */\r
280         nl = psSilk_VAD->NL[ k ];\r
281         SKP_assert( nl >= 0 );\r
282         \r
283         /* Add bias */\r
284         nrg = SKP_ADD_POS_SAT32( pX[ k ], psSilk_VAD->NoiseLevelBias[ k ] ); \r
285         SKP_assert( nrg > 0 );\r
286         \r
287         /* Invert energies */\r
288         inv_nrg = SKP_DIV32( SKP_int32_MAX, nrg );\r
289         SKP_assert( inv_nrg >= 0 );\r
290         \r
291         /* Less update when subband energy is high */\r
292         if( nrg > SKP_LSHIFT( nl, 3 ) ) {\r
293             coef = VAD_NOISE_LEVEL_SMOOTH_COEF_Q16 >> 3;\r
294         } else if( nrg < nl ) {\r
295             coef = VAD_NOISE_LEVEL_SMOOTH_COEF_Q16;\r
296         } else {\r
297             coef = SKP_SMULWB( SKP_SMULWW( inv_nrg, nl ), VAD_NOISE_LEVEL_SMOOTH_COEF_Q16 << 1 );\r
298         }\r
299 \r
300         /* Initially faster smoothing */\r
301         coef = SKP_max_int( coef, min_coef );\r
302 \r
303         /* Smooth inverse energies */\r
304         psSilk_VAD->inv_NL[ k ] = SKP_SMLAWB( psSilk_VAD->inv_NL[ k ], inv_nrg - psSilk_VAD->inv_NL[ k ], coef );\r
305         SKP_assert( psSilk_VAD->inv_NL[ k ] >= 0 );\r
306 \r
307         /* Compute noise level by inverting again */\r
308         nl = SKP_DIV32( SKP_int32_MAX, psSilk_VAD->inv_NL[ k ] );\r
309         SKP_assert( nl >= 0 );\r
310 \r
311         /* Limit noise levels (guarantee 7 bits of head room) */\r
312         nl = SKP_min( nl, 0x00FFFFFF );\r
313 \r
314         /* Store as part of state */\r
315         psSilk_VAD->NL[ k ] = nl;\r
316     }\r
317 \r
318     /* Increment frame counter */\r
319     psSilk_VAD->counter++;\r
320 }\r