Fixes an overflow in limit_warped_coefs()
[opus.git] / silk / fixed / noise_shape_analysis_FIX.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, are permitted provided that the following conditions
5 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 Internet Society, IETF or IETF Trust, nor the
12 names of specific contributors, may be used to endorse or promote
13 products derived from this software without specific prior written
14 permission.
15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
16 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
19 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
20 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
21 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
22 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
23 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
24 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
25 POSSIBILITY 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.                                    */
40 static 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 */
59 static OPUS_INLINE void limit_warped_coefs(
60     opus_int32           *coefs_Q24,
61     opus_int             lambda_Q16,
62     opus_int32           limit_Q24,
63     opus_int             order
64 ) {
65     opus_int   i, iter, ind = 0;
66     opus_int32 tmp, maxabs_Q24, chirp_Q16, gain_Q16;
67     opus_int32 nom_Q16, den_Q24;
68     opus_int32 limit_Q20, maxabs_Q20;
69
70     /* Convert to monic coefficients */
71     lambda_Q16 = -lambda_Q16;
72     for( i = order - 1; i > 0; i-- ) {
73         coefs_Q24[ i - 1 ] = silk_SMLAWB( coefs_Q24[ i - 1 ], coefs_Q24[ i ], lambda_Q16 );
74     }
75     lambda_Q16 = -lambda_Q16;
76     nom_Q16  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 16 ), -(opus_int32)lambda_Q16, lambda_Q16 );
77     den_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), coefs_Q24[ 0 ], lambda_Q16 );
78     gain_Q16 = silk_DIV32_varQ( nom_Q16, den_Q24, 24 );
79     for( i = 0; i < order; i++ ) {
80         coefs_Q24[ i ] = silk_SMULWW( gain_Q16, coefs_Q24[ i ] );
81     }
82     limit_Q20 = silk_RSHIFT(limit_Q24, 4);
83     for( iter = 0; iter < 10; iter++ ) {
84         /* Find maximum absolute value */
85         maxabs_Q24 = -1;
86         for( i = 0; i < order; i++ ) {
87             tmp = silk_abs_int32( coefs_Q24[ i ] );
88             if( tmp > maxabs_Q24 ) {
89                 maxabs_Q24 = tmp;
90                 ind = i;
91             }
92         }
93         /* Use Q20 to avoid any overflow when multiplying by (ind + 1) later. */
94         maxabs_Q20 = silk_RSHIFT(maxabs_Q24, 4);
95         if( maxabs_Q20 <= limit_Q20 ) {
96             /* Coefficients are within range - done */
97             return;
98         }
99
100         /* Convert back to true warped coefficients */
101         for( i = 1; i < order; i++ ) {
102             coefs_Q24[ i - 1 ] = silk_SMLAWB( coefs_Q24[ i - 1 ], coefs_Q24[ i ], lambda_Q16 );
103         }
104         gain_Q16 = silk_INVERSE32_varQ( gain_Q16, 32 );
105         for( i = 0; i < order; i++ ) {
106             coefs_Q24[ i ] = silk_SMULWW( gain_Q16, coefs_Q24[ i ] );
107         }
108
109         /* Apply bandwidth expansion */
110         chirp_Q16 = SILK_FIX_CONST( 0.99, 16 ) - silk_DIV32_varQ(
111             silk_SMULWB( maxabs_Q20 - limit_Q20, silk_SMLABB( SILK_FIX_CONST( 0.8, 10 ), SILK_FIX_CONST( 0.1, 10 ), iter ) ),
112             silk_MUL( maxabs_Q20, ind + 1 ), 22 );
113         silk_bwexpander_32( coefs_Q24, order, chirp_Q16 );
114
115         /* Convert to monic warped coefficients */
116         lambda_Q16 = -lambda_Q16;
117         for( i = order - 1; i > 0; i-- ) {
118             coefs_Q24[ i - 1 ] = silk_SMLAWB( coefs_Q24[ i - 1 ], coefs_Q24[ i ], lambda_Q16 );
119         }
120         lambda_Q16 = -lambda_Q16;
121         nom_Q16  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 16 ), -(opus_int32)lambda_Q16,        lambda_Q16 );
122         den_Q24  = silk_SMLAWB( SILK_FIX_CONST( 1.0, 24 ), coefs_Q24[ 0 ], lambda_Q16 );
123         gain_Q16 = silk_DIV32_varQ( nom_Q16, den_Q24, 24 );
124         for( i = 0; i < order; i++ ) {
125             coefs_Q24[ i ] = silk_SMULWW( gain_Q16, coefs_Q24[ i ] );
126         }
127     }
128     silk_assert( 0 );
129 }
130
131 /* Disable MIPS version until it's updated. */
132 #if 0 && defined(MIPSr1_ASM)
133 #include "mips/noise_shape_analysis_FIX_mipsr1.h"
134 #endif
135
136 /**************************************************************/
137 /* Compute noise shaping coefficients and initial gain values */
138 /**************************************************************/
139 #ifndef OVERRIDE_silk_noise_shape_analysis_FIX
140 void silk_noise_shape_analysis_FIX(
141     silk_encoder_state_FIX          *psEnc,                                 /* I/O  Encoder state FIX                                                           */
142     silk_encoder_control_FIX        *psEncCtrl,                             /* I/O  Encoder control FIX                                                         */
143     const opus_int16                *pitch_res,                             /* I    LPC residual from pitch analysis                                            */
144     const opus_int16                *x,                                     /* I    Input signal [ frame_length + la_shape ]                                    */
145     int                              arch                                   /* I    Run-time architecture                                                       */
146 )
147 {
148     silk_shape_state_FIX *psShapeSt = &psEnc->sShape;
149     opus_int     k, i, nSamples, nSegs, Qnrg, b_Q14, warping_Q16, scale = 0;
150     opus_int32   SNR_adj_dB_Q7, HarmShapeGain_Q16, Tilt_Q16, tmp32;
151     opus_int32   nrg, log_energy_Q7, log_energy_prev_Q7, energy_variation_Q7;
152     opus_int32   BWExp_Q16, gain_mult_Q16, gain_add_Q16, strength_Q16, b_Q8;
153     opus_int32   auto_corr[     MAX_SHAPE_LPC_ORDER + 1 ];
154     opus_int32   refl_coef_Q16[ MAX_SHAPE_LPC_ORDER ];
155     opus_int32   AR_Q24[       MAX_SHAPE_LPC_ORDER ];
156     VARDECL( opus_int16, x_windowed );
157     const opus_int16 *x_ptr, *pitch_res_ptr;
158     SAVE_STACK;
159
160     /* Point to start of first LPC analysis block */
161     x_ptr = x - psEnc->sCmn.la_shape;
162
163     /****************/
164     /* GAIN CONTROL */
165     /****************/
166     SNR_adj_dB_Q7 = psEnc->sCmn.SNR_dB_Q7;
167
168     /* Input quality is the average of the quality in the lowest two VAD bands */
169     psEncCtrl->input_quality_Q14 = ( opus_int )silk_RSHIFT( (opus_int32)psEnc->sCmn.input_quality_bands_Q15[ 0 ]
170         + psEnc->sCmn.input_quality_bands_Q15[ 1 ], 2 );
171
172     /* Coding quality level, between 0.0_Q0 and 1.0_Q0, but in Q14 */
173     psEncCtrl->coding_quality_Q14 = silk_RSHIFT( silk_sigm_Q15( silk_RSHIFT_ROUND( SNR_adj_dB_Q7 -
174         SILK_FIX_CONST( 20.0, 7 ), 4 ) ), 1 );
175
176     /* Reduce coding SNR during low speech activity */
177     if( psEnc->sCmn.useCBR == 0 ) {
178         b_Q8 = SILK_FIX_CONST( 1.0, 8 ) - psEnc->sCmn.speech_activity_Q8;
179         b_Q8 = silk_SMULWB( silk_LSHIFT( b_Q8, 8 ), b_Q8 );
180         SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7,
181             silk_SMULBB( SILK_FIX_CONST( -BG_SNR_DECR_dB, 7 ) >> ( 4 + 1 ), b_Q8 ),                                       /* Q11*/
182             silk_SMULWB( SILK_FIX_CONST( 1.0, 14 ) + psEncCtrl->input_quality_Q14, psEncCtrl->coding_quality_Q14 ) );     /* Q12*/
183     }
184
185     if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
186         /* Reduce gains for periodic signals */
187         SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7, SILK_FIX_CONST( HARM_SNR_INCR_dB, 8 ), psEnc->LTPCorr_Q15 );
188     } else {
189         /* For unvoiced signals and low-quality input, adjust the quality slower than SNR_dB setting */
190         SNR_adj_dB_Q7 = silk_SMLAWB( SNR_adj_dB_Q7,
191             silk_SMLAWB( SILK_FIX_CONST( 6.0, 9 ), -SILK_FIX_CONST( 0.4, 18 ), psEnc->sCmn.SNR_dB_Q7 ),
192             SILK_FIX_CONST( 1.0, 14 ) - psEncCtrl->input_quality_Q14 );
193     }
194
195     /*************************/
196     /* SPARSENESS PROCESSING */
197     /*************************/
198     /* Set quantizer offset */
199     if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
200         /* Initially set to 0; may be overruled in process_gains(..) */
201         psEnc->sCmn.indices.quantOffsetType = 0;
202     } else {
203         /* Sparseness measure, based on relative fluctuations of energy per 2 milliseconds */
204         nSamples = silk_LSHIFT( psEnc->sCmn.fs_kHz, 1 );
205         energy_variation_Q7 = 0;
206         log_energy_prev_Q7  = 0;
207         pitch_res_ptr = pitch_res;
208         nSegs = silk_SMULBB( SUB_FRAME_LENGTH_MS, psEnc->sCmn.nb_subfr ) / 2;
209         for( k = 0; k < nSegs; k++ ) {
210             silk_sum_sqr_shift( &nrg, &scale, pitch_res_ptr, nSamples );
211             nrg += silk_RSHIFT( nSamples, scale );           /* Q(-scale)*/
212
213             log_energy_Q7 = silk_lin2log( nrg );
214             if( k > 0 ) {
215                 energy_variation_Q7 += silk_abs( log_energy_Q7 - log_energy_prev_Q7 );
216             }
217             log_energy_prev_Q7 = log_energy_Q7;
218             pitch_res_ptr += nSamples;
219         }
220
221         /* Set quantization offset depending on sparseness measure */
222         if( energy_variation_Q7 > SILK_FIX_CONST( ENERGY_VARIATION_THRESHOLD_QNT_OFFSET, 7 ) * (nSegs-1) ) {
223             psEnc->sCmn.indices.quantOffsetType = 0;
224         } else {
225             psEnc->sCmn.indices.quantOffsetType = 1;
226         }
227     }
228
229     /*******************************/
230     /* Control bandwidth expansion */
231     /*******************************/
232     /* More BWE for signals with high prediction gain */
233     strength_Q16 = silk_SMULWB( psEncCtrl->predGain_Q16, SILK_FIX_CONST( FIND_PITCH_WHITE_NOISE_FRACTION, 16 ) );
234     BWExp_Q16 = silk_DIV32_varQ( SILK_FIX_CONST( BANDWIDTH_EXPANSION, 16 ),
235         silk_SMLAWW( SILK_FIX_CONST( 1.0, 16 ), strength_Q16, strength_Q16 ), 16 );
236
237     if( psEnc->sCmn.warping_Q16 > 0 ) {
238         /* Slightly more warping in analysis will move quantization noise up in frequency, where it's better masked */
239         warping_Q16 = silk_SMLAWB( psEnc->sCmn.warping_Q16, (opus_int32)psEncCtrl->coding_quality_Q14, SILK_FIX_CONST( 0.01, 18 ) );
240     } else {
241         warping_Q16 = 0;
242     }
243
244     /********************************************/
245     /* Compute noise shaping AR coefs and gains */
246     /********************************************/
247     ALLOC( x_windowed, psEnc->sCmn.shapeWinLength, opus_int16 );
248     for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
249         /* Apply window: sine slope followed by flat part followed by cosine slope */
250         opus_int shift, slope_part, flat_part;
251         flat_part = psEnc->sCmn.fs_kHz * 3;
252         slope_part = silk_RSHIFT( psEnc->sCmn.shapeWinLength - flat_part, 1 );
253
254         silk_apply_sine_window( x_windowed, x_ptr, 1, slope_part );
255         shift = slope_part;
256         silk_memcpy( x_windowed + shift, x_ptr + shift, flat_part * sizeof(opus_int16) );
257         shift += flat_part;
258         silk_apply_sine_window( x_windowed + shift, x_ptr + shift, 2, slope_part );
259
260         /* Update pointer: next LPC analysis block */
261         x_ptr += psEnc->sCmn.subfr_length;
262
263         if( psEnc->sCmn.warping_Q16 > 0 ) {
264             /* Calculate warped auto correlation */
265             silk_warped_autocorrelation_FIX( auto_corr, &scale, x_windowed, warping_Q16, psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder );
266         } else {
267             /* Calculate regular auto correlation */
268             silk_autocorr( auto_corr, &scale, x_windowed, psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder + 1, arch );
269         }
270
271         /* Add white noise, as a fraction of energy */
272         auto_corr[0] = silk_ADD32( auto_corr[0], silk_max_32( silk_SMULWB( silk_RSHIFT( auto_corr[ 0 ], 4 ),
273             SILK_FIX_CONST( SHAPE_WHITE_NOISE_FRACTION, 20 ) ), 1 ) );
274
275         /* Calculate the reflection coefficients using schur */
276         nrg = silk_schur64( refl_coef_Q16, auto_corr, psEnc->sCmn.shapingLPCOrder );
277         silk_assert( nrg >= 0 );
278
279         /* Convert reflection coefficients to prediction coefficients */
280         silk_k2a_Q16( AR_Q24, refl_coef_Q16, psEnc->sCmn.shapingLPCOrder );
281
282         Qnrg = -scale;          /* range: -12...30*/
283         silk_assert( Qnrg >= -12 );
284         silk_assert( Qnrg <=  30 );
285
286         /* Make sure that Qnrg is an even number */
287         if( Qnrg & 1 ) {
288             Qnrg -= 1;
289             nrg >>= 1;
290         }
291
292         tmp32 = silk_SQRT_APPROX( nrg );
293         Qnrg >>= 1;             /* range: -6...15*/
294
295         psEncCtrl->Gains_Q16[ k ] = silk_LSHIFT_SAT32( tmp32, 16 - Qnrg );
296
297         if( psEnc->sCmn.warping_Q16 > 0 ) {
298             /* Adjust gain for warping */
299             gain_mult_Q16 = warped_gain( AR_Q24, warping_Q16, psEnc->sCmn.shapingLPCOrder );
300             silk_assert( psEncCtrl->Gains_Q16[ k ] > 0 );
301             if( psEncCtrl->Gains_Q16[ k ] < SILK_FIX_CONST( 0.25, 16 ) ) {
302                 psEncCtrl->Gains_Q16[ k ] = silk_SMULWW( psEncCtrl->Gains_Q16[ k ], gain_mult_Q16 ); 
303             } else {
304                 psEncCtrl->Gains_Q16[ k ] = silk_SMULWW( silk_RSHIFT_ROUND( psEncCtrl->Gains_Q16[ k ], 1 ), gain_mult_Q16 );
305                 if ( psEncCtrl->Gains_Q16[ k ] >= ( silk_int32_MAX >> 1 ) ) {
306                     psEncCtrl->Gains_Q16[ k ] = silk_int32_MAX;
307                 } else {
308                     psEncCtrl->Gains_Q16[ k ] = silk_LSHIFT32( psEncCtrl->Gains_Q16[ k ], 1 );
309                 }
310             }
311             silk_assert( psEncCtrl->Gains_Q16[ k ] > 0 );
312         }
313
314         /* Bandwidth expansion */
315         silk_bwexpander_32( AR_Q24, psEnc->sCmn.shapingLPCOrder, BWExp_Q16 );
316
317         if( psEnc->sCmn.warping_Q16 > 0 ) {
318             /* Convert to monic warped prediction coefficients and limit absolute values */
319             limit_warped_coefs( AR_Q24, warping_Q16, SILK_FIX_CONST( 3.999, 24 ), psEnc->sCmn.shapingLPCOrder );
320
321             /* Convert from Q24 to Q13 and store in int16 */
322             for( i = 0; i < psEnc->sCmn.shapingLPCOrder; i++ ) {
323                 psEncCtrl->AR_Q13[ k * MAX_SHAPE_LPC_ORDER + i ] = (opus_int16)silk_SAT16( silk_RSHIFT_ROUND( AR_Q24[ i ], 11 ) );
324             }
325         } else {
326             silk_LPC_fit( &psEncCtrl->AR_Q13[ k * MAX_SHAPE_LPC_ORDER ], AR_Q24, 13, 24, psEnc->sCmn.shapingLPCOrder );
327         }
328     }
329
330     /*****************/
331     /* Gain tweaking */
332     /*****************/
333     /* Increase gains during low speech activity and put lower limit on gains */
334     gain_mult_Q16 = silk_log2lin( -silk_SMLAWB( -SILK_FIX_CONST( 16.0, 7 ), SNR_adj_dB_Q7, SILK_FIX_CONST( 0.16, 16 ) ) );
335     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 ) ) );
336     silk_assert( gain_mult_Q16 > 0 );
337     for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
338         psEncCtrl->Gains_Q16[ k ] = silk_SMULWW( psEncCtrl->Gains_Q16[ k ], gain_mult_Q16 );
339         silk_assert( psEncCtrl->Gains_Q16[ k ] >= 0 );
340         psEncCtrl->Gains_Q16[ k ] = silk_ADD_POS_SAT32( psEncCtrl->Gains_Q16[ k ], gain_add_Q16 );
341     }
342
343
344     /************************************************/
345     /* Control low-frequency shaping and noise tilt */
346     /************************************************/
347     /* Less low frequency shaping for noisy inputs */
348     strength_Q16 = silk_MUL( SILK_FIX_CONST( LOW_FREQ_SHAPING, 4 ), silk_SMLAWB( SILK_FIX_CONST( 1.0, 12 ),
349         SILK_FIX_CONST( LOW_QUALITY_LOW_FREQ_SHAPING_DECR, 13 ), psEnc->sCmn.input_quality_bands_Q15[ 0 ] - SILK_FIX_CONST( 1.0, 15 ) ) );
350     strength_Q16 = silk_RSHIFT( silk_MUL( strength_Q16, psEnc->sCmn.speech_activity_Q8 ), 8 );
351     if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
352         /* Reduce low frequencies quantization noise for periodic signals, depending on pitch lag */
353         /*f = 400; freqz([1, -0.98 + 2e-4 * f], [1, -0.97 + 7e-4 * f], 2^12, Fs); axis([0, 1000, -10, 1])*/
354         opus_int fs_kHz_inv = silk_DIV32_16( SILK_FIX_CONST( 0.2, 14 ), psEnc->sCmn.fs_kHz );
355         for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
356             b_Q14 = fs_kHz_inv + silk_DIV32_16( SILK_FIX_CONST( 3.0, 14 ), psEncCtrl->pitchL[ k ] );
357             /* Pack two coefficients in one int32 */
358             psEncCtrl->LF_shp_Q14[ k ]  = silk_LSHIFT( SILK_FIX_CONST( 1.0, 14 ) - b_Q14 - silk_SMULWB( strength_Q16, b_Q14 ), 16 );
359             psEncCtrl->LF_shp_Q14[ k ] |= (opus_uint16)( b_Q14 - SILK_FIX_CONST( 1.0, 14 ) );
360         }
361         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*/
362         Tilt_Q16 = - SILK_FIX_CONST( HP_NOISE_COEF, 16 ) -
363             silk_SMULWB( SILK_FIX_CONST( 1.0, 16 ) - SILK_FIX_CONST( HP_NOISE_COEF, 16 ),
364                 silk_SMULWB( SILK_FIX_CONST( HARM_HP_NOISE_COEF, 24 ), psEnc->sCmn.speech_activity_Q8 ) );
365     } else {
366         b_Q14 = silk_DIV32_16( 21299, psEnc->sCmn.fs_kHz ); /* 1.3_Q0 = 21299_Q14*/
367         /* Pack two coefficients in one int32 */
368         psEncCtrl->LF_shp_Q14[ 0 ]  = silk_LSHIFT( SILK_FIX_CONST( 1.0, 14 ) - b_Q14 -
369             silk_SMULWB( strength_Q16, silk_SMULWB( SILK_FIX_CONST( 0.6, 16 ), b_Q14 ) ), 16 );
370         psEncCtrl->LF_shp_Q14[ 0 ] |= (opus_uint16)( b_Q14 - SILK_FIX_CONST( 1.0, 14 ) );
371         for( k = 1; k < psEnc->sCmn.nb_subfr; k++ ) {
372             psEncCtrl->LF_shp_Q14[ k ] = psEncCtrl->LF_shp_Q14[ 0 ];
373         }
374         Tilt_Q16 = -SILK_FIX_CONST( HP_NOISE_COEF, 16 );
375     }
376
377     /****************************/
378     /* HARMONIC SHAPING CONTROL */
379     /****************************/
380     if( USE_HARM_SHAPING && psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
381         /* More harmonic noise shaping for high bitrates or noisy input */
382         HarmShapeGain_Q16 = silk_SMLAWB( SILK_FIX_CONST( HARMONIC_SHAPING, 16 ),
383                 SILK_FIX_CONST( 1.0, 16 ) - silk_SMULWB( SILK_FIX_CONST( 1.0, 18 ) - silk_LSHIFT( psEncCtrl->coding_quality_Q14, 4 ),
384                 psEncCtrl->input_quality_Q14 ), SILK_FIX_CONST( HIGH_RATE_OR_LOW_QUALITY_HARMONIC_SHAPING, 16 ) );
385
386         /* Less harmonic noise shaping for less periodic signals */
387         HarmShapeGain_Q16 = silk_SMULWB( silk_LSHIFT( HarmShapeGain_Q16, 1 ),
388             silk_SQRT_APPROX( silk_LSHIFT( psEnc->LTPCorr_Q15, 15 ) ) );
389     } else {
390         HarmShapeGain_Q16 = 0;
391     }
392
393     /*************************/
394     /* Smooth over subframes */
395     /*************************/
396     for( k = 0; k < MAX_NB_SUBFR; k++ ) {
397         psShapeSt->HarmShapeGain_smth_Q16 =
398             silk_SMLAWB( psShapeSt->HarmShapeGain_smth_Q16, HarmShapeGain_Q16 - psShapeSt->HarmShapeGain_smth_Q16, SILK_FIX_CONST( SUBFR_SMTH_COEF, 16 ) );
399         psShapeSt->Tilt_smth_Q16 =
400             silk_SMLAWB( psShapeSt->Tilt_smth_Q16,          Tilt_Q16          - psShapeSt->Tilt_smth_Q16,          SILK_FIX_CONST( SUBFR_SMTH_COEF, 16 ) );
401
402         psEncCtrl->HarmShapeGain_Q14[ k ] = ( opus_int )silk_RSHIFT_ROUND( psShapeSt->HarmShapeGain_smth_Q16, 2 );
403         psEncCtrl->Tilt_Q14[ k ]          = ( opus_int )silk_RSHIFT_ROUND( psShapeSt->Tilt_smth_Q16,          2 );
404     }
405     RESTORE_STACK;
406 }
407 #endif /* OVERRIDE_silk_noise_shape_analysis_FIX */