1/***********************************************************************
2Copyright (c) 2006-2011, Skype Limited. All rights reserved.
3Redistribution and use in source and binary forms, with or without
4modification, are permitted provided that the following conditions
5are met:
6- Redistributions of source code must retain the above copyright notice,
7this list of conditions and the following disclaimer.
8- Redistributions in binary form must reproduce the above copyright
9notice, this list of conditions and the following disclaimer in the
10documentation and/or other materials provided with the distribution.
11- Neither the name of Internet Society, IETF or IETF Trust, nor the
12names of specific contributors, may be used to endorse or promote
13products derived from this software without specific prior written
14permission.
15THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
19LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
20CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
21SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
23CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
24ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
25POSSIBILITY OF SUCH DAMAGE.
26***********************************************************************/
27
28#ifdef HAVE_CONFIG_H
29#include "config.h"
30#endif
31
32#include "main_FIX.h"
33#include "stack_alloc.h"
34#include "tuning_parameters.h"
35
36/* Compute gain to make warped filter coefficients have a zero mean log frequency response on a   */
37/* non-warped frequency scale. (So that it can be implemented with a minimum-phase monic filter.) */
38/* Note: A monic filter is one with the first coefficient equal to 1.0. In Silk we omit the first */
39/* coefficient in an array of coefficients, for monic filters.                                    */
40static OPUS_INLINE opus_int32 warped_gain( /* gain in Q16*/
41    const opus_int32     *coefs_Q24,
42    opus_int             lambda_Q16,
43    opus_int             order
44) {
45    opus_int   i;
46    opus_int32 gain_Q24;
47
48    lambda_Q16 = -lambda_Q16;
49    gain_Q24 = coefs_Q24[ order - 1 ];
50    for( i = order - 2; i >= 0; i-- ) {
51        gain_Q24 = silk_SMLAWB( coefs_Q24[ i ], gain_Q24, lambda_Q16 );
52    }
53    gain_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), gain_Q24, -lambda_Q16 );
54    return silk_INVERSE32_varQ( gain_Q24, 40 );
55}
56
57/* Convert warped filter coefficients to monic pseudo-warped coefficients and limit maximum     */
58/* amplitude of monic warped coefficients by using bandwidth expansion on the true coefficients */
59static OPUS_INLINE void limit_warped_coefs(
60    opus_int32           *coefs_syn_Q24,
61    opus_int32           *coefs_ana_Q24,
62    opus_int             lambda_Q16,
63    opus_int32           limit_Q24,
64    opus_int             order
65) {
66    opus_int   i, iter, ind = 0;
67    opus_int32 tmp, maxabs_Q24, chirp_Q16, gain_syn_Q16, gain_ana_Q16;
68    opus_int32 nom_Q16, den_Q24;
69
70    /* Convert to monic coefficients */
71    lambda_Q16 = -lambda_Q16;
72    for( i = order - 1; i > 0; i-- ) {
73        coefs_syn_Q24[ i - 1 ] = silk_SMLAWB( coefs_syn_Q24[ i - 1 ], coefs_syn_Q24[ i ], lambda_Q16 );
74        coefs_ana_Q24[ i - 1 ] = silk_SMLAWB( coefs_ana_Q24[ i - 1 ], coefs_ana_Q24[ i ], lambda_Q16 );
75    }
76    lambda_Q16 = -lambda_Q16;
77    nom_Q16  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 16 ), -(opus_int32)lambda_Q16,        lambda_Q16 );
78    den_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), coefs_syn_Q24[ 0 ], lambda_Q16 );
79    gain_syn_Q16 = silk_DIV32_varQ( nom_Q16, den_Q24, 24 );
80    den_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), coefs_ana_Q24[ 0 ], lambda_Q16 );
81    gain_ana_Q16 = silk_DIV32_varQ( nom_Q16, den_Q24, 24 );
82    for( i = 0; i < order; i++ ) {
83        coefs_syn_Q24[ i ] = silk_SMULWW( gain_syn_Q16, coefs_syn_Q24[ i ] );
84        coefs_ana_Q24[ i ] = silk_SMULWW( gain_ana_Q16, coefs_ana_Q24[ i ] );
85    }
86
87    for( iter = 0; iter < 10; iter++ ) {
88        /* Find maximum absolute value */
89        maxabs_Q24 = -1;
90        for( i = 0; i < order; i++ ) {
91            tmp = silk_max( silk_abs_int32( coefs_syn_Q24[ i ] ), silk_abs_int32( coefs_ana_Q24[ i ] ) );
92            if( tmp > maxabs_Q24 ) {
93                maxabs_Q24 = tmp;
94                ind = i;
95            }
96        }
97        if( maxabs_Q24 <= limit_Q24 ) {
98            /* Coefficients are within range - done */
99            return;
100        }
101
102        /* Convert back to true warped coefficients */
103        for( i = 1; i < order; i++ ) {
104            coefs_syn_Q24[ i - 1 ] = silk_SMLAWB( coefs_syn_Q24[ i - 1 ], coefs_syn_Q24[ i ], lambda_Q16 );
105            coefs_ana_Q24[ i - 1 ] = silk_SMLAWB( coefs_ana_Q24[ i - 1 ], coefs_ana_Q24[ i ], lambda_Q16 );
106        }
107        gain_syn_Q16 = silk_INVERSE32_varQ( gain_syn_Q16, 32 );
108        gain_ana_Q16 = silk_INVERSE32_varQ( gain_ana_Q16, 32 );
109        for( i = 0; i < order; i++ ) {
110            coefs_syn_Q24[ i ] = silk_SMULWW( gain_syn_Q16, coefs_syn_Q24[ i ] );
111            coefs_ana_Q24[ i ] = silk_SMULWW( gain_ana_Q16, coefs_ana_Q24[ i ] );
112        }
113
114        /* Apply bandwidth expansion */
115        chirp_Q16 = SILK_FIX_CONST( 0.99, 16 ) - silk_DIV32_varQ(
116            silk_SMULWB( maxabs_Q24 - limit_Q24, silk_SMLABB( SILK_FIX_CONST( 0.8, 10 ), SILK_FIX_CONST( 0.1, 10 ), iter ) ),
117            silk_MUL( maxabs_Q24, ind + 1 ), 22 );
118        silk_bwexpander_32( coefs_syn_Q24, order, chirp_Q16 );
119        silk_bwexpander_32( coefs_ana_Q24, order, chirp_Q16 );
120
121        /* Convert to monic warped coefficients */
122        lambda_Q16 = -lambda_Q16;
123        for( i = order - 1; i > 0; i-- ) {
124            coefs_syn_Q24[ i - 1 ] = silk_SMLAWB( coefs_syn_Q24[ i - 1 ], coefs_syn_Q24[ i ], lambda_Q16 );
125            coefs_ana_Q24[ i - 1 ] = silk_SMLAWB( coefs_ana_Q24[ i - 1 ], coefs_ana_Q24[ i ], lambda_Q16 );
126        }
127        lambda_Q16 = -lambda_Q16;
128        nom_Q16  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 16 ), -(opus_int32)lambda_Q16,        lambda_Q16 );
129        den_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), coefs_syn_Q24[ 0 ], lambda_Q16 );
130        gain_syn_Q16 = silk_DIV32_varQ( nom_Q16, den_Q24, 24 );
131        den_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), coefs_ana_Q24[ 0 ], lambda_Q16 );
132        gain_ana_Q16 = silk_DIV32_varQ( nom_Q16, den_Q24, 24 );
133        for( i = 0; i < order; i++ ) {
134            coefs_syn_Q24[ i ] = silk_SMULWW( gain_syn_Q16, coefs_syn_Q24[ i ] );
135            coefs_ana_Q24[ i ] = silk_SMULWW( gain_ana_Q16, coefs_ana_Q24[ i ] );
136        }
137    }
138    silk_assert( 0 );
139}
140
141/**************************************************************/
142/* Compute noise shaping coefficients and initial gain values */
143/**************************************************************/
144void silk_noise_shape_analysis_FIX(
145    silk_encoder_state_FIX          *psEnc,                                 /* I/O  Encoder state FIX                                                           */
146    silk_encoder_control_FIX        *psEncCtrl,                             /* I/O  Encoder control FIX                                                         */
147    const opus_int16                *pitch_res,                             /* I    LPC residual from pitch analysis                                            */
148    const opus_int16                *x,                                     /* I    Input signal [ frame_length + la_shape ]                                    */
149    int                              arch                                   /* I    Run-time architecture                                                       */
150)
151{
152    silk_shape_state_FIX *psShapeSt = &psEnc->sShape;
153    opus_int     k, i, nSamples, Qnrg, b_Q14, warping_Q16, scale = 0;
154    opus_int32   SNR_adj_dB_Q7, HarmBoost_Q16, HarmShapeGain_Q16, Tilt_Q16, tmp32;
155    opus_int32   nrg, pre_nrg_Q30, log_energy_Q7, log_energy_prev_Q7, energy_variation_Q7;
156    opus_int32   delta_Q16, BWExp1_Q16, BWExp2_Q16, gain_mult_Q16, gain_add_Q16, strength_Q16, b_Q8;
157    opus_int32   auto_corr[     MAX_SHAPE_LPC_ORDER + 1 ];
158    opus_int32   refl_coef_Q16[ MAX_SHAPE_LPC_ORDER ];
159    opus_int32   AR1_Q24[       MAX_SHAPE_LPC_ORDER ];
160    opus_int32   AR2_Q24[       MAX_SHAPE_LPC_ORDER ];
161    VARDECL( opus_int16, x_windowed );
162    const opus_int16 *x_ptr, *pitch_res_ptr;
163    SAVE_STACK;
164
165    /* Point to start of first LPC analysis block */
166    x_ptr = x - psEnc->sCmn.la_shape;
167
168    /****************/
169    /* GAIN CONTROL */
170    /****************/
171    SNR_adj_dB_Q7 = psEnc->sCmn.SNR_dB_Q7;
172
173    /* Input quality is the average of the quality in the lowest two VAD bands */
174    psEncCtrl->input_quality_Q14 = ( opus_int )silk_RSHIFT( (opus_int32)psEnc->sCmn.input_quality_bands_Q15[ 0 ]
175        + psEnc->sCmn.input_quality_bands_Q15[ 1 ], 2 );
176
177    /* Coding quality level, between 0.0_Q0 and 1.0_Q0, but in Q14 */
178    psEncCtrl->coding_quality_Q14 = silk_RSHIFT( silk_sigm_Q15( silk_RSHIFT_ROUND( SNR_adj_dB_Q7 -
179        SILK_FIX_CONST( 20.0, 7 ), 4 ) ), 1 );
180
181    /* Reduce coding SNR during low speech activity */
182    if( psEnc->sCmn.useCBR == 0 ) {
183        b_Q8 = SILK_FIX_CONST( 1.0, 8 ) - psEnc->sCmn.speech_activity_Q8;
184        b_Q8 = silk_SMULWB( silk_LSHIFT( b_Q8, 8 ), b_Q8 );
185        SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7,
186            silk_SMULBB( SILK_FIX_CONST( -BG_SNR_DECR_dB, 7 ) >> ( 4 + 1 ), b_Q8 ),                                       /* Q11*/
187            silk_SMULWB( SILK_FIX_CONST( 1.0, 14 ) + psEncCtrl->input_quality_Q14, psEncCtrl->coding_quality_Q14 ) );     /* Q12*/
188    }
189
190    if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
191        /* Reduce gains for periodic signals */
192        SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7, SILK_FIX_CONST( HARM_SNR_INCR_dB, 8 ), psEnc->LTPCorr_Q15 );
193    } else {
194        /* For unvoiced signals and low-quality input, adjust the quality slower than SNR_dB setting */
195        SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7,
196            silk_SMLAWB( SILK_FIX_CONST( 6.0, 9 ), -SILK_FIX_CONST( 0.4, 18 ), psEnc->sCmn.SNR_dB_Q7 ),
197            SILK_FIX_CONST( 1.0, 14 ) - psEncCtrl->input_quality_Q14 );
198    }
199
200    /*************************/
201    /* SPARSENESS PROCESSING */
202    /*************************/
203    /* Set quantizer offset */
204    if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
205        /* Initially set to 0; may be overruled in process_gains(..) */
206        psEnc->sCmn.indices.quantOffsetType = 0;
207        psEncCtrl->sparseness_Q8 = 0;
208    } else {
209        /* Sparseness measure, based on relative fluctuations of energy per 2 milliseconds */
210        nSamples = silk_LSHIFT( psEnc->sCmn.fs_kHz, 1 );
211        energy_variation_Q7 = 0;
212        log_energy_prev_Q7  = 0;
213        pitch_res_ptr = pitch_res;
214        for( k = 0; k < silk_SMULBB( SUB_FRAME_LENGTH_MS, psEnc->sCmn.nb_subfr ) / 2; k++ ) {
215            silk_sum_sqr_shift( &nrg, &scale, pitch_res_ptr, nSamples );
216            nrg += silk_RSHIFT( nSamples, scale );           /* Q(-scale)*/
217
218            log_energy_Q7 = silk_lin2log( nrg );
219            if( k > 0 ) {
220                energy_variation_Q7 += silk_abs( log_energy_Q7 - log_energy_prev_Q7 );
221            }
222            log_energy_prev_Q7 = log_energy_Q7;
223            pitch_res_ptr += nSamples;
224        }
225
226        psEncCtrl->sparseness_Q8 = silk_RSHIFT( silk_sigm_Q15( silk_SMULWB( energy_variation_Q7 -
227            SILK_FIX_CONST( 5.0, 7 ), SILK_FIX_CONST( 0.1, 16 ) ) ), 7 );
228
229        /* Set quantization offset depending on sparseness measure */
230        if( psEncCtrl->sparseness_Q8 > SILK_FIX_CONST( SPARSENESS_THRESHOLD_QNT_OFFSET, 8 ) ) {
231            psEnc->sCmn.indices.quantOffsetType = 0;
232        } else {
233            psEnc->sCmn.indices.quantOffsetType = 1;
234        }
235
236        /* Increase coding SNR for sparse signals */
237        SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7, SILK_FIX_CONST( SPARSE_SNR_INCR_dB, 15 ), psEncCtrl->sparseness_Q8 - SILK_FIX_CONST( 0.5, 8 ) );
238    }
239
240    /*******************************/
241    /* Control bandwidth expansion */
242    /*******************************/
243    /* More BWE for signals with high prediction gain */
244    strength_Q16 = silk_SMULWB( psEncCtrl->predGain_Q16, SILK_FIX_CONST( FIND_PITCH_WHITE_NOISE_FRACTION, 16 ) );
245    BWExp1_Q16 = BWExp2_Q16 = silk_DIV32_varQ( SILK_FIX_CONST( BANDWIDTH_EXPANSION, 16 ),
246        silk_SMLAWW( SILK_FIX_CONST( 1.0, 16 ), strength_Q16, strength_Q16 ), 16 );
247    delta_Q16  = silk_SMULWB( SILK_FIX_CONST( 1.0, 16 ) - silk_SMULBB( 3, psEncCtrl->coding_quality_Q14 ),
248        SILK_FIX_CONST( LOW_RATE_BANDWIDTH_EXPANSION_DELTA, 16 ) );
249    BWExp1_Q16 = silk_SUB32( BWExp1_Q16, delta_Q16 );
250    BWExp2_Q16 = silk_ADD32( BWExp2_Q16, delta_Q16 );
251    /* BWExp1 will be applied after BWExp2, so make it relative */
252    BWExp1_Q16 = silk_DIV32_16( silk_LSHIFT( BWExp1_Q16, 14 ), silk_RSHIFT( BWExp2_Q16, 2 ) );
253
254    if( psEnc->sCmn.warping_Q16 > 0 ) {
255        /* Slightly more warping in analysis will move quantization noise up in frequency, where it's better masked */
256        warping_Q16 = silk_SMLAWB( psEnc->sCmn.warping_Q16, (opus_int32)psEncCtrl->coding_quality_Q14, SILK_FIX_CONST( 0.01, 18 ) );
257    } else {
258        warping_Q16 = 0;
259    }
260
261    /********************************************/
262    /* Compute noise shaping AR coefs and gains */
263    /********************************************/
264    ALLOC( x_windowed, psEnc->sCmn.shapeWinLength, opus_int16 );
265    for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
266        /* Apply window: sine slope followed by flat part followed by cosine slope */
267        opus_int shift, slope_part, flat_part;
268        flat_part = psEnc->sCmn.fs_kHz * 3;
269        slope_part = silk_RSHIFT( psEnc->sCmn.shapeWinLength - flat_part, 1 );
270
271        silk_apply_sine_window( x_windowed, x_ptr, 1, slope_part );
272        shift = slope_part;
273        silk_memcpy( x_windowed + shift, x_ptr + shift, flat_part * sizeof(opus_int16) );
274        shift += flat_part;
275        silk_apply_sine_window( x_windowed + shift, x_ptr + shift, 2, slope_part );
276
277        /* Update pointer: next LPC analysis block */
278        x_ptr += psEnc->sCmn.subfr_length;
279
280        if( psEnc->sCmn.warping_Q16 > 0 ) {
281            /* Calculate warped auto correlation */
282            silk_warped_autocorrelation_FIX( auto_corr, &scale, x_windowed, warping_Q16, psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder );
283        } else {
284            /* Calculate regular auto correlation */
285            silk_autocorr( auto_corr, &scale, x_windowed, psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder + 1, arch );
286        }
287
288        /* Add white noise, as a fraction of energy */
289        auto_corr[0] = silk_ADD32( auto_corr[0], silk_max_32( silk_SMULWB( silk_RSHIFT( auto_corr[ 0 ], 4 ),
290            SILK_FIX_CONST( SHAPE_WHITE_NOISE_FRACTION, 20 ) ), 1 ) );
291
292        /* Calculate the reflection coefficients using schur */
293        nrg = silk_schur64( refl_coef_Q16, auto_corr, psEnc->sCmn.shapingLPCOrder );
294        silk_assert( nrg >= 0 );
295
296        /* Convert reflection coefficients to prediction coefficients */
297        silk_k2a_Q16( AR2_Q24, refl_coef_Q16, psEnc->sCmn.shapingLPCOrder );
298
299        Qnrg = -scale;          /* range: -12...30*/
300        silk_assert( Qnrg >= -12 );
301        silk_assert( Qnrg <=  30 );
302
303        /* Make sure that Qnrg is an even number */
304        if( Qnrg & 1 ) {
305            Qnrg -= 1;
306            nrg >>= 1;
307        }
308
309        tmp32 = silk_SQRT_APPROX( nrg );
310        Qnrg >>= 1;             /* range: -6...15*/
311
312        psEncCtrl->Gains_Q16[ k ] = silk_LSHIFT_SAT32( tmp32, 16 - Qnrg );
313
314        if( psEnc->sCmn.warping_Q16 > 0 ) {
315            /* Adjust gain for warping */
316            gain_mult_Q16 = warped_gain( AR2_Q24, warping_Q16, psEnc->sCmn.shapingLPCOrder );
317            silk_assert( psEncCtrl->Gains_Q16[ k ] >= 0 );
318            if ( silk_SMULWW( silk_RSHIFT_ROUND( psEncCtrl->Gains_Q16[ k ], 1 ), gain_mult_Q16 ) >= ( silk_int32_MAX >> 1 ) ) {
319               psEncCtrl->Gains_Q16[ k ] = silk_int32_MAX;
320            } else {
321               psEncCtrl->Gains_Q16[ k ] = silk_SMULWW( psEncCtrl->Gains_Q16[ k ], gain_mult_Q16 );
322            }
323        }
324
325        /* Bandwidth expansion for synthesis filter shaping */
326        silk_bwexpander_32( AR2_Q24, psEnc->sCmn.shapingLPCOrder, BWExp2_Q16 );
327
328        /* Compute noise shaping filter coefficients */
329        silk_memcpy( AR1_Q24, AR2_Q24, psEnc->sCmn.shapingLPCOrder * sizeof( opus_int32 ) );
330
331        /* Bandwidth expansion for analysis filter shaping */
332        silk_assert( BWExp1_Q16 <= SILK_FIX_CONST( 1.0, 16 ) );
333        silk_bwexpander_32( AR1_Q24, psEnc->sCmn.shapingLPCOrder, BWExp1_Q16 );
334
335        /* Ratio of prediction gains, in energy domain */
336        pre_nrg_Q30 = silk_LPC_inverse_pred_gain_Q24( AR2_Q24, psEnc->sCmn.shapingLPCOrder );
337        nrg         = silk_LPC_inverse_pred_gain_Q24( AR1_Q24, psEnc->sCmn.shapingLPCOrder );
338
339        /*psEncCtrl->GainsPre[ k ] = 1.0f - 0.7f * ( 1.0f - pre_nrg / nrg ) = 0.3f + 0.7f * pre_nrg / nrg;*/
340        pre_nrg_Q30 = silk_LSHIFT32( silk_SMULWB( pre_nrg_Q30, SILK_FIX_CONST( 0.7, 15 ) ), 1 );
341        psEncCtrl->GainsPre_Q14[ k ] = ( opus_int ) SILK_FIX_CONST( 0.3, 14 ) + silk_DIV32_varQ( pre_nrg_Q30, nrg, 14 );
342
343        /* Convert to monic warped prediction coefficients and limit absolute values */
344        limit_warped_coefs( AR2_Q24, AR1_Q24, warping_Q16, SILK_FIX_CONST( 3.999, 24 ), psEnc->sCmn.shapingLPCOrder );
345
346        /* Convert from Q24 to Q13 and store in int16 */
347        for( i = 0; i < psEnc->sCmn.shapingLPCOrder; i++ ) {
348            psEncCtrl->AR1_Q13[ k * MAX_SHAPE_LPC_ORDER + i ] = (opus_int16)silk_SAT16( silk_RSHIFT_ROUND( AR1_Q24[ i ], 11 ) );
349            psEncCtrl->AR2_Q13[ k * MAX_SHAPE_LPC_ORDER + i ] = (opus_int16)silk_SAT16( silk_RSHIFT_ROUND( AR2_Q24[ i ], 11 ) );
350        }
351    }
352
353    /*****************/
354    /* Gain tweaking */
355    /*****************/
356    /* Increase gains during low speech activity and put lower limit on gains */
357    gain_mult_Q16 = silk_log2lin( -silk_SMLAWB( -SILK_FIX_CONST( 16.0, 7 ), SNR_adj_dB_Q7, SILK_FIX_CONST( 0.16, 16 ) ) );
358    gain_add_Q16  = silk_log2lin(  silk_SMLAWB(  SILK_FIX_CONST( 16.0, 7 ), SILK_FIX_CONST( MIN_QGAIN_DB, 7 ), SILK_FIX_CONST( 0.16, 16 ) ) );
359    silk_assert( gain_mult_Q16 > 0 );
360    for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
361        psEncCtrl->Gains_Q16[ k ] = silk_SMULWW( psEncCtrl->Gains_Q16[ k ], gain_mult_Q16 );
362        silk_assert( psEncCtrl->Gains_Q16[ k ] >= 0 );
363        psEncCtrl->Gains_Q16[ k ] = silk_ADD_POS_SAT32( psEncCtrl->Gains_Q16[ k ], gain_add_Q16 );
364    }
365
366    gain_mult_Q16 = SILK_FIX_CONST( 1.0, 16 ) + silk_RSHIFT_ROUND( silk_MLA( SILK_FIX_CONST( INPUT_TILT, 26 ),
367        psEncCtrl->coding_quality_Q14, SILK_FIX_CONST( HIGH_RATE_INPUT_TILT, 12 ) ), 10 );
368    for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
369        psEncCtrl->GainsPre_Q14[ k ] = silk_SMULWB( gain_mult_Q16, psEncCtrl->GainsPre_Q14[ k ] );
370    }
371
372    /************************************************/
373    /* Control low-frequency shaping and noise tilt */
374    /************************************************/
375    /* Less low frequency shaping for noisy inputs */
376    strength_Q16 = silk_MUL( SILK_FIX_CONST( LOW_FREQ_SHAPING, 4 ), silk_SMLAWB( SILK_FIX_CONST( 1.0, 12 ),
377        SILK_FIX_CONST( LOW_QUALITY_LOW_FREQ_SHAPING_DECR, 13 ), psEnc->sCmn.input_quality_bands_Q15[ 0 ] - SILK_FIX_CONST( 1.0, 15 ) ) );
378    strength_Q16 = silk_RSHIFT( silk_MUL( strength_Q16, psEnc->sCmn.speech_activity_Q8 ), 8 );
379    if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
380        /* Reduce low frequencies quantization noise for periodic signals, depending on pitch lag */
381        /*f = 400; freqz([1, -0.98 + 2e-4 * f], [1, -0.97 + 7e-4 * f], 2^12, Fs); axis([0, 1000, -10, 1])*/
382        opus_int fs_kHz_inv = silk_DIV32_16( SILK_FIX_CONST( 0.2, 14 ), psEnc->sCmn.fs_kHz );
383        for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
384            b_Q14 = fs_kHz_inv + silk_DIV32_16( SILK_FIX_CONST( 3.0, 14 ), psEncCtrl->pitchL[ k ] );
385            /* Pack two coefficients in one int32 */
386            psEncCtrl->LF_shp_Q14[ k ]  = silk_LSHIFT( SILK_FIX_CONST( 1.0, 14 ) - b_Q14 - silk_SMULWB( strength_Q16, b_Q14 ), 16 );
387            psEncCtrl->LF_shp_Q14[ k ] |= (opus_uint16)( b_Q14 - SILK_FIX_CONST( 1.0, 14 ) );
388        }
389        silk_assert( SILK_FIX_CONST( HARM_HP_NOISE_COEF, 24 ) < SILK_FIX_CONST( 0.5, 24 ) ); /* Guarantees that second argument to SMULWB() is within range of an opus_int16*/
390        Tilt_Q16 = - SILK_FIX_CONST( HP_NOISE_COEF, 16 ) -
391            silk_SMULWB( SILK_FIX_CONST( 1.0, 16 ) - SILK_FIX_CONST( HP_NOISE_COEF, 16 ),
392                silk_SMULWB( SILK_FIX_CONST( HARM_HP_NOISE_COEF, 24 ), psEnc->sCmn.speech_activity_Q8 ) );
393    } else {
394        b_Q14 = silk_DIV32_16( 21299, psEnc->sCmn.fs_kHz ); /* 1.3_Q0 = 21299_Q14*/
395        /* Pack two coefficients in one int32 */
396        psEncCtrl->LF_shp_Q14[ 0 ]  = silk_LSHIFT( SILK_FIX_CONST( 1.0, 14 ) - b_Q14 -
397            silk_SMULWB( strength_Q16, silk_SMULWB( SILK_FIX_CONST( 0.6, 16 ), b_Q14 ) ), 16 );
398        psEncCtrl->LF_shp_Q14[ 0 ] |= (opus_uint16)( b_Q14 - SILK_FIX_CONST( 1.0, 14 ) );
399        for( k = 1; k < psEnc->sCmn.nb_subfr; k++ ) {
400            psEncCtrl->LF_shp_Q14[ k ] = psEncCtrl->LF_shp_Q14[ 0 ];
401        }
402        Tilt_Q16 = -SILK_FIX_CONST( HP_NOISE_COEF, 16 );
403    }
404
405    /****************************/
406    /* HARMONIC SHAPING CONTROL */
407    /****************************/
408    /* Control boosting of harmonic frequencies */
409    HarmBoost_Q16 = silk_SMULWB( silk_SMULWB( SILK_FIX_CONST( 1.0, 17 ) - silk_LSHIFT( psEncCtrl->coding_quality_Q14, 3 ),
410        psEnc->LTPCorr_Q15 ), SILK_FIX_CONST( LOW_RATE_HARMONIC_BOOST, 16 ) );
411
412    /* More harmonic boost for noisy input signals */
413    HarmBoost_Q16 = silk_SMLAWB( HarmBoost_Q16,
414        SILK_FIX_CONST( 1.0, 16 ) - silk_LSHIFT( psEncCtrl->input_quality_Q14, 2 ), SILK_FIX_CONST( LOW_INPUT_QUALITY_HARMONIC_BOOST, 16 ) );
415
416    if( USE_HARM_SHAPING && psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
417        /* More harmonic noise shaping for high bitrates or noisy input */
418        HarmShapeGain_Q16 = silk_SMLAWB( SILK_FIX_CONST( HARMONIC_SHAPING, 16 ),
419                SILK_FIX_CONST( 1.0, 16 ) - silk_SMULWB( SILK_FIX_CONST( 1.0, 18 ) - silk_LSHIFT( psEncCtrl->coding_quality_Q14, 4 ),
420                psEncCtrl->input_quality_Q14 ), SILK_FIX_CONST( HIGH_RATE_OR_LOW_QUALITY_HARMONIC_SHAPING, 16 ) );
421
422        /* Less harmonic noise shaping for less periodic signals */
423        HarmShapeGain_Q16 = silk_SMULWB( silk_LSHIFT( HarmShapeGain_Q16, 1 ),
424            silk_SQRT_APPROX( silk_LSHIFT( psEnc->LTPCorr_Q15, 15 ) ) );
425    } else {
426        HarmShapeGain_Q16 = 0;
427    }
428
429    /*************************/
430    /* Smooth over subframes */
431    /*************************/
432    for( k = 0; k < MAX_NB_SUBFR; k++ ) {
433        psShapeSt->HarmBoost_smth_Q16 =
434            silk_SMLAWB( psShapeSt->HarmBoost_smth_Q16,     HarmBoost_Q16     - psShapeSt->HarmBoost_smth_Q16,     SILK_FIX_CONST( SUBFR_SMTH_COEF, 16 ) );
435        psShapeSt->HarmShapeGain_smth_Q16 =
436            silk_SMLAWB( psShapeSt->HarmShapeGain_smth_Q16, HarmShapeGain_Q16 - psShapeSt->HarmShapeGain_smth_Q16, SILK_FIX_CONST( SUBFR_SMTH_COEF, 16 ) );
437        psShapeSt->Tilt_smth_Q16 =
438            silk_SMLAWB( psShapeSt->Tilt_smth_Q16,          Tilt_Q16          - psShapeSt->Tilt_smth_Q16,          SILK_FIX_CONST( SUBFR_SMTH_COEF, 16 ) );
439
440        psEncCtrl->HarmBoost_Q14[ k ]     = ( opus_int )silk_RSHIFT_ROUND( psShapeSt->HarmBoost_smth_Q16,     2 );
441        psEncCtrl->HarmShapeGain_Q14[ k ] = ( opus_int )silk_RSHIFT_ROUND( psShapeSt->HarmShapeGain_smth_Q16, 2 );
442        psEncCtrl->Tilt_Q14[ k ]          = ( opus_int )silk_RSHIFT_ROUND( psShapeSt->Tilt_smth_Q16,          2 );
443    }
444    RESTORE_STACK;
445}
446