Fix stack corruption with high rate SILK encoding added by the recent CBR changes...
[opus.git] / silk / float / noise_shape_analysis_FLP.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 "main_FLP.h"
33 #include "tuning_parameters.h"
34
35 /* Compute gain to make warped filter coefficients have a zero mean log frequency response on a     */
36 /* non-warped frequency scale. (So that it can be implemented with a minimum-phase monic filter.)   */
37 static inline silk_float warped_gain(
38     const silk_float     *coefs,
39     silk_float           lambda,
40     opus_int             order
41 ) {
42     opus_int   i;
43     silk_float gain;
44
45     lambda = -lambda;
46     gain = coefs[ order - 1 ];
47     for( i = order - 2; i >= 0; i-- ) {
48         gain = lambda * gain + coefs[ i ];
49     }
50     return (silk_float)( 1.0f / ( 1.0f - lambda * gain ) );
51 }
52
53 /* Convert warped filter coefficients to monic pseudo-warped coefficients and limit maximum     */
54 /* amplitude of monic warped coefficients by using bandwidth expansion on the true coefficients */
55 static inline void warped_true2monic_coefs(
56     silk_float           *coefs_syn,
57     silk_float           *coefs_ana,
58     silk_float           lambda,
59     silk_float           limit,
60     opus_int             order
61 ) {
62     opus_int   i, iter, ind = 0;
63     silk_float tmp, maxabs, chirp, gain_syn, gain_ana;
64
65     /* Convert to monic coefficients */
66     for( i = order - 1; i > 0; i-- ) {
67         coefs_syn[ i - 1 ] -= lambda * coefs_syn[ i ];
68         coefs_ana[ i - 1 ] -= lambda * coefs_ana[ i ];
69     }
70     gain_syn = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_syn[ 0 ] );
71     gain_ana = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_ana[ 0 ] );
72     for( i = 0; i < order; i++ ) {
73         coefs_syn[ i ] *= gain_syn;
74         coefs_ana[ i ] *= gain_ana;
75     }
76
77     /* Limit */
78     for( iter = 0; iter < 10; iter++ ) {
79         /* Find maximum absolute value */
80         maxabs = -1.0f;
81         for( i = 0; i < order; i++ ) {
82             tmp = silk_max( silk_abs_float( coefs_syn[ i ] ), silk_abs_float( coefs_ana[ i ] ) );
83             if( tmp > maxabs ) {
84                 maxabs = tmp;
85                 ind = i;
86             }
87         }
88         if( maxabs <= limit ) {
89             /* Coefficients are within range - done */
90             return;
91         }
92
93         /* Convert back to true warped coefficients */
94         for( i = 1; i < order; i++ ) {
95             coefs_syn[ i - 1 ] += lambda * coefs_syn[ i ];
96             coefs_ana[ i - 1 ] += lambda * coefs_ana[ i ];
97         }
98         gain_syn = 1.0f / gain_syn;
99         gain_ana = 1.0f / gain_ana;
100         for( i = 0; i < order; i++ ) {
101             coefs_syn[ i ] *= gain_syn;
102             coefs_ana[ i ] *= gain_ana;
103         }
104
105         /* Apply bandwidth expansion */
106         chirp = 0.99f - ( 0.8f + 0.1f * iter ) * ( maxabs - limit ) / ( maxabs * ( ind + 1 ) );
107         silk_bwexpander_FLP( coefs_syn, order, chirp );
108         silk_bwexpander_FLP( coefs_ana, order, chirp );
109
110         /* Convert to monic warped coefficients */
111         for( i = order - 1; i > 0; i-- ) {
112             coefs_syn[ i - 1 ] -= lambda * coefs_syn[ i ];
113             coefs_ana[ i - 1 ] -= lambda * coefs_ana[ i ];
114         }
115         gain_syn = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_syn[ 0 ] );
116         gain_ana = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_ana[ 0 ] );
117         for( i = 0; i < order; i++ ) {
118             coefs_syn[ i ] *= gain_syn;
119             coefs_ana[ i ] *= gain_ana;
120         }
121     }
122     silk_assert( 0 );
123 }
124
125 /* Compute noise shaping coefficients and initial gain values */
126 void silk_noise_shape_analysis_FLP(
127     silk_encoder_state_FLP          *psEnc,             /* I/O  Encoder state FLP                       */
128     silk_encoder_control_FLP        *psEncCtrl,         /* I/O  Encoder control FLP                     */
129     const silk_float                 *pitch_res,         /* I    LPC residual from pitch analysis        */
130     const silk_float                 *x                  /* I    Input signal [frame_length + la_shape]  */
131 )
132 {
133     silk_shape_state_FLP *psShapeSt = &psEnc->sShape;
134     opus_int     k, nSamples;
135     silk_float   SNR_adj_dB, HarmBoost, HarmShapeGain, Tilt;
136     silk_float   nrg, pre_nrg, log_energy, log_energy_prev, energy_variation;
137     silk_float   delta, BWExp1, BWExp2, gain_mult, gain_add, strength, b, warping;
138     silk_float   x_windowed[ SHAPE_LPC_WIN_MAX ];
139     silk_float   auto_corr[ MAX_SHAPE_LPC_ORDER + 1 ];
140     const silk_float *x_ptr, *pitch_res_ptr;
141
142     /* Point to start of first LPC analysis block */
143     x_ptr = x - psEnc->sCmn.la_shape;
144
145     /****************/
146     /* GAIN CONTROL */
147     /****************/
148     SNR_adj_dB = psEnc->sCmn.SNR_dB_Q7 * ( 1 / 128.0f );
149
150     /* Input quality is the average of the quality in the lowest two VAD bands */
151     psEncCtrl->input_quality = 0.5f * ( psEnc->sCmn.input_quality_bands_Q15[ 0 ] + psEnc->sCmn.input_quality_bands_Q15[ 1 ] ) * ( 1.0f / 32768.0f );
152
153     /* Coding quality level, between 0.0 and 1.0 */
154     psEncCtrl->coding_quality = silk_sigmoid( 0.25f * ( SNR_adj_dB - 18.0f ) );
155
156     if( psEnc->sCmn.useCBR == 0 ) {
157         /* Reduce coding SNR during low speech activity */
158         b = 1.0f - psEnc->sCmn.speech_activity_Q8 * ( 1.0f /  256.0f );
159         SNR_adj_dB -= BG_SNR_DECR_dB * psEncCtrl->coding_quality * ( 0.5f + 0.5f * psEncCtrl->input_quality ) * b * b;
160     }
161
162     if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
163         /* Reduce gains for periodic signals */
164         SNR_adj_dB += HARM_SNR_INCR_dB * psEnc->LTPCorr;
165     } else {
166         /* For unvoiced signals and low-quality input, adjust the quality slower than SNR_dB setting */
167         SNR_adj_dB += ( -0.4f * psEnc->sCmn.SNR_dB_Q7 * ( 1 / 128.0f ) + 6.0f ) * ( 1.0f - psEncCtrl->input_quality );
168     }
169
170     /*************************/
171     /* SPARSENESS PROCESSING */
172     /*************************/
173     /* Set quantizer offset */
174     if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
175         /* Initally set to 0; may be overruled in process_gains(..) */
176         psEnc->sCmn.indices.quantOffsetType = 0;
177         psEncCtrl->sparseness = 0.0f;
178     } else {
179         /* Sparseness measure, based on relative fluctuations of energy per 2 milliseconds */
180         nSamples = 2 * psEnc->sCmn.fs_kHz;
181         energy_variation = 0.0f;
182         log_energy_prev  = 0.0f;
183         pitch_res_ptr = pitch_res;
184         for( k = 0; k < silk_SMULBB( SUB_FRAME_LENGTH_MS, psEnc->sCmn.nb_subfr ) / 2; k++ ) {
185             nrg = ( silk_float )nSamples + ( silk_float )silk_energy_FLP( pitch_res_ptr, nSamples );
186             log_energy = silk_log2( nrg );
187             if( k > 0 ) {
188                 energy_variation += silk_abs_float( log_energy - log_energy_prev );
189             }
190             log_energy_prev = log_energy;
191             pitch_res_ptr += nSamples;
192         }
193         psEncCtrl->sparseness = silk_sigmoid( 0.4f * ( energy_variation - 5.0f ) );
194
195         /* Set quantization offset depending on sparseness measure */
196         if( psEncCtrl->sparseness > SPARSENESS_THRESHOLD_QNT_OFFSET ) {
197             psEnc->sCmn.indices.quantOffsetType = 0;
198         } else {
199             psEnc->sCmn.indices.quantOffsetType = 1;
200         }
201
202         /* Increase coding SNR for sparse signals */
203         SNR_adj_dB += SPARSE_SNR_INCR_dB * ( psEncCtrl->sparseness - 0.5f );
204     }
205
206     /*******************************/
207     /* Control bandwidth expansion */
208     /*******************************/
209     /* More BWE for signals with high prediction gain */
210     strength = FIND_PITCH_WHITE_NOISE_FRACTION * psEncCtrl->predGain;           /* between 0.0 and 1.0 */
211     BWExp1 = BWExp2 = BANDWIDTH_EXPANSION / ( 1.0f + strength * strength );
212     delta  = LOW_RATE_BANDWIDTH_EXPANSION_DELTA * ( 1.0f - 0.75f * psEncCtrl->coding_quality );
213     BWExp1 -= delta;
214     BWExp2 += delta;
215     /* BWExp1 will be applied after BWExp2, so make it relative */
216     BWExp1 /= BWExp2;
217
218     if( psEnc->sCmn.warping_Q16 > 0 ) {
219         /* Slightly more warping in analysis will move quantization noise up in frequency, where it's better masked */
220         warping = (silk_float)psEnc->sCmn.warping_Q16 / 65536.0f + 0.01f * psEncCtrl->coding_quality;
221     } else {
222         warping = 0.0f;
223     }
224
225     /********************************************/
226     /* Compute noise shaping AR coefs and gains */
227     /********************************************/
228     for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
229         /* Apply window: sine slope followed by flat part followed by cosine slope */
230         opus_int shift, slope_part, flat_part;
231         flat_part = psEnc->sCmn.fs_kHz * 3;
232         slope_part = ( psEnc->sCmn.shapeWinLength - flat_part ) / 2;
233
234         silk_apply_sine_window_FLP( x_windowed, x_ptr, 1, slope_part );
235         shift = slope_part;
236         silk_memcpy( x_windowed + shift, x_ptr + shift, flat_part * sizeof(silk_float) );
237         shift += flat_part;
238         silk_apply_sine_window_FLP( x_windowed + shift, x_ptr + shift, 2, slope_part );
239
240         /* Update pointer: next LPC analysis block */
241         x_ptr += psEnc->sCmn.subfr_length;
242
243         if( psEnc->sCmn.warping_Q16 > 0 ) {
244             /* Calculate warped auto correlation */
245             silk_warped_autocorrelation_FLP( auto_corr, x_windowed, warping,
246                 psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder );
247         } else {
248             /* Calculate regular auto correlation */
249             silk_autocorrelation_FLP( auto_corr, x_windowed, psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder + 1 );
250         }
251
252         /* Add white noise, as a fraction of energy */
253         auto_corr[ 0 ] += auto_corr[ 0 ] * SHAPE_WHITE_NOISE_FRACTION;
254
255         /* Convert correlations to prediction coefficients, and compute residual energy */
256         nrg = silk_levinsondurbin_FLP( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], auto_corr, psEnc->sCmn.shapingLPCOrder );
257         psEncCtrl->Gains[ k ] = ( silk_float )sqrt( nrg );
258
259         if( psEnc->sCmn.warping_Q16 > 0 ) {
260             /* Adjust gain for warping */
261             psEncCtrl->Gains[ k ] *= warped_gain( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], warping, psEnc->sCmn.shapingLPCOrder );
262         }
263
264         /* Bandwidth expansion for synthesis filter shaping */
265         silk_bwexpander_FLP( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder, BWExp2 );
266
267         /* Compute noise shaping filter coefficients */
268         silk_memcpy(
269             &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ],
270             &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ],
271             psEnc->sCmn.shapingLPCOrder * sizeof( silk_float ) );
272
273         /* Bandwidth expansion for analysis filter shaping */
274         silk_bwexpander_FLP( &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder, BWExp1 );
275
276         /* Ratio of prediction gains, in energy domain */
277         silk_LPC_inverse_pred_gain_FLP( &pre_nrg, &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder );
278         silk_LPC_inverse_pred_gain_FLP( &nrg,     &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder );
279         psEncCtrl->GainsPre[ k ] = 1.0f - 0.7f * ( 1.0f - pre_nrg / nrg );
280
281         /* Convert to monic warped prediction coefficients and limit absolute values */
282         warped_true2monic_coefs( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ],
283             warping, 3.999f, psEnc->sCmn.shapingLPCOrder );
284     }
285
286     /*****************/
287     /* Gain tweaking */
288     /*****************/
289     /* Increase gains during low speech activity */
290     gain_mult = (silk_float)pow( 2.0f, -0.16f * SNR_adj_dB );
291     gain_add  = (silk_float)pow( 2.0f,  0.16f * MIN_QGAIN_DB );
292     for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
293         psEncCtrl->Gains[ k ] *= gain_mult;
294         psEncCtrl->Gains[ k ] += gain_add;
295     }
296
297     gain_mult = 1.0f + INPUT_TILT + psEncCtrl->coding_quality * HIGH_RATE_INPUT_TILT;
298     for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
299         psEncCtrl->GainsPre[ k ] *= gain_mult;
300     }
301
302     /************************************************/
303     /* Control low-frequency shaping and noise tilt */
304     /************************************************/
305     /* Less low frequency shaping for noisy inputs */
306     strength = LOW_FREQ_SHAPING * ( 1.0f + LOW_QUALITY_LOW_FREQ_SHAPING_DECR * ( psEnc->sCmn.input_quality_bands_Q15[ 0 ] * ( 1.0f / 32768.0f ) - 1.0f ) );
307     strength *= psEnc->sCmn.speech_activity_Q8 * ( 1.0f /  256.0f );
308     if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
309         /* Reduce low frequencies quantization noise for periodic signals, depending on pitch lag */
310         /*f = 400; freqz([1, -0.98 + 2e-4 * f], [1, -0.97 + 7e-4 * f], 2^12, Fs); axis([0, 1000, -10, 1])*/
311         for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
312             b = 0.2f / psEnc->sCmn.fs_kHz + 3.0f / psEncCtrl->pitchL[ k ];
313             psEncCtrl->LF_MA_shp[ k ] = -1.0f + b;
314             psEncCtrl->LF_AR_shp[ k ] =  1.0f - b - b * strength;
315         }
316         Tilt = - HP_NOISE_COEF -
317             (1 - HP_NOISE_COEF) * HARM_HP_NOISE_COEF * psEnc->sCmn.speech_activity_Q8 * ( 1.0f /  256.0f );
318     } else {
319         b = 1.3f / psEnc->sCmn.fs_kHz;
320         psEncCtrl->LF_MA_shp[ 0 ] = -1.0f + b;
321         psEncCtrl->LF_AR_shp[ 0 ] =  1.0f - b - b * strength * 0.6f;
322         for( k = 1; k < psEnc->sCmn.nb_subfr; k++ ) {
323             psEncCtrl->LF_MA_shp[ k ] = psEncCtrl->LF_MA_shp[ 0 ];
324             psEncCtrl->LF_AR_shp[ k ] = psEncCtrl->LF_AR_shp[ 0 ];
325         }
326         Tilt = -HP_NOISE_COEF;
327     }
328
329     /****************************/
330     /* HARMONIC SHAPING CONTROL */
331     /****************************/
332     /* Control boosting of harmonic frequencies */
333     HarmBoost = LOW_RATE_HARMONIC_BOOST * ( 1.0f - psEncCtrl->coding_quality ) * psEnc->LTPCorr;
334
335     /* More harmonic boost for noisy input signals */
336     HarmBoost += LOW_INPUT_QUALITY_HARMONIC_BOOST * ( 1.0f - psEncCtrl->input_quality );
337
338     if( USE_HARM_SHAPING && psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
339         /* Harmonic noise shaping */
340         HarmShapeGain = HARMONIC_SHAPING;
341
342         /* More harmonic noise shaping for high bitrates or noisy input */
343         HarmShapeGain += HIGH_RATE_OR_LOW_QUALITY_HARMONIC_SHAPING *
344             ( 1.0f - ( 1.0f - psEncCtrl->coding_quality ) * psEncCtrl->input_quality );
345
346         /* Less harmonic noise shaping for less periodic signals */
347         HarmShapeGain *= ( silk_float )sqrt( psEnc->LTPCorr );
348     } else {
349         HarmShapeGain = 0.0f;
350     }
351
352     /*************************/
353     /* Smooth over subframes */
354     /*************************/
355     for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
356         psShapeSt->HarmBoost_smth     += SUBFR_SMTH_COEF * ( HarmBoost - psShapeSt->HarmBoost_smth );
357         psEncCtrl->HarmBoost[ k ]      = psShapeSt->HarmBoost_smth;
358         psShapeSt->HarmShapeGain_smth += SUBFR_SMTH_COEF * ( HarmShapeGain - psShapeSt->HarmShapeGain_smth );
359         psEncCtrl->HarmShapeGain[ k ]  = psShapeSt->HarmShapeGain_smth;
360         psShapeSt->Tilt_smth          += SUBFR_SMTH_COEF * ( Tilt - psShapeSt->Tilt_smth );
361         psEncCtrl->Tilt[ k ]           = psShapeSt->Tilt_smth;
362     }
363 }