/* * Copyright (c) 2012 Andrew D'Addesio * Copyright (c) 2013-2014 Mozilla Corporation * * This file is part of FFmpeg. * * FFmpeg is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * FFmpeg is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with FFmpeg; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ /** * @file * Opus SILK decoder */ #include <stdint.h> #include "libavutil/mem.h"#include "mathops.h"#include "opus.h"#include "rc.h"#include "silk.h"#include "tab.h" #define ROUND_MULL(a,b,s) (((MUL64(a, b) >> ((s) - 1)) + 1) >> 1) typedef struct SilkFrame {    int coded;    int log_gain;    int16_t nlsf[16];    float    lpc[16];     float output     [2 * SILK_HISTORY];    float lpc_history[2 * SILK_HISTORY];    int primarylag;     int prev_voiced;} SilkFrame; struct SilkContext {    void *logctx;    int output_channels;     int midonly;    int subframes;    int sflength;    int flength;    int nlsf_interp_factor;     enum OpusBandwidth bandwidth;    int wb;     SilkFrame frame[2];    float prev_stereo_weights[2];    float stereo_weights[2];     int prev_coded_channels;}; static inline void silk_stabilize_lsf(int16_t nlsf[16], int order, const uint16_t min_delta[17]){    int pass, i;    for (pass = 0; pass < 20; pass++) {        int k, min_diff = 0;        for (i = 0; i < order+1; i++) {            int low  = i != 0     ? nlsf[i-1] : 0;            int high = i != order ? nlsf[i]   : 32768;            int diff = (high - low) - (min_delta[i]);             if (diff < min_diff) {                min_diff = diff;                k = i;                 if (pass == 20)                    break;            }        }        if (min_diff == 0) /* no issues; stabilized */            return;         /* wiggle one or two LSFs */        if (k == 0) {            /* repel away from lower bound */            nlsf[0] = min_delta[0];        } else if (k == order) {            /* repel away from higher bound */            nlsf[order-1] = 32768 - min_delta[order];        } else {            /* repel away from current position */            int min_center = 0, max_center = 32768, center_val;             /* lower extent */            for (i = 0; i < k; i++)                min_center += min_delta[i];            min_center += min_delta[k] >> 1;             /* upper extent */            for (i = order; i > k; i--)                max_center -= min_delta[i];            max_center -= min_delta[k] >> 1;             /* move apart */            center_val = nlsf[k - 1] + nlsf[k];            center_val = (center_val >> 1) + (center_val & 1); // rounded divide by 2            center_val = FFMIN(max_center, FFMAX(min_center, center_val));             nlsf[k - 1] = center_val - (min_delta[k] >> 1);            nlsf[k]     = nlsf[k - 1] + min_delta[k];        }    }     /* resort to the fall-back method, the standard method for LSF stabilization */     /* sort; as the LSFs should be nearly sorted, use insertion sort */    for (i = 1; i < order; i++) {        int j, value = nlsf[i];        for (j = i - 1; j >= 0 && nlsf[j] > value; j--)            nlsf[j + 1] = nlsf[j];        nlsf[j + 1] = value;    }     /* push forwards to increase distance */    if (nlsf[0] < min_delta[0])        nlsf[0] = min_delta[0];    for (i = 1; i < order; i++)        nlsf[i] = FFMAX(nlsf[i], FFMIN(nlsf[i - 1] + min_delta[i], 32767));     /* push backwards to increase distance */    if (nlsf[order-1] > 32768 - min_delta[order])        nlsf[order-1] = 32768 - min_delta[order];    for (i = order-2; i >= 0; i--)        if (nlsf[i] > nlsf[i + 1] - min_delta[i+1])            nlsf[i] = nlsf[i + 1] - min_delta[i+1];     return;} static inline int silk_is_lpc_stable(const int16_t lpc[16], int order){    int k, j, DC_resp = 0;    int32_t lpc32[2][16];       // Q24    int totalinvgain = 1 << 30; // 1.0 in Q30    int32_t *row = lpc32[0], *prevrow;     /* initialize the first row for the Levinson recursion */    for (k = 0; k < order; k++) {        DC_resp += lpc[k];        row[k] = lpc[k] * 4096;    }     if (DC_resp >= 4096)        return 0;     /* check if prediction gain pushes any coefficients too far */    for (k = order - 1; 1; k--) {        int rc;      // Q31; reflection coefficient        int gaindiv; // Q30; inverse of the gain (the divisor)        int gain;    // gain for this reflection coefficient        int fbits;   // fractional bits used for the gain        int error;   // Q29; estimate of the error of our partial estimate of 1/gaindiv         if (FFABS(row[k]) > 16773022)            return 0;         rc      = -(row[k] * 128);        gaindiv = (1 << 30) - MULH(rc, rc);         totalinvgain = MULH(totalinvgain, gaindiv) << 2;        if (k == 0)            return (totalinvgain >= 107374);         /* approximate 1.0/gaindiv */        fbits = opus_ilog(gaindiv);        gain  = ((1 << 29) - 1) / (gaindiv >> (fbits + 1 - 16)); // Q<fbits-16>        error = (1 << 29) - MULL(gaindiv << (15 + 16 - fbits), gain, 16);        gain  = ((gain << 16) + (error * gain >> 13));         /* switch to the next row of the LPC coefficients */        prevrow = row;        row = lpc32[k & 1];         for (j = 0; j < k; j++) {            int x = av_sat_sub32(prevrow[j], ROUND_MULL(prevrow[k - j - 1], rc, 31));            int64_t tmp = ROUND_MULL(x, gain, fbits);             /* per RFC 8251 section 6, if this calculation overflows, the filter               is considered unstable. */            if (tmp < INT32_MIN || tmp > INT32_MAX)                return 0;             row[j] = (int32_t)tmp;        }    }} static void silk_lsp2poly(const int32_t lsp[/* 2 * half_order - 1 */],                          int32_t pol[/* half_order + 1 */], int half_order){    int i, j;     pol[0] = 65536; // 1.0 in Q16    pol[1] = -lsp[0];     for (i = 1; i < half_order; i++) {        pol[i + 1] = pol[i - 1] * 2 - ROUND_MULL(lsp[2 * i], pol[i], 16);        for (j = i; j > 1; j--)            pol[j] += pol[j - 2] - ROUND_MULL(lsp[2 * i], pol[j - 1], 16);         pol[1] -= lsp[2 * i];    }} static void silk_lsf2lpc(const int16_t nlsf[16], float lpcf[16], int order){    int i, k;    int32_t lsp[16];     // Q17; 2*cos(LSF)    int32_t p[9], q[9];  // Q16    int32_t lpc32[16];   // Q17    int16_t lpc[16];     // Q12     /* convert the LSFs to LSPs, i.e. 2*cos(LSF) */    for (k = 0; k < order; k++) {        int index = nlsf[k] >> 8;        int offset = nlsf[k] & 255;        int k2 = (order == 10) ? ff_silk_lsf_ordering_nbmb[k] : ff_silk_lsf_ordering_wb[k];         /* interpolate and round */        lsp[k2]  = ff_silk_cosine[index] * 256;        lsp[k2] += (ff_silk_cosine[index + 1] - ff_silk_cosine[index]) * offset;        lsp[k2]  = (lsp[k2] + 4) >> 3;    }     silk_lsp2poly(lsp    , p, order >> 1);    silk_lsp2poly(lsp + 1, q, order >> 1);     /* reconstruct A(z) */    for (k = 0; k < order>>1; k++) {        int32_t p_tmp = p[k + 1] + p[k];        int32_t q_tmp = q[k + 1] - q[k];        lpc32[k]         = -q_tmp - p_tmp;        lpc32[order-k-1] =  q_tmp - p_tmp;    }     /* limit the range of the LPC coefficients to each fit within an int16_t */    for (i = 0; i < 10; i++) {        int j;        unsigned int maxabs = 0;        for (j = 0, k = 0; j < order; j++) {            unsigned int x = FFABS(lpc32[k]);            if (x > maxabs) {                maxabs = x; // Q17                k      = j;            }        }         maxabs = (maxabs + 16) >> 5; // convert to Q12         if (maxabs > 32767) {            /* perform bandwidth expansion */            unsigned int chirp, chirp_base; // Q16            maxabs = FFMIN(maxabs, 163838); // anything above this overflows chirp's numerator            chirp_base = chirp = 65470 - ((maxabs - 32767) << 14) / ((maxabs * (k+1)) >> 2);             for (k = 0; k < order; k++) {                lpc32[k] = ROUND_MULL(lpc32[k], chirp, 16);                chirp    = (chirp_base * chirp + 32768) >> 16;            }        } else break;    }     if (i == 10) {        /* time's up: just clamp */        for (k = 0; k < order; k++) {            int x = (lpc32[k] + 16) >> 5;            lpc[k] = av_clip_int16(x);            lpc32[k] = lpc[k] << 5; // shortcut mandated by the spec; drops lower 5 bits        }    } else {        for (k = 0; k < order; k++)            lpc[k] = (lpc32[k] + 16) >> 5;    }     /* if the prediction gain causes the LPC filter to become unstable,       apply further bandwidth expansion on the Q17 coefficients */    for (i = 1; i <= 16 && !silk_is_lpc_stable(lpc, order); i++) {        unsigned int chirp, chirp_base;        chirp_base = chirp = 65536 - (1 << i);         for (k = 0; k < order; k++) {            lpc32[k] = ROUND_MULL(lpc32[k], chirp, 16);            lpc[k]   = (lpc32[k] + 16) >> 5;            chirp    = (chirp_base * chirp + 32768) >> 16;        }    }     for (i = 0; i < order; i++)        lpcf[i] = lpc[i] / 4096.0f;} static inline void silk_decode_lpc(SilkContext *s, SilkFrame *frame,                                   OpusRangeCoder *rc,                                   float lpc_leadin[16], float lpc[16],                                   int *lpc_order, int *has_lpc_leadin, int voiced){    int i;    int order;                   // order of the LP polynomial; 10 for NB/MB and 16 for WB    int8_t  lsf_i1, lsf_i2[16];  // stage-1 and stage-2 codebook indices    int16_t lsf_res[16];         // residual as a Q10 value    int16_t nlsf[16];            // Q15     *lpc_order = order = s->wb ? 16 : 10;     /* obtain LSF stage-1 and stage-2 indices */    lsf_i1 = ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s1[s->wb][voiced]);    for (i = 0; i < order; i++) {        int index = s->wb ? ff_silk_lsf_s2_model_sel_wb  [lsf_i1][i] :                            ff_silk_lsf_s2_model_sel_nbmb[lsf_i1][i];        lsf_i2[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s2[index]) - 4;        if (lsf_i2[i] == -4)            lsf_i2[i] -= ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s2_ext);        else if (lsf_i2[i] == 4)            lsf_i2[i] += ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s2_ext);    }     /* reverse the backwards-prediction step */    for (i = order - 1; i >= 0; i--) {        int qstep = s->wb ? 9830 : 11796;         lsf_res[i] = lsf_i2[i] * 1024;        if (lsf_i2[i] < 0)      lsf_res[i] += 102;        else if (lsf_i2[i] > 0) lsf_res[i] -= 102;        lsf_res[i] = (lsf_res[i] * qstep) >> 16;         if (i + 1 < order) {            int weight = s->wb ? ff_silk_lsf_pred_weights_wb  [ff_silk_lsf_weight_sel_wb  [lsf_i1][i]][i] :                                 ff_silk_lsf_pred_weights_nbmb[ff_silk_lsf_weight_sel_nbmb[lsf_i1][i]][i];            lsf_res[i] += (lsf_res[i+1] * weight) >> 8;        }    }     /* reconstruct the NLSF coefficients from the supplied indices */    for (i = 0; i < order; i++) {        const uint8_t * codebook = s->wb ? ff_silk_lsf_codebook_wb  [lsf_i1] :                                           ff_silk_lsf_codebook_nbmb[lsf_i1];        int cur, prev, next, weight_sq, weight, ipart, fpart, y, value;         /* find the weight of the residual */        /* TODO: precompute */        cur = codebook[i];        prev = i ? codebook[i - 1] : 0;        next = i + 1 < order ? codebook[i + 1] : 256;        weight_sq = (1024 / (cur - prev) + 1024 / (next - cur)) << 16;         /* approximate square-root with mandated fixed-point arithmetic */        ipart = opus_ilog(weight_sq);        fpart = (weight_sq >> (ipart-8)) & 127;        y = ((ipart & 1) ? 32768 : 46214) >> ((32 - ipart)>>1);        weight = y + ((213 * fpart * y) >> 16);         value = cur * 128 + (lsf_res[i] * 16384) / weight;        nlsf[i] = av_clip_uintp2(value, 15);    }     /* stabilize the NLSF coefficients */    silk_stabilize_lsf(nlsf, order, s->wb ? ff_silk_lsf_min_spacing_wb :                                            ff_silk_lsf_min_spacing_nbmb);     /* produce an interpolation for the first 2 subframes, */    /* and then convert both sets of NLSFs to LPC coefficients */    *has_lpc_leadin = 0;    if (s->subframes == 4) {        int offset = ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_interpolation_offset);        if (offset != 4 && frame->coded) {            *has_lpc_leadin = 1;            if (offset != 0) {                int16_t nlsf_leadin[16];                for (i = 0; i < order; i++)                    nlsf_leadin[i] = frame->nlsf[i] +                        ((nlsf[i] - frame->nlsf[i]) * offset >> 2);                silk_lsf2lpc(nlsf_leadin, lpc_leadin, order);            } else  /* avoid re-computation for a (roughly) 1-in-4 occurrence */                memcpy(lpc_leadin, frame->lpc, 16 * sizeof(float));        } else            offset = 4;        s->nlsf_interp_factor = offset;         silk_lsf2lpc(nlsf, lpc, order);    } else {        s->nlsf_interp_factor = 4;        silk_lsf2lpc(nlsf, lpc, order);    }     memcpy(frame->nlsf, nlsf, order * sizeof(nlsf[0]));    memcpy(frame->lpc,  lpc,  order * sizeof(lpc[0]));} static inline void silk_count_children(OpusRangeCoder *rc, int model, int32_t total,                                       int32_t child[2]){    if (total != 0) {        child[0] = ff_opus_rc_dec_cdf(rc,                       ff_silk_model_pulse_location[model] + (((total - 1 + 5) * (total - 1)) >> 1));        child[1] = total - child[0];    } else {        child[0] = 0;        child[1] = 0;    }} static inline void silk_decode_excitation(SilkContext *s, OpusRangeCoder *rc,                                          float* excitationf,                                          int qoffset_high, int active, int voiced){    int i;    uint32_t seed;    int shellblocks;    int ratelevel;    uint8_t pulsecount[20];     // total pulses in each shell block    uint8_t lsbcount[20] = {0}; // raw lsbits defined for each pulse in each shell block    int32_t excitation[320];    // Q23     /* excitation parameters */    seed = ff_opus_rc_dec_cdf(rc, ff_silk_model_lcg_seed);    shellblocks = ff_silk_shell_blocks[s->bandwidth][s->subframes >> 2];    ratelevel = ff_opus_rc_dec_cdf(rc, ff_silk_model_exc_rate[voiced]);     for (i = 0; i < shellblocks; i++) {        pulsecount[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_pulse_count[ratelevel]);        if (pulsecount[i] == 17) {            while (pulsecount[i] == 17 && ++lsbcount[i] != 10)                pulsecount[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_pulse_count[9]);            if (lsbcount[i] == 10)                pulsecount[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_pulse_count[10]);        }    }     /* decode pulse locations using PVQ */    for (i = 0; i < shellblocks; i++) {        if (pulsecount[i] != 0) {            int a, b, c, d;            int32_t * location = excitation + 16*i;            int32_t branch[4][2];            branch[0][0] = pulsecount[i];             /* unrolled tail recursion */            for (a = 0; a < 1; a++) {                silk_count_children(rc, 0, branch[0][a], branch[1]);                for (b = 0; b < 2; b++) {                    silk_count_children(rc, 1, branch[1][b], branch[2]);                    for (c = 0; c < 2; c++) {                        silk_count_children(rc, 2, branch[2][c], branch[3]);                        for (d = 0; d < 2; d++) {                            silk_count_children(rc, 3, branch[3][d], location);                            location += 2;                        }                    }                }            }        } else            memset(excitation + 16*i, 0, 16*sizeof(int32_t));    }     /* decode least significant bits */    for (i = 0; i < shellblocks << 4; i++) {        int bit;        for (bit = 0; bit < lsbcount[i >> 4]; bit++)            excitation[i] = (excitation[i] << 1) |                            ff_opus_rc_dec_cdf(rc, ff_silk_model_excitation_lsb);    }     /* decode signs */    for (i = 0; i < shellblocks << 4; i++) {        if (excitation[i] != 0) {            int sign = ff_opus_rc_dec_cdf(rc, ff_silk_model_excitation_sign[active +                                         voiced][qoffset_high][FFMIN(pulsecount[i >> 4], 6)]);            if (sign == 0)                excitation[i] *= -1;        }    }     /* assemble the excitation */    for (i = 0; i < shellblocks << 4; i++) {        int value = excitation[i];        excitation[i] = value * 256 | ff_silk_quant_offset[voiced][qoffset_high];        if (value < 0)      excitation[i] += 20;        else if (value > 0) excitation[i] -= 20;         /* invert samples pseudorandomly */        seed = 196314165 * seed + 907633515;

Lines 1–500 of 905. Load more