1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/silk/float/noise_shape_analysis_FLP.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,365 @@
1.4 +/***********************************************************************
1.5 +Copyright (c) 2006-2011, Skype Limited. All rights reserved.
1.6 +Redistribution and use in source and binary forms, with or without
1.7 +modification, are permitted provided that the following conditions
1.8 +are met:
1.9 +- Redistributions of source code must retain the above copyright notice,
1.10 +this list of conditions and the following disclaimer.
1.11 +- Redistributions in binary form must reproduce the above copyright
1.12 +notice, this list of conditions and the following disclaimer in the
1.13 +documentation and/or other materials provided with the distribution.
1.14 +- Neither the name of Internet Society, IETF or IETF Trust, nor the
1.15 +names of specific contributors, may be used to endorse or promote
1.16 +products derived from this software without specific prior written
1.17 +permission.
1.18 +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
1.19 +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1.20 +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
1.21 +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
1.22 +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
1.23 +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
1.24 +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
1.25 +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
1.26 +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
1.27 +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
1.28 +POSSIBILITY OF SUCH DAMAGE.
1.29 +***********************************************************************/
1.30 +
1.31 +#ifdef HAVE_CONFIG_H
1.32 +#include "config.h"
1.33 +#endif
1.34 +
1.35 +#include "main_FLP.h"
1.36 +#include "tuning_parameters.h"
1.37 +
1.38 +/* Compute gain to make warped filter coefficients have a zero mean log frequency response on a */
1.39 +/* non-warped frequency scale. (So that it can be implemented with a minimum-phase monic filter.) */
1.40 +/* Note: A monic filter is one with the first coefficient equal to 1.0. In Silk we omit the first */
1.41 +/* coefficient in an array of coefficients, for monic filters. */
1.42 +static OPUS_INLINE silk_float warped_gain(
1.43 + const silk_float *coefs,
1.44 + silk_float lambda,
1.45 + opus_int order
1.46 +) {
1.47 + opus_int i;
1.48 + silk_float gain;
1.49 +
1.50 + lambda = -lambda;
1.51 + gain = coefs[ order - 1 ];
1.52 + for( i = order - 2; i >= 0; i-- ) {
1.53 + gain = lambda * gain + coefs[ i ];
1.54 + }
1.55 + return (silk_float)( 1.0f / ( 1.0f - lambda * gain ) );
1.56 +}
1.57 +
1.58 +/* Convert warped filter coefficients to monic pseudo-warped coefficients and limit maximum */
1.59 +/* amplitude of monic warped coefficients by using bandwidth expansion on the true coefficients */
1.60 +static OPUS_INLINE void warped_true2monic_coefs(
1.61 + silk_float *coefs_syn,
1.62 + silk_float *coefs_ana,
1.63 + silk_float lambda,
1.64 + silk_float limit,
1.65 + opus_int order
1.66 +) {
1.67 + opus_int i, iter, ind = 0;
1.68 + silk_float tmp, maxabs, chirp, gain_syn, gain_ana;
1.69 +
1.70 + /* Convert to monic coefficients */
1.71 + for( i = order - 1; i > 0; i-- ) {
1.72 + coefs_syn[ i - 1 ] -= lambda * coefs_syn[ i ];
1.73 + coefs_ana[ i - 1 ] -= lambda * coefs_ana[ i ];
1.74 + }
1.75 + gain_syn = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_syn[ 0 ] );
1.76 + gain_ana = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_ana[ 0 ] );
1.77 + for( i = 0; i < order; i++ ) {
1.78 + coefs_syn[ i ] *= gain_syn;
1.79 + coefs_ana[ i ] *= gain_ana;
1.80 + }
1.81 +
1.82 + /* Limit */
1.83 + for( iter = 0; iter < 10; iter++ ) {
1.84 + /* Find maximum absolute value */
1.85 + maxabs = -1.0f;
1.86 + for( i = 0; i < order; i++ ) {
1.87 + tmp = silk_max( silk_abs_float( coefs_syn[ i ] ), silk_abs_float( coefs_ana[ i ] ) );
1.88 + if( tmp > maxabs ) {
1.89 + maxabs = tmp;
1.90 + ind = i;
1.91 + }
1.92 + }
1.93 + if( maxabs <= limit ) {
1.94 + /* Coefficients are within range - done */
1.95 + return;
1.96 + }
1.97 +
1.98 + /* Convert back to true warped coefficients */
1.99 + for( i = 1; i < order; i++ ) {
1.100 + coefs_syn[ i - 1 ] += lambda * coefs_syn[ i ];
1.101 + coefs_ana[ i - 1 ] += lambda * coefs_ana[ i ];
1.102 + }
1.103 + gain_syn = 1.0f / gain_syn;
1.104 + gain_ana = 1.0f / gain_ana;
1.105 + for( i = 0; i < order; i++ ) {
1.106 + coefs_syn[ i ] *= gain_syn;
1.107 + coefs_ana[ i ] *= gain_ana;
1.108 + }
1.109 +
1.110 + /* Apply bandwidth expansion */
1.111 + chirp = 0.99f - ( 0.8f + 0.1f * iter ) * ( maxabs - limit ) / ( maxabs * ( ind + 1 ) );
1.112 + silk_bwexpander_FLP( coefs_syn, order, chirp );
1.113 + silk_bwexpander_FLP( coefs_ana, order, chirp );
1.114 +
1.115 + /* Convert to monic warped coefficients */
1.116 + for( i = order - 1; i > 0; i-- ) {
1.117 + coefs_syn[ i - 1 ] -= lambda * coefs_syn[ i ];
1.118 + coefs_ana[ i - 1 ] -= lambda * coefs_ana[ i ];
1.119 + }
1.120 + gain_syn = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_syn[ 0 ] );
1.121 + gain_ana = ( 1.0f - lambda * lambda ) / ( 1.0f + lambda * coefs_ana[ 0 ] );
1.122 + for( i = 0; i < order; i++ ) {
1.123 + coefs_syn[ i ] *= gain_syn;
1.124 + coefs_ana[ i ] *= gain_ana;
1.125 + }
1.126 + }
1.127 + silk_assert( 0 );
1.128 +}
1.129 +
1.130 +/* Compute noise shaping coefficients and initial gain values */
1.131 +void silk_noise_shape_analysis_FLP(
1.132 + silk_encoder_state_FLP *psEnc, /* I/O Encoder state FLP */
1.133 + silk_encoder_control_FLP *psEncCtrl, /* I/O Encoder control FLP */
1.134 + const silk_float *pitch_res, /* I LPC residual from pitch analysis */
1.135 + const silk_float *x /* I Input signal [frame_length + la_shape] */
1.136 +)
1.137 +{
1.138 + silk_shape_state_FLP *psShapeSt = &psEnc->sShape;
1.139 + opus_int k, nSamples;
1.140 + silk_float SNR_adj_dB, HarmBoost, HarmShapeGain, Tilt;
1.141 + silk_float nrg, pre_nrg, log_energy, log_energy_prev, energy_variation;
1.142 + silk_float delta, BWExp1, BWExp2, gain_mult, gain_add, strength, b, warping;
1.143 + silk_float x_windowed[ SHAPE_LPC_WIN_MAX ];
1.144 + silk_float auto_corr[ MAX_SHAPE_LPC_ORDER + 1 ];
1.145 + const silk_float *x_ptr, *pitch_res_ptr;
1.146 +
1.147 + /* Point to start of first LPC analysis block */
1.148 + x_ptr = x - psEnc->sCmn.la_shape;
1.149 +
1.150 + /****************/
1.151 + /* GAIN CONTROL */
1.152 + /****************/
1.153 + SNR_adj_dB = psEnc->sCmn.SNR_dB_Q7 * ( 1 / 128.0f );
1.154 +
1.155 + /* Input quality is the average of the quality in the lowest two VAD bands */
1.156 + psEncCtrl->input_quality = 0.5f * ( psEnc->sCmn.input_quality_bands_Q15[ 0 ] + psEnc->sCmn.input_quality_bands_Q15[ 1 ] ) * ( 1.0f / 32768.0f );
1.157 +
1.158 + /* Coding quality level, between 0.0 and 1.0 */
1.159 + psEncCtrl->coding_quality = silk_sigmoid( 0.25f * ( SNR_adj_dB - 20.0f ) );
1.160 +
1.161 + if( psEnc->sCmn.useCBR == 0 ) {
1.162 + /* Reduce coding SNR during low speech activity */
1.163 + b = 1.0f - psEnc->sCmn.speech_activity_Q8 * ( 1.0f / 256.0f );
1.164 + SNR_adj_dB -= BG_SNR_DECR_dB * psEncCtrl->coding_quality * ( 0.5f + 0.5f * psEncCtrl->input_quality ) * b * b;
1.165 + }
1.166 +
1.167 + if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
1.168 + /* Reduce gains for periodic signals */
1.169 + SNR_adj_dB += HARM_SNR_INCR_dB * psEnc->LTPCorr;
1.170 + } else {
1.171 + /* For unvoiced signals and low-quality input, adjust the quality slower than SNR_dB setting */
1.172 + SNR_adj_dB += ( -0.4f * psEnc->sCmn.SNR_dB_Q7 * ( 1 / 128.0f ) + 6.0f ) * ( 1.0f - psEncCtrl->input_quality );
1.173 + }
1.174 +
1.175 + /*************************/
1.176 + /* SPARSENESS PROCESSING */
1.177 + /*************************/
1.178 + /* Set quantizer offset */
1.179 + if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
1.180 + /* Initially set to 0; may be overruled in process_gains(..) */
1.181 + psEnc->sCmn.indices.quantOffsetType = 0;
1.182 + psEncCtrl->sparseness = 0.0f;
1.183 + } else {
1.184 + /* Sparseness measure, based on relative fluctuations of energy per 2 milliseconds */
1.185 + nSamples = 2 * psEnc->sCmn.fs_kHz;
1.186 + energy_variation = 0.0f;
1.187 + log_energy_prev = 0.0f;
1.188 + pitch_res_ptr = pitch_res;
1.189 + for( k = 0; k < silk_SMULBB( SUB_FRAME_LENGTH_MS, psEnc->sCmn.nb_subfr ) / 2; k++ ) {
1.190 + nrg = ( silk_float )nSamples + ( silk_float )silk_energy_FLP( pitch_res_ptr, nSamples );
1.191 + log_energy = silk_log2( nrg );
1.192 + if( k > 0 ) {
1.193 + energy_variation += silk_abs_float( log_energy - log_energy_prev );
1.194 + }
1.195 + log_energy_prev = log_energy;
1.196 + pitch_res_ptr += nSamples;
1.197 + }
1.198 + psEncCtrl->sparseness = silk_sigmoid( 0.4f * ( energy_variation - 5.0f ) );
1.199 +
1.200 + /* Set quantization offset depending on sparseness measure */
1.201 + if( psEncCtrl->sparseness > SPARSENESS_THRESHOLD_QNT_OFFSET ) {
1.202 + psEnc->sCmn.indices.quantOffsetType = 0;
1.203 + } else {
1.204 + psEnc->sCmn.indices.quantOffsetType = 1;
1.205 + }
1.206 +
1.207 + /* Increase coding SNR for sparse signals */
1.208 + SNR_adj_dB += SPARSE_SNR_INCR_dB * ( psEncCtrl->sparseness - 0.5f );
1.209 + }
1.210 +
1.211 + /*******************************/
1.212 + /* Control bandwidth expansion */
1.213 + /*******************************/
1.214 + /* More BWE for signals with high prediction gain */
1.215 + strength = FIND_PITCH_WHITE_NOISE_FRACTION * psEncCtrl->predGain; /* between 0.0 and 1.0 */
1.216 + BWExp1 = BWExp2 = BANDWIDTH_EXPANSION / ( 1.0f + strength * strength );
1.217 + delta = LOW_RATE_BANDWIDTH_EXPANSION_DELTA * ( 1.0f - 0.75f * psEncCtrl->coding_quality );
1.218 + BWExp1 -= delta;
1.219 + BWExp2 += delta;
1.220 + /* BWExp1 will be applied after BWExp2, so make it relative */
1.221 + BWExp1 /= BWExp2;
1.222 +
1.223 + if( psEnc->sCmn.warping_Q16 > 0 ) {
1.224 + /* Slightly more warping in analysis will move quantization noise up in frequency, where it's better masked */
1.225 + warping = (silk_float)psEnc->sCmn.warping_Q16 / 65536.0f + 0.01f * psEncCtrl->coding_quality;
1.226 + } else {
1.227 + warping = 0.0f;
1.228 + }
1.229 +
1.230 + /********************************************/
1.231 + /* Compute noise shaping AR coefs and gains */
1.232 + /********************************************/
1.233 + for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
1.234 + /* Apply window: sine slope followed by flat part followed by cosine slope */
1.235 + opus_int shift, slope_part, flat_part;
1.236 + flat_part = psEnc->sCmn.fs_kHz * 3;
1.237 + slope_part = ( psEnc->sCmn.shapeWinLength - flat_part ) / 2;
1.238 +
1.239 + silk_apply_sine_window_FLP( x_windowed, x_ptr, 1, slope_part );
1.240 + shift = slope_part;
1.241 + silk_memcpy( x_windowed + shift, x_ptr + shift, flat_part * sizeof(silk_float) );
1.242 + shift += flat_part;
1.243 + silk_apply_sine_window_FLP( x_windowed + shift, x_ptr + shift, 2, slope_part );
1.244 +
1.245 + /* Update pointer: next LPC analysis block */
1.246 + x_ptr += psEnc->sCmn.subfr_length;
1.247 +
1.248 + if( psEnc->sCmn.warping_Q16 > 0 ) {
1.249 + /* Calculate warped auto correlation */
1.250 + silk_warped_autocorrelation_FLP( auto_corr, x_windowed, warping,
1.251 + psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder );
1.252 + } else {
1.253 + /* Calculate regular auto correlation */
1.254 + silk_autocorrelation_FLP( auto_corr, x_windowed, psEnc->sCmn.shapeWinLength, psEnc->sCmn.shapingLPCOrder + 1 );
1.255 + }
1.256 +
1.257 + /* Add white noise, as a fraction of energy */
1.258 + auto_corr[ 0 ] += auto_corr[ 0 ] * SHAPE_WHITE_NOISE_FRACTION;
1.259 +
1.260 + /* Convert correlations to prediction coefficients, and compute residual energy */
1.261 + nrg = silk_levinsondurbin_FLP( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], auto_corr, psEnc->sCmn.shapingLPCOrder );
1.262 + psEncCtrl->Gains[ k ] = ( silk_float )sqrt( nrg );
1.263 +
1.264 + if( psEnc->sCmn.warping_Q16 > 0 ) {
1.265 + /* Adjust gain for warping */
1.266 + psEncCtrl->Gains[ k ] *= warped_gain( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], warping, psEnc->sCmn.shapingLPCOrder );
1.267 + }
1.268 +
1.269 + /* Bandwidth expansion for synthesis filter shaping */
1.270 + silk_bwexpander_FLP( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder, BWExp2 );
1.271 +
1.272 + /* Compute noise shaping filter coefficients */
1.273 + silk_memcpy(
1.274 + &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ],
1.275 + &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ],
1.276 + psEnc->sCmn.shapingLPCOrder * sizeof( silk_float ) );
1.277 +
1.278 + /* Bandwidth expansion for analysis filter shaping */
1.279 + silk_bwexpander_FLP( &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder, BWExp1 );
1.280 +
1.281 + /* Ratio of prediction gains, in energy domain */
1.282 + pre_nrg = silk_LPC_inverse_pred_gain_FLP( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder );
1.283 + nrg = silk_LPC_inverse_pred_gain_FLP( &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ], psEnc->sCmn.shapingLPCOrder );
1.284 + psEncCtrl->GainsPre[ k ] = 1.0f - 0.7f * ( 1.0f - pre_nrg / nrg );
1.285 +
1.286 + /* Convert to monic warped prediction coefficients and limit absolute values */
1.287 + warped_true2monic_coefs( &psEncCtrl->AR2[ k * MAX_SHAPE_LPC_ORDER ], &psEncCtrl->AR1[ k * MAX_SHAPE_LPC_ORDER ],
1.288 + warping, 3.999f, psEnc->sCmn.shapingLPCOrder );
1.289 + }
1.290 +
1.291 + /*****************/
1.292 + /* Gain tweaking */
1.293 + /*****************/
1.294 + /* Increase gains during low speech activity */
1.295 + gain_mult = (silk_float)pow( 2.0f, -0.16f * SNR_adj_dB );
1.296 + gain_add = (silk_float)pow( 2.0f, 0.16f * MIN_QGAIN_DB );
1.297 + for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
1.298 + psEncCtrl->Gains[ k ] *= gain_mult;
1.299 + psEncCtrl->Gains[ k ] += gain_add;
1.300 + }
1.301 +
1.302 + gain_mult = 1.0f + INPUT_TILT + psEncCtrl->coding_quality * HIGH_RATE_INPUT_TILT;
1.303 + for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
1.304 + psEncCtrl->GainsPre[ k ] *= gain_mult;
1.305 + }
1.306 +
1.307 + /************************************************/
1.308 + /* Control low-frequency shaping and noise tilt */
1.309 + /************************************************/
1.310 + /* Less low frequency shaping for noisy inputs */
1.311 + 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 ) );
1.312 + strength *= psEnc->sCmn.speech_activity_Q8 * ( 1.0f / 256.0f );
1.313 + if( psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
1.314 + /* Reduce low frequencies quantization noise for periodic signals, depending on pitch lag */
1.315 + /*f = 400; freqz([1, -0.98 + 2e-4 * f], [1, -0.97 + 7e-4 * f], 2^12, Fs); axis([0, 1000, -10, 1])*/
1.316 + for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
1.317 + b = 0.2f / psEnc->sCmn.fs_kHz + 3.0f / psEncCtrl->pitchL[ k ];
1.318 + psEncCtrl->LF_MA_shp[ k ] = -1.0f + b;
1.319 + psEncCtrl->LF_AR_shp[ k ] = 1.0f - b - b * strength;
1.320 + }
1.321 + Tilt = - HP_NOISE_COEF -
1.322 + (1 - HP_NOISE_COEF) * HARM_HP_NOISE_COEF * psEnc->sCmn.speech_activity_Q8 * ( 1.0f / 256.0f );
1.323 + } else {
1.324 + b = 1.3f / psEnc->sCmn.fs_kHz;
1.325 + psEncCtrl->LF_MA_shp[ 0 ] = -1.0f + b;
1.326 + psEncCtrl->LF_AR_shp[ 0 ] = 1.0f - b - b * strength * 0.6f;
1.327 + for( k = 1; k < psEnc->sCmn.nb_subfr; k++ ) {
1.328 + psEncCtrl->LF_MA_shp[ k ] = psEncCtrl->LF_MA_shp[ 0 ];
1.329 + psEncCtrl->LF_AR_shp[ k ] = psEncCtrl->LF_AR_shp[ 0 ];
1.330 + }
1.331 + Tilt = -HP_NOISE_COEF;
1.332 + }
1.333 +
1.334 + /****************************/
1.335 + /* HARMONIC SHAPING CONTROL */
1.336 + /****************************/
1.337 + /* Control boosting of harmonic frequencies */
1.338 + HarmBoost = LOW_RATE_HARMONIC_BOOST * ( 1.0f - psEncCtrl->coding_quality ) * psEnc->LTPCorr;
1.339 +
1.340 + /* More harmonic boost for noisy input signals */
1.341 + HarmBoost += LOW_INPUT_QUALITY_HARMONIC_BOOST * ( 1.0f - psEncCtrl->input_quality );
1.342 +
1.343 + if( USE_HARM_SHAPING && psEnc->sCmn.indices.signalType == TYPE_VOICED ) {
1.344 + /* Harmonic noise shaping */
1.345 + HarmShapeGain = HARMONIC_SHAPING;
1.346 +
1.347 + /* More harmonic noise shaping for high bitrates or noisy input */
1.348 + HarmShapeGain += HIGH_RATE_OR_LOW_QUALITY_HARMONIC_SHAPING *
1.349 + ( 1.0f - ( 1.0f - psEncCtrl->coding_quality ) * psEncCtrl->input_quality );
1.350 +
1.351 + /* Less harmonic noise shaping for less periodic signals */
1.352 + HarmShapeGain *= ( silk_float )sqrt( psEnc->LTPCorr );
1.353 + } else {
1.354 + HarmShapeGain = 0.0f;
1.355 + }
1.356 +
1.357 + /*************************/
1.358 + /* Smooth over subframes */
1.359 + /*************************/
1.360 + for( k = 0; k < psEnc->sCmn.nb_subfr; k++ ) {
1.361 + psShapeSt->HarmBoost_smth += SUBFR_SMTH_COEF * ( HarmBoost - psShapeSt->HarmBoost_smth );
1.362 + psEncCtrl->HarmBoost[ k ] = psShapeSt->HarmBoost_smth;
1.363 + psShapeSt->HarmShapeGain_smth += SUBFR_SMTH_COEF * ( HarmShapeGain - psShapeSt->HarmShapeGain_smth );
1.364 + psEncCtrl->HarmShapeGain[ k ] = psShapeSt->HarmShapeGain_smth;
1.365 + psShapeSt->Tilt_smth += SUBFR_SMTH_COEF * ( Tilt - psShapeSt->Tilt_smth );
1.366 + psEncCtrl->Tilt[ k ] = psShapeSt->Tilt_smth;
1.367 + }
1.368 +}
