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