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