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