flacenc.c

ffmpeg的完整源代码和作者自己写的文档。不但有在Linux的工程哦

    res = &data[pred_order];
    res_end = &data[n >> pmax];
    for(i=0; i<parts; i++) {
        uint32_t sum = 0;
        while(res < res_end){
            sum += *(res++);
        }
        sums[pmax][i] = sum;
        res_end+= n >> pmax;
    }
    /* sums for lower levels */
    for(i=pmax-1; i>=pmin; i--) {
        parts = (1 << i);
        for(j=0; j<parts; j++) {
            sums[i][j] = sums[i+1][2*j] + sums[i+1][2*j+1];
        }
    }
}

static uint32_t calc_rice_params(RiceContext *rc, int pmin, int pmax,
                                 int32_t *data, int n, int pred_order)
{
    int i;
    uint32_t bits[MAX_PARTITION_ORDER+1];
    int opt_porder;
    RiceContext tmp_rc;
    uint32_t *udata;
    uint32_t sums[MAX_PARTITION_ORDER+1][MAX_PARTITIONS];

    assert(pmin >= 0 && pmin <= MAX_PARTITION_ORDER);
    assert(pmax >= 0 && pmax <= MAX_PARTITION_ORDER);
    assert(pmin <= pmax);

    udata = av_malloc(n * sizeof(uint32_t));
    for(i=0; i<n; i++) {
        udata[i] = (2*data[i]) ^ (data[i]>>31);
    }

    calc_sums(pmin, pmax, udata, n, pred_order, sums);

    opt_porder = pmin;
    bits[pmin] = UINT32_MAX;
    for(i=pmin; i<=pmax; i++) {
        bits[i] = calc_optimal_rice_params(&tmp_rc, i, sums[i], n, pred_order);
        if(bits[i] <= bits[opt_porder]) {
            opt_porder = i;
            *rc= tmp_rc;
        }
    }

    av_freep(&udata);
    return bits[opt_porder];
}

static int get_max_p_order(int max_porder, int n, int order)
{
    int porder = FFMIN(max_porder, av_log2(n^(n-1)));
    if(order > 0)
        porder = FFMIN(porder, av_log2(n/order));
    return porder;
}

static uint32_t calc_rice_params_fixed(RiceContext *rc, int pmin, int pmax,
                                       int32_t *data, int n, int pred_order,
                                       int bps)
{
    uint32_t bits;
    pmin = get_max_p_order(pmin, n, pred_order);
    pmax = get_max_p_order(pmax, n, pred_order);
    bits = pred_order*bps + 6;
    bits += calc_rice_params(rc, pmin, pmax, data, n, pred_order);
    return bits;
}

static uint32_t calc_rice_params_lpc(RiceContext *rc, int pmin, int pmax,
                                     int32_t *data, int n, int pred_order,
                                     int bps, int precision)
{
    uint32_t bits;
    pmin = get_max_p_order(pmin, n, pred_order);
    pmax = get_max_p_order(pmax, n, pred_order);
    bits = pred_order*bps + 4 + 5 + pred_order*precision + 6;
    bits += calc_rice_params(rc, pmin, pmax, data, n, pred_order);
    return bits;
}

/**
 * Apply Welch window function to audio block
 */
static void apply_welch_window(const int32_t *data, int len, double *w_data)
{
    int i, n2;
    double w;
    double c;

    n2 = (len >> 1);
    c = 2.0 / (len - 1.0);
    for(i=0; i<n2; i++) {
        w = c - i - 1.0;
        w = 1.0 - (w * w);
        w_data[i] = data[i] * w;
        w_data[len-1-i] = data[len-1-i] * w;
    }
}

/**
 * Calculates autocorrelation data from audio samples
 * A Welch window function is applied before calculation.
 */
void ff_flac_compute_autocorr(const int32_t *data, int len, int lag,
                              double *autoc)
{
    int i, j;
    double tmp[len + lag + 1];
    double *data1= tmp + lag;

    apply_welch_window(data, len, data1);

    for(j=0; j<lag; j++)
        data1[j-lag]= 0.0;
    data1[len] = 0.0;

    for(j=0; j<lag; j+=2){
        double sum0 = 1.0, sum1 = 1.0;
        for(i=0; i<len; i++){
            sum0 += data1[i] * data1[i-j];
            sum1 += data1[i] * data1[i-j-1];
        }
        autoc[j  ] = sum0;
        autoc[j+1] = sum1;
    }

    if(j==lag){
        double sum = 1.0;
        for(i=0; i<len; i+=2){
            sum += data1[i  ] * data1[i-j  ]
                 + data1[i+1] * data1[i-j+1];
        }
        autoc[j] = sum;
    }
}

/**
 * Levinson-Durbin recursion.
 * Produces LPC coefficients from autocorrelation data.
 */
static void compute_lpc_coefs(const double *autoc, int max_order,
                              double lpc[][MAX_LPC_ORDER], double *ref)
{
   int i, j, i2;
   double r, err, tmp;
   double lpc_tmp[MAX_LPC_ORDER];

   for(i=0; i<max_order; i++) lpc_tmp[i] = 0;
   err = autoc[0];

   for(i=0; i<max_order; i++) {
      r = -autoc[i+1];
      for(j=0; j<i; j++) {
          r -= lpc_tmp[j] * autoc[i-j];
      }
      r /= err;
      ref[i] = fabs(r);

      err *= 1.0 - (r * r);

      i2 = (i >> 1);
      lpc_tmp[i] = r;
      for(j=0; j<i2; j++) {
         tmp = lpc_tmp[j];
         lpc_tmp[j] += r * lpc_tmp[i-1-j];
         lpc_tmp[i-1-j] += r * tmp;
      }
      if(i & 1) {
          lpc_tmp[j] += lpc_tmp[j] * r;
      }

      for(j=0; j<=i; j++) {
          lpc[i][j] = -lpc_tmp[j];
      }
   }
}

/**
 * Quantize LPC coefficients
 */
static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
                               int32_t *lpc_out, int *shift)
{
    int i;
    double cmax, error;
    int32_t qmax;
    int sh;

    /* define maximum levels */
    qmax = (1 << (precision - 1)) - 1;

    /* find maximum coefficient value */
    cmax = 0.0;
    for(i=0; i<order; i++) {
        cmax= FFMAX(cmax, fabs(lpc_in[i]));
    }

    /* if maximum value quantizes to zero, return all zeros */
    if(cmax * (1 << MAX_LPC_SHIFT) < 1.0) {
        *shift = 0;
        memset(lpc_out, 0, sizeof(int32_t) * order);
        return;
    }

    /* calculate level shift which scales max coeff to available bits */
    sh = MAX_LPC_SHIFT;
    while((cmax * (1 << sh) > qmax) && (sh > 0)) {
        sh--;
    }

    /* since negative shift values are unsupported in decoder, scale down
       coefficients instead */
    if(sh == 0 && cmax > qmax) {
        double scale = ((double)qmax) / cmax;
        for(i=0; i<order; i++) {
            lpc_in[i] *= scale;
        }
    }

    /* output quantized coefficients and level shift */
    error=0;
    for(i=0; i<order; i++) {
        error += lpc_in[i] * (1 << sh);
        lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
        error -= lpc_out[i];
    }
    *shift = sh;
}

static int estimate_best_order(double *ref, int max_order)
{
    int i, est;

    est = 1;
    for(i=max_order-1; i>=0; i--) {
        if(ref[i] > 0.10) {
            est = i+1;
            break;
        }
    }
    return est;
}

/**
 * Calculate LPC coefficients for multiple orders
 */
static int lpc_calc_coefs(FlacEncodeContext *s,
                          const int32_t *samples, int blocksize, int max_order,
                          int precision, int32_t coefs[][MAX_LPC_ORDER],
                          int *shift, int use_lpc, int omethod)
{
    double autoc[MAX_LPC_ORDER+1];
    double ref[MAX_LPC_ORDER];
    double lpc[MAX_LPC_ORDER][MAX_LPC_ORDER];
    int i, j, pass;
    int opt_order;

    assert(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER);

    if(use_lpc == 1){
        s->dsp.flac_compute_autocorr(samples, blocksize, max_order, autoc);

        compute_lpc_coefs(autoc, max_order, lpc, ref);
    }else{
        LLSModel m[2];
        double var[MAX_LPC_ORDER+1], weight;

        for(pass=0; pass<use_lpc-1; pass++){
            av_init_lls(&m[pass&1], max_order);

            weight=0;
            for(i=max_order; i<blocksize; i++){
                for(j=0; j<=max_order; j++)
                    var[j]= samples[i-j];

                if(pass){
                    double eval, inv, rinv;
                    eval= av_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
                    eval= (512>>pass) + fabs(eval - var[0]);
                    inv = 1/eval;
                    rinv = sqrt(inv);
                    for(j=0; j<=max_order; j++)
                        var[j] *= rinv;
                    weight += inv;
                }else
                    weight++;

                av_update_lls(&m[pass&1], var, 1.0);
            }
            av_solve_lls(&m[pass&1], 0.001, 0);
        }

        for(i=0; i<max_order; i++){
            for(j=0; j<max_order; j++)
                lpc[i][j]= m[(pass-1)&1].coeff[i][j];
            ref[i]= sqrt(m[(pass-1)&1].variance[i] / weight) * (blocksize - max_order) / 4000;
        }
        for(i=max_order-1; i>0; i--)
            ref[i] = ref[i-1] - ref[i];
    }
    opt_order = max_order;

    if(omethod == ORDER_METHOD_EST) {
        opt_order = estimate_best_order(ref, max_order);
        i = opt_order-1;
        quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
    } else {
        for(i=0; i<max_order; i++) {
            quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i]);
        }
    }

    return opt_order;
}


static void encode_residual_verbatim(int32_t *res, int32_t *smp, int n)
{
    assert(n > 0);
    memcpy(res, smp, n * sizeof(int32_t));
}

static void encode_residual_fixed(int32_t *res, const int32_t *smp, int n,
                                  int order)
{
    int i;

    for(i=0; i<order; i++) {
        res[i] = smp[i];
    }

    if(order==0){
        for(i=order; i<n; i++)
            res[i]= smp[i];
    }else if(order==1){
        for(i=order; i<n; i++)
            res[i]= smp[i] - smp[i-1];
    }else if(order==2){
        int a = smp[order-1] - smp[order-2];
        for(i=order; i<n; i+=2) {
            int b = smp[i] - smp[i-1];
            res[i]= b - a;
            a = smp[i+1] - smp[i];
            res[i+1]= a - b;
        }
    }else if(order==3){
        int a = smp[order-1] - smp[order-2];
        int c = smp[order-1] - 2*smp[order-2] + smp[order-3];
        for(i=order; i<n; i+=2) {
            int b = smp[i] - smp[i-1];
            int d = b - a;
            res[i]= d - c;
            a = smp[i+1] - smp[i];
            c = a - b;
            res[i+1]= c - d;
        }
    }else{
        int a = smp[order-1] - smp[order-2];
        int c = smp[order-1] - 2*smp[order-2] + smp[order-3];
        int e = smp[order-1] - 3*smp[order-2] + 3*smp[order-3] - smp[order-4];
        for(i=order; i<n; i+=2) {
            int b = smp[i] - smp[i-1];
            int d = b - a;
            int f = d - c;
            res[i]= f - e;
            a = smp[i+1] - smp[i];
            c = a - b;
            e = c - d;
            res[i+1]= e - f;
        }
    }
}

#define LPC1(x) {\
    int c = coefs[(x)-1];\
    p0 += c*s;\
    s = smp[i-(x)+1];\
    p1 += c*s;\
}

static av_always_inline void encode_residual_lpc_unrolled(
    int32_t *res, const int32_t *smp, int n,
    int order, const int32_t *coefs, int shift, int big)
{
    int i;
    for(i=order; i<n; i+=2) {
        int s = smp[i-order];
        int p0 = 0, p1 = 0;
        if(big) {
            switch(order) {
                case 32: LPC1(32)
                case 31: LPC1(31)
                case 30: LPC1(30)
                case 29: LPC1(29)
                case 28: LPC1(28)
                case 27: LPC1(27)
                case 26: LPC1(26)
                case 25: LPC1(25)
                case 24: LPC1(24)
                case 23: LPC1(23)
                case 22: LPC1(22)
                case 21: LPC1(21)
                case 20: LPC1(20)
                case 19: LPC1(19)
                case 18: LPC1(18)
                case 17: LPC1(17)
                case 16: LPC1(16)
                case 15: LPC1(15)
                case 14: LPC1(14)
                case 13: LPC1(13)
                case 12: LPC1(12)
                case 11: LPC1(11)
                case 10: LPC1(10)
                case  9: LPC1( 9)
                         LPC1( 8)
                         LPC1( 7)
                         LPC1( 6)
                         LPC1( 5)
                         LPC1( 4)
                         LPC1( 3)
                         LPC1( 2)
                         LPC1( 1)
            }
        } else {
            switch(order) {
                case  8: LPC1( 8)
                case  7: LPC1( 7)
                case  6: LPC1( 6)
                case  5: LPC1( 5)
                case  4: LPC1( 4)
                case  3: LPC1( 3)
                case  2: LPC1( 2)
                case  1: LPC1( 1)
            }
        }
        res[i  ] = smp[i  ] - (p0 >> shift);
        res[i+1] = smp[i+1] - (p1 >> shift);
    }
}

static void encode_residual_lpc(int32_t *res, const int32_t *smp, int n,
                                int order, const int32_t *coefs, int shift)
{
    int i;
    for(i=0; i<order; i++) {
        res[i] = smp[i];
    }
#ifdef CONFIG_SMALL
    for(i=order; i<n; i+=2) {
        int j;
        int s = smp[i];
        int p0 = 0, p1 = 0;
        for(j=0; j<order; j++) {
            int c = coefs[j];
            p1 += c*s;
            s = smp[i-j-1];
            p0 += c*s;
        }
        res[i  ] = smp[i  ] - (p0 >> shift);
        res[i+1] = smp[i+1] - (p1 >> shift);
    }
#else
    switch(order) {
        case  1: encode_residual_lpc_unrolled(res, smp, n, 1, coefs, shift, 0); break;
        case  2: encode_residual_lpc_unrolled(res, smp, n, 2, coefs, shift, 0); break;
        case  3: encode_residual_lpc_unrolled(res, smp, n, 3, coefs, shift, 0); break;
        case  4: encode_residual_lpc_unrolled(res, smp, n, 4, coefs, shift, 0); break;
        case  5: encode_residual_lpc_unrolled(res, smp, n, 5, coefs, shift, 0); break;
        case  6: encode_residual_lpc_unrolled(res, smp, n, 6, coefs, shift, 0); break;
        case  7: encode_residual_lpc_unrolled(res, smp, n, 7, coefs, shift, 0); break;
        case  8: encode_residual_lpc_unrolled(res, smp, n, 8, coefs, shift, 0); break;
        default: encode_residual_lpc_unrolled(res, smp, n, order, coefs, shift, 1); break;
    }
#endif
}

static int encode_residual(FlacEncodeContext *ctx, int ch)
{
    int i, n;
    int min_order, max_order, opt_order, precision, omethod;
    int min_porder, max_porder;
    FlacFrame *frame;
    FlacSubframe *sub;
    int32_t coefs[MAX_LPC_ORDER][MAX_LPC_ORDER];
    int shift[MAX_LPC_ORDER];
    int32_t *res, *smp;

    frame = &ctx->frame;
    sub = &frame->subframes[ch];
    res = sub->residual;
    smp = sub->samples;
    n = frame->blocksize;

    /* CONSTANT */
    for(i=1; i<n; i++) {