LCOV - code coverage report
Current view: top level - libavcodec - opus_silk.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 391 428 91.4 %
Date: 2017-12-14 19:11:59 Functions: 14 14 100.0 %

          Line data    Source code
       1             : /*
       2             :  * Copyright (c) 2012 Andrew D'Addesio
       3             :  * Copyright (c) 2013-2014 Mozilla Corporation
       4             :  *
       5             :  * This file is part of FFmpeg.
       6             :  *
       7             :  * FFmpeg is free software; you can redistribute it and/or
       8             :  * modify it under the terms of the GNU Lesser General Public
       9             :  * License as published by the Free Software Foundation; either
      10             :  * version 2.1 of the License, or (at your option) any later version.
      11             :  *
      12             :  * FFmpeg is distributed in the hope that it will be useful,
      13             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      14             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      15             :  * Lesser General Public License for more details.
      16             :  *
      17             :  * You should have received a copy of the GNU Lesser General Public
      18             :  * License along with FFmpeg; if not, write to the Free Software
      19             :  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
      20             :  */
      21             : 
      22             : /**
      23             :  * @file
      24             :  * Opus SILK decoder
      25             :  */
      26             : 
      27             : #include <stdint.h>
      28             : 
      29             : #include "opus.h"
      30             : #include "opustab.h"
      31             : 
      32             : typedef struct SilkFrame {
      33             :     int coded;
      34             :     int log_gain;
      35             :     int16_t nlsf[16];
      36             :     float    lpc[16];
      37             : 
      38             :     float output     [2 * SILK_HISTORY];
      39             :     float lpc_history[2 * SILK_HISTORY];
      40             :     int primarylag;
      41             : 
      42             :     int prev_voiced;
      43             : } SilkFrame;
      44             : 
      45             : struct SilkContext {
      46             :     AVCodecContext *avctx;
      47             :     int output_channels;
      48             : 
      49             :     int midonly;
      50             :     int subframes;
      51             :     int sflength;
      52             :     int flength;
      53             :     int nlsf_interp_factor;
      54             : 
      55             :     enum OpusBandwidth bandwidth;
      56             :     int wb;
      57             : 
      58             :     SilkFrame frame[2];
      59             :     float prev_stereo_weights[2];
      60             :     float stereo_weights[2];
      61             : 
      62             :     int prev_coded_channels;
      63             : };
      64             : 
      65       13297 : static inline void silk_stabilize_lsf(int16_t nlsf[16], int order, const uint16_t min_delta[17])
      66             : {
      67             :     int pass, i;
      68       13683 :     for (pass = 0; pass < 20; pass++) {
      69       13683 :         int k, min_diff = 0;
      70      218526 :         for (i = 0; i < order+1; i++) {
      71      204843 :             int low  = i != 0     ? nlsf[i-1] : 0;
      72      204843 :             int high = i != order ? nlsf[i]   : 32768;
      73      204843 :             int diff = (high - low) - (min_delta[i]);
      74             : 
      75      204843 :             if (diff < min_diff) {
      76         387 :                 min_diff = diff;
      77         387 :                 k = i;
      78             : 
      79         387 :                 if (pass == 20)
      80           0 :                     break;
      81             :             }
      82             :         }
      83       13683 :         if (min_diff == 0) /* no issues; stabilized */
      84       13297 :             return;
      85             : 
      86             :         /* wiggle one or two LSFs */
      87         386 :         if (k == 0) {
      88             :             /* repel away from lower bound */
      89          74 :             nlsf[0] = min_delta[0];
      90         312 :         } else if (k == order) {
      91             :             /* repel away from higher bound */
      92           0 :             nlsf[order-1] = 32768 - min_delta[order];
      93             :         } else {
      94             :             /* repel away from current position */
      95         312 :             int min_center = 0, max_center = 32768, center_val;
      96             : 
      97             :             /* lower extent */
      98        1429 :             for (i = 0; i < k; i++)
      99        1117 :                 min_center += min_delta[i];
     100         312 :             min_center += min_delta[k] >> 1;
     101             : 
     102             :             /* upper extent */
     103        3689 :             for (i = order; i > k; i--)
     104        3377 :                 max_center -= min_delta[i];
     105         312 :             max_center -= min_delta[k] >> 1;
     106             : 
     107             :             /* move apart */
     108         312 :             center_val = nlsf[k - 1] + nlsf[k];
     109         312 :             center_val = (center_val >> 1) + (center_val & 1); // rounded divide by 2
     110         312 :             center_val = FFMIN(max_center, FFMAX(min_center, center_val));
     111             : 
     112         312 :             nlsf[k - 1] = center_val - (min_delta[k] >> 1);
     113         312 :             nlsf[k]     = nlsf[k - 1] + min_delta[k];
     114             :         }
     115             :     }
     116             : 
     117             :     /* resort to the fall-back method, the standard method for LSF stabilization */
     118             : 
     119             :     /* sort; as the LSFs should be nearly sorted, use insertion sort */
     120           0 :     for (i = 1; i < order; i++) {
     121           0 :         int j, value = nlsf[i];
     122           0 :         for (j = i - 1; j >= 0 && nlsf[j] > value; j--)
     123           0 :             nlsf[j + 1] = nlsf[j];
     124           0 :         nlsf[j + 1] = value;
     125             :     }
     126             : 
     127             :     /* push forwards to increase distance */
     128           0 :     if (nlsf[0] < min_delta[0])
     129           0 :         nlsf[0] = min_delta[0];
     130           0 :     for (i = 1; i < order; i++)
     131           0 :         nlsf[i] = FFMAX(nlsf[i], FFMIN(nlsf[i - 1] + min_delta[i], 32767));
     132             : 
     133             :     /* push backwards to increase distance */
     134           0 :     if (nlsf[order-1] > 32768 - min_delta[order])
     135           0 :         nlsf[order-1] = 32768 - min_delta[order];
     136           0 :     for (i = order-2; i >= 0; i--)
     137           0 :         if (nlsf[i] > nlsf[i + 1] - min_delta[i+1])
     138           0 :             nlsf[i] = nlsf[i + 1] - min_delta[i+1];
     139             : 
     140           0 :     return;
     141             : }
     142             : 
     143       15533 : static inline int silk_is_lpc_stable(const int16_t lpc[16], int order)
     144             : {
     145       15533 :     int k, j, DC_resp = 0;
     146             :     int32_t lpc32[2][16];       // Q24
     147       15533 :     int totalinvgain = 1 << 30; // 1.0 in Q30
     148       15533 :     int32_t *row = lpc32[0], *prevrow;
     149             : 
     150             :     /* initialize the first row for the Levinson recursion */
     151      231301 :     for (k = 0; k < order; k++) {
     152      215768 :         DC_resp += lpc[k];
     153      215768 :         row[k] = lpc[k] * 4096;
     154             :     }
     155             : 
     156       15533 :     if (DC_resp >= 4096)
     157           0 :         return 0;
     158             : 
     159             :     /* check if prediction gain pushes any coefficients too far */
     160      215699 :     for (k = order - 1; 1; k--) {
     161             :         int rc;      // Q31; reflection coefficient
     162             :         int gaindiv; // Q30; inverse of the gain (the divisor)
     163             :         int gain;    // gain for this reflection coefficient
     164             :         int fbits;   // fractional bits used for the gain
     165             :         int error;   // Q29; estimate of the error of our partial estimate of 1/gaindiv
     166             : 
     167      215699 :         if (FFABS(row[k]) > 16773022)
     168          69 :             return 0;
     169             : 
     170      215630 :         rc      = -(row[k] * 128);
     171      215630 :         gaindiv = (1 << 30) - MULH(rc, rc);
     172             : 
     173      215630 :         totalinvgain = MULH(totalinvgain, gaindiv) << 2;
     174      215630 :         if (k == 0)
     175       15464 :             return (totalinvgain >= 107374);
     176             : 
     177             :         /* approximate 1.0/gaindiv */
     178      200166 :         fbits = opus_ilog(gaindiv);
     179      200166 :         gain  = ((1 << 29) - 1) / (gaindiv >> (fbits + 1 - 16)); // Q<fbits-16>
     180      200166 :         error = (1 << 29) - MULL(gaindiv << (15 + 16 - fbits), gain, 16);
     181      200166 :         gain  = ((gain << 16) + (error * gain >> 13));
     182             : 
     183             :         /* switch to the next row of the LPC coefficients */
     184      200166 :         prevrow = row;
     185      200166 :         row = lpc32[k & 1];
     186             : 
     187     1654557 :         for (j = 0; j < k; j++) {
     188     1454391 :             int x = av_sat_sub32(prevrow[j], ROUND_MULL(prevrow[k - j - 1], rc, 31));
     189     1454391 :             int64_t tmp = ROUND_MULL(x, gain, fbits);
     190             : 
     191             :             /* per RFC 8251 section 6, if this calculation overflows, the filter
     192             :                is considered unstable. */
     193     1454391 :             if (tmp < INT32_MIN || tmp > INT32_MAX)
     194           0 :                 return 0;
     195             : 
     196     1454391 :             row[j] = (int32_t)tmp;
     197             :         }
     198             :     }
     199             : }
     200             : 
     201       30668 : static void silk_lsp2poly(const int32_t lsp[16], int32_t pol[16], int half_order)
     202             : {
     203             :     int i, j;
     204             : 
     205       30668 :     pol[0] = 65536; // 1.0 in Q16
     206       30668 :     pol[1] = -lsp[0];
     207             : 
     208      212698 :     for (i = 1; i < half_order; i++) {
     209      182030 :         pol[i + 1] = pol[i - 1] * 2 - ROUND_MULL(lsp[2 * i], pol[i], 16);
     210      662828 :         for (j = i; j > 1; j--)
     211      480798 :             pol[j] += pol[j - 2] - ROUND_MULL(lsp[2 * i], pol[j - 1], 16);
     212             : 
     213      182030 :         pol[1] -= lsp[2 * i];
     214             :     }
     215       30668 : }
     216             : 
     217       15334 : static void silk_lsf2lpc(const int16_t nlsf[16], float lpcf[16], int order)
     218             : {
     219             :     int i, k;
     220             :     int32_t lsp[16];     // Q17; 2*cos(LSF)
     221             :     int32_t p[9], q[9];  // Q16
     222             :     int32_t lpc32[16];   // Q17
     223             :     int16_t lpc[16];     // Q12
     224             : 
     225             :     /* convert the LSFs to LSPs, i.e. 2*cos(LSF) */
     226      228032 :     for (k = 0; k < order; k++) {
     227      212698 :         int index = nlsf[k] >> 8;
     228      212698 :         int offset = nlsf[k] & 255;
     229      212698 :         int k2 = (order == 10) ? ff_silk_lsf_ordering_nbmb[k] : ff_silk_lsf_ordering_wb[k];
     230             : 
     231             :         /* interpolate and round */
     232      212698 :         lsp[k2]  = ff_silk_cosine[index] * 256;
     233      212698 :         lsp[k2] += (ff_silk_cosine[index + 1] - ff_silk_cosine[index]) * offset;
     234      212698 :         lsp[k2]  = (lsp[k2] + 4) >> 3;
     235             :     }
     236             : 
     237       15334 :     silk_lsp2poly(lsp    , p, order >> 1);
     238       15334 :     silk_lsp2poly(lsp + 1, q, order >> 1);
     239             : 
     240             :     /* reconstruct A(z) */
     241      121683 :     for (k = 0; k < order>>1; k++) {
     242      106349 :         lpc32[k]         = -p[k + 1] - p[k] - q[k + 1] + q[k];
     243      106349 :         lpc32[order-k-1] = -p[k + 1] - p[k] + q[k + 1] - q[k];
     244             :     }
     245             : 
     246             :     /* limit the range of the LPC coefficients to each fit within an int16_t */
     247       15334 :     for (i = 0; i < 10; i++) {
     248             :         int j;
     249       15334 :         unsigned int maxabs = 0;
     250      228032 :         for (j = 0, k = 0; j < order; j++) {
     251      212698 :             unsigned int x = FFABS(lpc32[k]);
     252      212698 :             if (x > maxabs) {
     253       15334 :                 maxabs = x; // Q17
     254       15334 :                 k      = j;
     255             :             }
     256             :         }
     257             : 
     258       15334 :         maxabs = (maxabs + 16) >> 5; // convert to Q12
     259             : 
     260       15334 :         if (maxabs > 32767) {
     261             :             /* perform bandwidth expansion */
     262             :             unsigned int chirp, chirp_base; // Q16
     263           0 :             maxabs = FFMIN(maxabs, 163838); // anything above this overflows chirp's numerator
     264           0 :             chirp_base = chirp = 65470 - ((maxabs - 32767) << 14) / ((maxabs * (k+1)) >> 2);
     265             : 
     266           0 :             for (k = 0; k < order; k++) {
     267           0 :                 lpc32[k] = ROUND_MULL(lpc32[k], chirp, 16);
     268           0 :                 chirp    = (chirp_base * chirp + 32768) >> 16;
     269             :             }
     270       15334 :         } else break;
     271             :     }
     272             : 
     273       15334 :     if (i == 10) {
     274             :         /* time's up: just clamp */
     275           0 :         for (k = 0; k < order; k++) {
     276           0 :             int x = (lpc32[k] + 16) >> 5;
     277           0 :             lpc[k] = av_clip_int16(x);
     278           0 :             lpc32[k] = lpc[k] << 5; // shortcut mandated by the spec; drops lower 5 bits
     279             :         }
     280             :     } else {
     281      228032 :         for (k = 0; k < order; k++)
     282      212698 :             lpc[k] = (lpc32[k] + 16) >> 5;
     283             :     }
     284             : 
     285             :     /* if the prediction gain causes the LPC filter to become unstable,
     286             :        apply further bandwidth expansion on the Q17 coefficients */
     287       15533 :     for (i = 1; i <= 16 && !silk_is_lpc_stable(lpc, order); i++) {
     288             :         unsigned int chirp, chirp_base;
     289         199 :         chirp_base = chirp = 65536 - (1 << i);
     290             : 
     291        3269 :         for (k = 0; k < order; k++) {
     292        3070 :             lpc32[k] = ROUND_MULL(lpc32[k], chirp, 16);
     293        3070 :             lpc[k]   = (lpc32[k] + 16) >> 5;
     294        3070 :             chirp    = (chirp_base * chirp + 32768) >> 16;
     295             :         }
     296             :     }
     297             : 
     298      228032 :     for (i = 0; i < order; i++)
     299      212698 :         lpcf[i] = lpc[i] / 4096.0f;
     300       15334 : }
     301             : 
     302       13297 : static inline void silk_decode_lpc(SilkContext *s, SilkFrame *frame,
     303             :                                    OpusRangeCoder *rc,
     304             :                                    float lpc_leadin[16], float lpc[16],
     305             :                                    int *lpc_order, int *has_lpc_leadin, int voiced)
     306             : {
     307             :     int i;
     308             :     int order;                   // order of the LP polynomial; 10 for NB/MB and 16 for WB
     309             :     int8_t  lsf_i1, lsf_i2[16];  // stage-1 and stage-2 codebook indices
     310             :     int16_t lsf_res[16];         // residual as a Q10 value
     311             :     int16_t nlsf[16];            // Q15
     312             : 
     313       13297 :     *lpc_order = order = s->wb ? 16 : 10;
     314             : 
     315             :     /* obtain LSF stage-1 and stage-2 indices */
     316       13297 :     lsf_i1 = ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s1[s->wb][voiced]);
     317      198977 :     for (i = 0; i < order; i++) {
     318      230800 :         int index = s->wb ? ff_silk_lsf_s2_model_sel_wb  [lsf_i1][i] :
     319       45120 :                             ff_silk_lsf_s2_model_sel_nbmb[lsf_i1][i];
     320      185680 :         lsf_i2[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s2[index]) - 4;
     321      185680 :         if (lsf_i2[i] == -4)
     322         460 :             lsf_i2[i] -= ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s2_ext);
     323      185220 :         else if (lsf_i2[i] == 4)
     324          15 :             lsf_i2[i] += ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_s2_ext);
     325             :     }
     326             : 
     327             :     /* reverse the backwards-prediction step */
     328      198977 :     for (i = order - 1; i >= 0; i--) {
     329      185680 :         int qstep = s->wb ? 9830 : 11796;
     330             : 
     331      185680 :         lsf_res[i] = lsf_i2[i] * 1024;
     332      185680 :         if (lsf_i2[i] < 0)      lsf_res[i] += 102;
     333      141071 :         else if (lsf_i2[i] > 0) lsf_res[i] -= 102;
     334      185680 :         lsf_res[i] = (lsf_res[i] * qstep) >> 16;
     335             : 
     336      185680 :         if (i + 1 < order) {
     337      212991 :             int weight = s->wb ? ff_silk_lsf_pred_weights_wb  [ff_silk_lsf_weight_sel_wb  [lsf_i1][i]][i] :
     338       40608 :                                  ff_silk_lsf_pred_weights_nbmb[ff_silk_lsf_weight_sel_nbmb[lsf_i1][i]][i];
     339      172383 :             lsf_res[i] += (lsf_res[i+1] * weight) >> 8;
     340             :         }
     341             :     }
     342             : 
     343             :     /* reconstruct the NLSF coefficients from the supplied indices */
     344      198977 :     for (i = 0; i < order; i++) {
     345      230800 :         const uint8_t * codebook = s->wb ? ff_silk_lsf_codebook_wb  [lsf_i1] :
     346       45120 :                                            ff_silk_lsf_codebook_nbmb[lsf_i1];
     347             :         int cur, prev, next, weight_sq, weight, ipart, fpart, y, value;
     348             : 
     349             :         /* find the weight of the residual */
     350             :         /* TODO: precompute */
     351      185680 :         cur = codebook[i];
     352      185680 :         prev = i ? codebook[i - 1] : 0;
     353      185680 :         next = i + 1 < order ? codebook[i + 1] : 256;
     354      185680 :         weight_sq = (1024 / (cur - prev) + 1024 / (next - cur)) << 16;
     355             : 
     356             :         /* approximate square-root with mandated fixed-point arithmetic */
     357      185680 :         ipart = opus_ilog(weight_sq);
     358      185680 :         fpart = (weight_sq >> (ipart-8)) & 127;
     359      185680 :         y = ((ipart & 1) ? 32768 : 46214) >> ((32 - ipart)>>1);
     360      185680 :         weight = y + ((213 * fpart * y) >> 16);
     361             : 
     362      185680 :         value = cur * 128 + (lsf_res[i] * 16384) / weight;
     363      185680 :         nlsf[i] = av_clip_uintp2(value, 15);
     364             :     }
     365             : 
     366             :     /* stabilize the NLSF coefficients */
     367       13297 :     silk_stabilize_lsf(nlsf, order, s->wb ? ff_silk_lsf_min_spacing_wb :
     368             :                                             ff_silk_lsf_min_spacing_nbmb);
     369             : 
     370             :     /* produce an interpolation for the first 2 subframes, */
     371             :     /* and then convert both sets of NLSFs to LPC coefficients */
     372       13297 :     *has_lpc_leadin = 0;
     373       13297 :     if (s->subframes == 4) {
     374        7616 :         int offset = ff_opus_rc_dec_cdf(rc, ff_silk_model_lsf_interpolation_offset);
     375        7616 :         if (offset != 4 && frame->coded) {
     376        2344 :             *has_lpc_leadin = 1;
     377        4688 :             if (offset != 0) {
     378             :                 int16_t nlsf_leadin[16];
     379       29055 :                 for (i = 0; i < order; i++)
     380       54036 :                     nlsf_leadin[i] = frame->nlsf[i] +
     381       27018 :                         ((nlsf[i] - frame->nlsf[i]) * offset >> 2);
     382        2037 :                 silk_lsf2lpc(nlsf_leadin, lpc_leadin, order);
     383             :             } else  /* avoid re-computation for a (roughly) 1-in-4 occurrence */
     384         307 :                 memcpy(lpc_leadin, frame->lpc, 16 * sizeof(float));
     385             :         } else
     386        5272 :             offset = 4;
     387        7616 :         s->nlsf_interp_factor = offset;
     388             : 
     389        7616 :         silk_lsf2lpc(nlsf, lpc, order);
     390             :     } else {
     391        5681 :         s->nlsf_interp_factor = 4;
     392        5681 :         silk_lsf2lpc(nlsf, lpc, order);
     393             :     }
     394             : 
     395       13297 :     memcpy(frame->nlsf, nlsf, order * sizeof(nlsf[0]));
     396       13297 :     memcpy(frame->lpc,  lpc,  order * sizeof(lpc[0]));
     397       13297 : }
     398             : 
     399     2319855 : static inline void silk_count_children(OpusRangeCoder *rc, int model, int32_t total,
     400             :                                        int32_t child[2])
     401             : {
     402     2319855 :     if (total != 0) {
     403     1295054 :         child[0] = ff_opus_rc_dec_cdf(rc,
     404     1295054 :                        ff_silk_model_pulse_location[model] + (((total - 1 + 5) * (total - 1)) >> 1));
     405     1295054 :         child[1] = total - child[0];
     406             :     } else {
     407     1024801 :         child[0] = 0;
     408     1024801 :         child[1] = 0;
     409             :     }
     410     2319855 : }
     411             : 
     412       13297 : static inline void silk_decode_excitation(SilkContext *s, OpusRangeCoder *rc,
     413             :                                           float* excitationf,
     414             :                                           int qoffset_high, int active, int voiced)
     415             : {
     416             :     int i;
     417             :     uint32_t seed;
     418             :     int shellblocks;
     419             :     int ratelevel;
     420             :     uint8_t pulsecount[20];     // total pulses in each shell block
     421       13297 :     uint8_t lsbcount[20] = {0}; // raw lsbits defined for each pulse in each shell block
     422             :     int32_t excitation[320];    // Q23
     423             : 
     424             :     /* excitation parameters */
     425       13297 :     seed = ff_opus_rc_dec_cdf(rc, ff_silk_model_lcg_seed);
     426       13297 :     shellblocks = ff_silk_shell_blocks[s->bandwidth][s->subframes >> 2];
     427       13297 :     ratelevel = ff_opus_rc_dec_cdf(rc, ff_silk_model_exc_rate[voiced]);
     428             : 
     429      192819 :     for (i = 0; i < shellblocks; i++) {
     430      179522 :         pulsecount[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_pulse_count[ratelevel]);
     431      179522 :         if (pulsecount[i] == 17) {
     432       33805 :             while (pulsecount[i] == 17 && ++lsbcount[i] != 10)
     433       13191 :                 pulsecount[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_pulse_count[9]);
     434       10307 :             if (lsbcount[i] == 10)
     435           0 :                 pulsecount[i] = ff_opus_rc_dec_cdf(rc, ff_silk_model_pulse_count[10]);
     436             :         }
     437             :     }
     438             : 
     439             :     /* decode pulse locations using PVQ */
     440      192819 :     for (i = 0; i < shellblocks; i++) {
     441      179522 :         if (pulsecount[i] != 0) {
     442             :             int a, b, c, d;
     443      154657 :             int32_t * location = excitation + 16*i;
     444             :             int32_t branch[4][2];
     445      154657 :             branch[0][0] = pulsecount[i];
     446             : 
     447             :             /* unrolled tail recursion */
     448      309314 :             for (a = 0; a < 1; a++) {
     449      154657 :                 silk_count_children(rc, 0, branch[0][a], branch[1]);
     450      463971 :                 for (b = 0; b < 2; b++) {
     451      309314 :                     silk_count_children(rc, 1, branch[1][b], branch[2]);
     452      927942 :                     for (c = 0; c < 2; c++) {
     453      618628 :                         silk_count_children(rc, 2, branch[2][c], branch[3]);
     454     1855884 :                         for (d = 0; d < 2; d++) {
     455     1237256 :                             silk_count_children(rc, 3, branch[3][d], location);
     456     1237256 :                             location += 2;
     457             :                         }
     458             :                     }
     459             :                 }
     460             :             }
     461             :         } else
     462       24865 :             memset(excitation + 16*i, 0, 16*sizeof(int32_t));
     463             :     }
     464             : 
     465             :     /* decode least significant bits */
     466     2885649 :     for (i = 0; i < shellblocks << 4; i++) {
     467             :         int bit;
     468     3083408 :         for (bit = 0; bit < lsbcount[i >> 4]; bit++)
     469      422112 :             excitation[i] = (excitation[i] << 1) |
     470      211056 :                             ff_opus_rc_dec_cdf(rc, ff_silk_model_excitation_lsb);
     471             :     }
     472             : 
     473             :     /* decode signs */
     474     2885649 :     for (i = 0; i < shellblocks << 4; i++) {
     475     2872352 :         if (excitation[i] != 0) {
     476     1329860 :             int sign = ff_opus_rc_dec_cdf(rc, ff_silk_model_excitation_sign[active +
     477      664930 :                                          voiced][qoffset_high][FFMIN(pulsecount[i >> 4], 6)]);
     478      664930 :             if (sign == 0)
     479      442065 :                 excitation[i] *= -1;
     480             :         }
     481             :     }
     482             : 
     483             :     /* assemble the excitation */
     484     2885649 :     for (i = 0; i < shellblocks << 4; i++) {
     485     2872352 :         int value = excitation[i];
     486     2872352 :         excitation[i] = value * 256 | ff_silk_quant_offset[voiced][qoffset_high];
     487     2872352 :         if (value < 0)      excitation[i] += 20;
     488     2430287 :         else if (value > 0) excitation[i] -= 20;
     489             : 
     490             :         /* invert samples pseudorandomly */
     491     2872352 :         seed = 196314165 * seed + 907633515;
     492     2872352 :         if (seed & 0x80000000)
     493     1421376 :             excitation[i] *= -1;
     494     2872352 :         seed += value;
     495             : 
     496     2872352 :         excitationf[i] = excitation[i] / 8388608.0f;
     497             :     }
     498       13297 : }
     499             : 
     500             : /** Maximum residual history according to 4.2.7.6.1 */
     501             : #define SILK_MAX_LAG  (288 + LTP_ORDER / 2)
     502             : 
     503             : /** Order of the LTP filter */
     504             : #define LTP_ORDER 5
     505             : 
     506       13297 : static void silk_decode_frame(SilkContext *s, OpusRangeCoder *rc,
     507             :                               int frame_num, int channel, int coded_channels, int active, int active1)
     508             : {
     509             :     /* per frame */
     510             :     int voiced;       // combines with active to indicate inactive, active, or active+voiced
     511             :     int qoffset_high;
     512             :     int order;                             // order of the LPC coefficients
     513             :     float lpc_leadin[16], lpc_body[16], residual[SILK_MAX_LAG + SILK_HISTORY];
     514             :     int has_lpc_leadin;
     515             :     float ltpscale;
     516             : 
     517             :     /* per subframe */
     518             :     struct {
     519             :         float gain;
     520             :         int pitchlag;
     521             :         float ltptaps[5];
     522             :     } sf[4];
     523             : 
     524       13297 :     SilkFrame * const frame = s->frame + channel;
     525             : 
     526             :     int i;
     527             : 
     528             :     /* obtain stereo weights */
     529       13297 :     if (coded_channels == 2 && channel == 0) {
     530             :         int n, wi[2], ws[2], w[2];
     531        4591 :         n     = ff_opus_rc_dec_cdf(rc, ff_silk_model_stereo_s1);
     532        4591 :         wi[0] = ff_opus_rc_dec_cdf(rc, ff_silk_model_stereo_s2) + 3 * (n / 5);
     533        4591 :         ws[0] = ff_opus_rc_dec_cdf(rc, ff_silk_model_stereo_s3);
     534        4591 :         wi[1] = ff_opus_rc_dec_cdf(rc, ff_silk_model_stereo_s2) + 3 * (n % 5);
     535        4591 :         ws[1] = ff_opus_rc_dec_cdf(rc, ff_silk_model_stereo_s3);
     536             : 
     537       13773 :         for (i = 0; i < 2; i++)
     538       18364 :             w[i] = ff_silk_stereo_weights[wi[i]] +
     539        9182 :                    (((ff_silk_stereo_weights[wi[i] + 1] - ff_silk_stereo_weights[wi[i]]) * 6554) >> 16)
     540        9182 :                     * (ws[i]*2 + 1);
     541             : 
     542        4591 :         s->stereo_weights[0] = (w[0] - w[1]) / 8192.0;
     543        4591 :         s->stereo_weights[1] = w[1]          / 8192.0;
     544             : 
     545             :         /* and read the mid-only flag */
     546        4591 :         s->midonly = active1 ? 0 : ff_opus_rc_dec_cdf(rc, ff_silk_model_mid_only);
     547             :     }
     548             : 
     549             :     /* obtain frame type */
     550       13297 :     if (!active) {
     551        3436 :         qoffset_high = ff_opus_rc_dec_cdf(rc, ff_silk_model_frame_type_inactive);
     552        3436 :         voiced = 0;
     553             :     } else {
     554        9861 :         int type = ff_opus_rc_dec_cdf(rc, ff_silk_model_frame_type_active);
     555        9861 :         qoffset_high = type & 1;
     556        9861 :         voiced = type >> 1;
     557             :     }
     558             : 
     559             :     /* obtain subframe quantization gains */
     560       55123 :     for (i = 0; i < s->subframes; i++) {
     561             :         int log_gain;     //Q7
     562             :         int ipart, fpart, lingain;
     563             : 
     564       53509 :         if (i == 0 && (frame_num == 0 || !frame->coded)) {
     565             :             /* gain is coded absolute */
     566       11683 :             int x = ff_opus_rc_dec_cdf(rc, ff_silk_model_gain_highbits[active + voiced]);
     567       11683 :             log_gain = (x<<3) | ff_opus_rc_dec_cdf(rc, ff_silk_model_gain_lowbits);
     568             : 
     569       11683 :             if (frame->coded)
     570       11642 :                 log_gain = FFMAX(log_gain, frame->log_gain - 16);
     571             :         } else {
     572             :             /* gain is coded relative */
     573       30143 :             int delta_gain = ff_opus_rc_dec_cdf(rc, ff_silk_model_gain_delta);
     574       30143 :             log_gain = av_clip_uintp2(FFMAX((delta_gain<<1) - 16,
     575             :                                      frame->log_gain + delta_gain - 4), 6);
     576             :         }
     577             : 
     578       41826 :         frame->log_gain = log_gain;
     579             : 
     580             :         /* approximate 2**(x/128) with a Q7 (i.e. non-integer) input */
     581       41826 :         log_gain = (log_gain * 0x1D1C71 >> 16) + 2090;
     582       41826 :         ipart = log_gain >> 7;
     583       41826 :         fpart = log_gain & 127;
     584       41826 :         lingain = (1 << ipart) + ((-174 * fpart * (128-fpart) >>16) + fpart) * ((1<<ipart) >> 7);
     585       41826 :         sf[i].gain = lingain / 65536.0f;
     586             :     }
     587             : 
     588             :     /* obtain LPC filter coefficients */
     589       13297 :     silk_decode_lpc(s, frame, rc, lpc_leadin, lpc_body, &order, &has_lpc_leadin, voiced);
     590             : 
     591             :     /* obtain pitch lags, if this is a voiced frame */
     592       13297 :     if (voiced) {
     593        4622 :         int lag_absolute = (!frame_num || !frame->prev_voiced);
     594             :         int primarylag;         // primary pitch lag for the entire SILK frame
     595             :         int ltpfilter;
     596             :         const int8_t * offsets;
     597             : 
     598        4622 :         if (!lag_absolute) {
     599         609 :             int delta = ff_opus_rc_dec_cdf(rc, ff_silk_model_pitch_delta);
     600         609 :             if (delta)
     601         538 :                 primarylag = frame->primarylag + delta - 9;
     602             :             else
     603          71 :                 lag_absolute = 1;
     604             :         }
     605             : 
     606        4622 :         if (lag_absolute) {
     607             :             /* primary lag is coded absolute */
     608             :             int highbits, lowbits;
     609             :             static const uint16_t * const model[] = {
     610             :                 ff_silk_model_pitch_lowbits_nb, ff_silk_model_pitch_lowbits_mb,
     611             :                 ff_silk_model_pitch_lowbits_wb
     612             :             };
     613        4084 :             highbits = ff_opus_rc_dec_cdf(rc, ff_silk_model_pitch_highbits);
     614        4084 :             lowbits  = ff_opus_rc_dec_cdf(rc, model[s->bandwidth]);
     615             : 
     616        8168 :             primarylag = ff_silk_pitch_min_lag[s->bandwidth] +
     617        4084 :                          highbits*ff_silk_pitch_scale[s->bandwidth] + lowbits;
     618             :         }
     619        4622 :         frame->primarylag = primarylag;
     620             : 
     621        4622 :         if (s->subframes == 2)
     622        3048 :             offsets = (s->bandwidth == OPUS_BANDWIDTH_NARROWBAND)
     623         239 :                      ? ff_silk_pitch_offset_nb10ms[ff_opus_rc_dec_cdf(rc,
     624             :                                                 ff_silk_model_pitch_contour_nb10ms)]
     625        1763 :                      : ff_silk_pitch_offset_mbwb10ms[ff_opus_rc_dec_cdf(rc,
     626             :                                                 ff_silk_model_pitch_contour_mbwb10ms)];
     627             :         else
     628        6196 :             offsets = (s->bandwidth == OPUS_BANDWIDTH_NARROWBAND)
     629         807 :                      ? ff_silk_pitch_offset_nb20ms[ff_opus_rc_dec_cdf(rc,
     630             :                                                 ff_silk_model_pitch_contour_nb20ms)]
     631        3905 :                      : ff_silk_pitch_offset_mbwb20ms[ff_opus_rc_dec_cdf(rc,
     632             :                                                 ff_silk_model_pitch_contour_mbwb20ms)];
     633             : 
     634       20062 :         for (i = 0; i < s->subframes; i++)
     635       30880 :             sf[i].pitchlag = av_clip(primarylag + offsets[i],
     636       15440 :                                      ff_silk_pitch_min_lag[s->bandwidth],
     637       15440 :                                      ff_silk_pitch_max_lag[s->bandwidth]);
     638             : 
     639             :         /* obtain LTP filter coefficients */
     640        4622 :         ltpfilter = ff_opus_rc_dec_cdf(rc, ff_silk_model_ltp_filter);
     641       20062 :         for (i = 0; i < s->subframes; i++) {
     642             :             int index, j;
     643             :             static const uint16_t * const filter_sel[] = {
     644             :                 ff_silk_model_ltp_filter0_sel, ff_silk_model_ltp_filter1_sel,
     645             :                 ff_silk_model_ltp_filter2_sel
     646             :             };
     647             :             static const int8_t (* const filter_taps[])[5] = {
     648             :                 ff_silk_ltp_filter0_taps, ff_silk_ltp_filter1_taps, ff_silk_ltp_filter2_taps
     649             :             };
     650       15440 :             index = ff_opus_rc_dec_cdf(rc, filter_sel[ltpfilter]);
     651       92640 :             for (j = 0; j < 5; j++)
     652       77200 :                 sf[i].ltptaps[j] = filter_taps[ltpfilter][index][j] / 128.0f;
     653             :         }
     654             :     }
     655             : 
     656             :     /* obtain LTP scale factor */
     657       13297 :     if (voiced && frame_num == 0)
     658        7852 :         ltpscale = ff_silk_ltp_scale_factor[ff_opus_rc_dec_cdf(rc,
     659        3926 :                                          ff_silk_model_ltp_scale_index)] / 16384.0f;
     660        9371 :     else ltpscale = 15565.0f/16384.0f;
     661             : 
     662             :     /* generate the excitation signal for the entire frame */
     663       13297 :     silk_decode_excitation(s, rc, residual + SILK_MAX_LAG, qoffset_high,
     664             :                            active, voiced);
     665             : 
     666             :     /* skip synthesising the side channel if we want mono-only */
     667       13297 :     if (s->output_channels == channel)
     668           0 :         return;
     669             : 
     670             :     /* generate the output signal */
     671       55123 :     for (i = 0; i < s->subframes; i++) {
     672       41826 :         const float * lpc_coeff = (i < 2 && has_lpc_leadin) ? lpc_leadin : lpc_body;
     673       41826 :         float *dst    = frame->output      + SILK_HISTORY + i * s->sflength;
     674       41826 :         float *resptr = residual           + SILK_MAX_LAG + i * s->sflength;
     675       41826 :         float *lpc    = frame->lpc_history + SILK_HISTORY + i * s->sflength;
     676             :         float sum;
     677             :         int j, k;
     678             : 
     679       41826 :         if (voiced) {
     680             :             int out_end;
     681             :             float scale;
     682             : 
     683       15440 :             if (i < 2 || s->nlsf_interp_factor == 4) {
     684       13554 :                 out_end = -i * s->sflength;
     685       13554 :                 scale   = ltpscale;
     686             :             } else {
     687        1886 :                 out_end = -(i - 2) * s->sflength;
     688        1886 :                 scale   = 1.0f;
     689             :             }
     690             : 
     691             :             /* when the LPC coefficients change, a re-whitening filter is used */
     692             :             /* to produce a residual that accounts for the change */
     693      666703 :             for (j = - sf[i].pitchlag - LTP_ORDER/2; j < out_end; j++) {
     694      651263 :                 sum = dst[j];
     695     9876529 :                 for (k = 0; k < order; k++)
     696     9225266 :                     sum -= lpc_coeff[k] * dst[j - k - 1];
     697      651263 :                 resptr[j] = av_clipf(sum, -1.0f, 1.0f) * scale / sf[i].gain;
     698             :             }
     699             : 
     700       15440 :             if (out_end) {
     701        9875 :                 float rescale = sf[i-1].gain / sf[i].gain;
     702     1094295 :                 for (j = out_end; j < 0; j++)
     703     1084420 :                     resptr[j] *= rescale;
     704             :             }
     705             : 
     706             :             /* LTP synthesis */
     707     1044440 :             for (j = 0; j < s->sflength; j++) {
     708     1029000 :                 sum = resptr[j];
     709     6174000 :                 for (k = 0; k < LTP_ORDER; k++)
     710     5145000 :                     sum += sf[i].ltptaps[k] * resptr[j - sf[i].pitchlag + LTP_ORDER/2 - k];
     711     1029000 :                 resptr[j] = sum;
     712             :             }
     713             :         }
     714             : 
     715             :         /* LPC synthesis */
     716     2909186 :         for (j = 0; j < s->sflength; j++) {
     717     2867360 :             sum = resptr[j] * sf[i].gain;
     718    44202400 :             for (k = 1; k <= order; k++)
     719    41335040 :                 sum += lpc_coeff[k - 1] * lpc[j - k];
     720             : 
     721     2867360 :             lpc[j] = sum;
     722     2867360 :             dst[j] = av_clipf(sum, -1.0f, 1.0f);
     723             :         }
     724             :     }
     725             : 
     726       13297 :     frame->prev_voiced = voiced;
     727       13297 :     memmove(frame->lpc_history, frame->lpc_history + s->flength, SILK_HISTORY * sizeof(float));
     728       13297 :     memmove(frame->output,      frame->output      + s->flength, SILK_HISTORY * sizeof(float));
     729             : 
     730       13297 :     frame->coded = 1;
     731             : }
     732             : 
     733        4591 : static void silk_unmix_ms(SilkContext *s, float *l, float *r)
     734             : {
     735        4591 :     float *mid    = s->frame[0].output + SILK_HISTORY - s->flength;
     736        4591 :     float *side   = s->frame[1].output + SILK_HISTORY - s->flength;
     737        4591 :     float w0_prev = s->prev_stereo_weights[0];
     738        4591 :     float w1_prev = s->prev_stereo_weights[1];
     739        4591 :     float w0      = s->stereo_weights[0];
     740        4591 :     float w1      = s->stereo_weights[1];
     741        4591 :     int n1        = ff_silk_stereo_interp_len[s->bandwidth];
     742             :     int i;
     743             : 
     744      521807 :     for (i = 0; i < n1; i++) {
     745      517216 :         float interp0 = w0_prev + i * (w0 - w0_prev) / n1;
     746      517216 :         float interp1 = w1_prev + i * (w1 - w1_prev) / n1;
     747      517216 :         float p0      = 0.25 * (mid[i - 2] + 2 * mid[i - 1] + mid[i]);
     748             : 
     749      517216 :         l[i] = av_clipf((1 + interp1) * mid[i - 1] + side[i - 1] + interp0 * p0, -1.0, 1.0);
     750      517216 :         r[i] = av_clipf((1 - interp1) * mid[i - 1] - side[i - 1] - interp0 * p0, -1.0, 1.0);
     751             :     }
     752             : 
     753      427255 :     for (; i < s->flength; i++) {
     754      422664 :         float p0 = 0.25 * (mid[i - 2] + 2 * mid[i - 1] + mid[i]);
     755             : 
     756      422664 :         l[i] = av_clipf((1 + w1) * mid[i - 1] + side[i - 1] + w0 * p0, -1.0, 1.0);
     757      422664 :         r[i] = av_clipf((1 - w1) * mid[i - 1] - side[i - 1] - w0 * p0, -1.0, 1.0);
     758             :     }
     759             : 
     760        4591 :     memcpy(s->prev_stereo_weights, s->stereo_weights, sizeof(s->stereo_weights));
     761        4591 : }
     762             : 
     763       49927 : static void silk_flush_frame(SilkFrame *frame)
     764             : {
     765       49927 :     if (!frame->coded)
     766       49894 :         return;
     767             : 
     768          33 :     memset(frame->output,      0, sizeof(frame->output));
     769          33 :     memset(frame->lpc_history, 0, sizeof(frame->lpc_history));
     770             : 
     771          33 :     memset(frame->lpc,  0, sizeof(frame->lpc));
     772          33 :     memset(frame->nlsf, 0, sizeof(frame->nlsf));
     773             : 
     774          33 :     frame->log_gain = 0;
     775             : 
     776          33 :     frame->primarylag  = 0;
     777          33 :     frame->prev_voiced = 0;
     778          33 :     frame->coded       = 0;
     779             : }
     780             : 
     781        9461 : int ff_silk_decode_superframe(SilkContext *s, OpusRangeCoder *rc,
     782             :                               float *output[2],
     783             :                               enum OpusBandwidth bandwidth,
     784             :                               int coded_channels,
     785             :                               int duration_ms)
     786             : {
     787             :     int active[2][6], redundancy[2];
     788             :     int nb_frames, i, j;
     789             : 
     790        9461 :     if (bandwidth > OPUS_BANDWIDTH_WIDEBAND ||
     791        9461 :         coded_channels > 2 || duration_ms > 60) {
     792           0 :         av_log(s->avctx, AV_LOG_ERROR, "Invalid parameters passed "
     793             :                "to the SILK decoder.\n");
     794           0 :         return AVERROR(EINVAL);
     795             :     }
     796             : 
     797        9461 :     nb_frames = 1 + (duration_ms > 20) + (duration_ms > 40);
     798        9461 :     s->subframes = duration_ms / nb_frames / 5;         // 5ms subframes
     799        9461 :     s->sflength  = 20 * (bandwidth + 2);
     800        9461 :     s->flength   = s->sflength * s->subframes;
     801        9461 :     s->bandwidth = bandwidth;
     802        9461 :     s->wb        = bandwidth == OPUS_BANDWIDTH_WIDEBAND;
     803             : 
     804             :     /* make sure to flush the side channel when switching from mono to stereo */
     805        9461 :     if (coded_channels > s->prev_coded_channels)
     806          24 :         silk_flush_frame(&s->frame[1]);
     807        9461 :     s->prev_coded_channels = coded_channels;
     808             : 
     809             :     /* read the LP-layer header bits */
     810       22975 :     for (i = 0; i < coded_channels; i++) {
     811       28642 :         for (j = 0; j < nb_frames; j++)
     812       15128 :             active[i][j] = ff_opus_rc_dec_log(rc, 1);
     813             : 
     814       13514 :         redundancy[i] = ff_opus_rc_dec_log(rc, 1);
     815       13514 :         if (redundancy[i]) {
     816           0 :             avpriv_report_missing_feature(s->avctx, "LBRR frames");
     817           0 :             return AVERROR_PATCHWELCOME;
     818             :         }
     819             :     }
     820             : 
     821       19998 :     for (i = 0; i < nb_frames; i++) {
     822       23834 :         for (j = 0; j < coded_channels && !s->midonly; j++)
     823       13297 :             silk_decode_frame(s, rc, i, j, coded_channels, active[j][i], active[1][i]);
     824             : 
     825             :         /* reset the side channel if it is not coded */
     826       10537 :         if (s->midonly && s->frame[1].coded)
     827          13 :             silk_flush_frame(&s->frame[1]);
     828             : 
     829       16483 :         if (coded_channels == 1 || s->output_channels == 1) {
     830       17838 :             for (j = 0; j < s->output_channels; j++) {
     831       23784 :                 memcpy(output[j] + i * s->flength,
     832       11892 :                        s->frame[0].output + SILK_HISTORY - s->flength - 2,
     833       11892 :                        s->flength * sizeof(float));
     834             :             }
     835             :         } else {
     836        4591 :             silk_unmix_ms(s, output[0] + i * s->flength, output[1] + i * s->flength);
     837             :         }
     838             : 
     839       10537 :         s->midonly        = 0;
     840             :     }
     841             : 
     842        9461 :     return nb_frames * s->flength;
     843             : }
     844             : 
     845          41 : void ff_silk_free(SilkContext **ps)
     846             : {
     847          41 :     av_freep(ps);
     848          41 : }
     849             : 
     850       24945 : void ff_silk_flush(SilkContext *s)
     851             : {
     852       24945 :     silk_flush_frame(&s->frame[0]);
     853       24945 :     silk_flush_frame(&s->frame[1]);
     854             : 
     855       24945 :     memset(s->prev_stereo_weights, 0, sizeof(s->prev_stereo_weights));
     856       24945 : }
     857             : 
     858          41 : int ff_silk_init(AVCodecContext *avctx, SilkContext **ps, int output_channels)
     859             : {
     860             :     SilkContext *s;
     861             : 
     862          41 :     if (output_channels != 1 && output_channels != 2) {
     863           0 :         av_log(avctx, AV_LOG_ERROR, "Invalid number of output channels: %d\n",
     864             :                output_channels);
     865           0 :         return AVERROR(EINVAL);
     866             :     }
     867             : 
     868          41 :     s = av_mallocz(sizeof(*s));
     869          41 :     if (!s)
     870           0 :         return AVERROR(ENOMEM);
     871             : 
     872          41 :     s->avctx           = avctx;
     873          41 :     s->output_channels = output_channels;
     874             : 
     875          41 :     ff_silk_flush(s);
     876             : 
     877          41 :     *ps = s;
     878             : 
     879          41 :     return 0;
     880             : }

Generated by: LCOV version 1.13