lpc.c 8.24 KB
Newer Older
1
/*
2
 * LPC utility code
3
 * Copyright (c) 2006  Justin Ruggles <justin.ruggles@gmail.com>
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
 *
 * 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
 */

22
#include "libavutil/common.h"
23
#include "libavutil/lls.h"
24

25 26
#define LPC_USE_DOUBLE
#include "lpc.h"
Michael Niedermayer's avatar
Michael Niedermayer committed
27
#include "libavutil/avassert.h"
28 29


30 31 32
/**
 * Apply Welch window function to audio block
 */
33 34
static void lpc_apply_welch_window_c(const int32_t *data, int len,
                                     double *w_data)
35 36 37 38 39
{
    int i, n2;
    double w;
    double c;

40 41
    /* The optimization in commit fa4ed8c does not support odd len.
     * If someone wants odd len extend that change. */
Michael Niedermayer's avatar
Michael Niedermayer committed
42
    av_assert2(!(len & 1));
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

    n2 = (len >> 1);
    c = 2.0 / (len - 1.0);

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

/**
58
 * Calculate autocorrelation data from audio samples
59 60
 * A Welch window function is applied before calculation.
 */
61
static void lpc_compute_autocorr_c(const double *data, int len, int lag,
62
                                   double *autoc)
63 64 65 66 67 68
{
    int i, j;

    for(j=0; j<lag; j+=2){
        double sum0 = 1.0, sum1 = 1.0;
        for(i=j; i<len; i++){
69 70
            sum0 += data[i] * data[i-j];
            sum1 += data[i] * data[i-j-1];
71 72 73 74 75 76 77 78
        }
        autoc[j  ] = sum0;
        autoc[j+1] = sum1;
    }

    if(j==lag){
        double sum = 1.0;
        for(i=j-1; i<len; i+=2){
79 80
            sum += data[i  ] * data[i-j  ]
                 + data[i+1] * data[i-j+1];
81 82 83 84 85
        }
        autoc[j] = sum;
    }
}

86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
/**
 * Quantize LPC coefficients
 */
static void quantize_lpc_coefs(double *lpc_in, int order, int precision,
                               int32_t *lpc_out, int *shift, int max_shift, int zero_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_shift) < 1.0) {
        *shift = zero_shift;
        memset(lpc_out, 0, sizeof(int32_t) * order);
        return;
    }

    /* calculate level shift which scales max coeff to available bits */
    sh = max_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++) {
131
        error -= lpc_in[i] * (1 << sh);
132 133 134 135 136 137
        lpc_out[i] = av_clip(lrintf(error), -qmax, qmax);
        error -= lpc_out[i];
    }
    *shift = sh;
}

138
static int estimate_best_order(double *ref, int min_order, int max_order)
139 140 141
{
    int i, est;

142 143
    est = min_order;
    for(i=max_order-1; i>=min_order-1; i--) {
144 145 146 147 148 149 150 151
        if(ref[i] > 0.10) {
            est = i+1;
            break;
        }
    }
    return est;
}

152 153 154 155 156 157 158 159 160 161 162 163
int ff_lpc_calc_ref_coefs(LPCContext *s,
                          const int32_t *samples, int order, double *ref)
{
    double autoc[MAX_LPC_ORDER + 1];

    s->lpc_apply_welch_window(samples, s->blocksize, s->windowed_samples);
    s->lpc_compute_autocorr(s->windowed_samples, s->blocksize, order, autoc);
    compute_ref_coefs(autoc, order, ref, NULL);

    return order;
}

164 165
/**
 * Calculate LPC coefficients for multiple orders
166
 *
167 168
 * @param lpc_type LPC method for determining coefficients,
 *                 see #FFLPCType for details
169
 */
170
int ff_lpc_calc_coefs(LPCContext *s,
171 172
                      const int32_t *samples, int blocksize, int min_order,
                      int max_order, int precision,
173
                      int32_t coefs[][MAX_LPC_ORDER], int *shift,
174
                      enum FFLPCType lpc_type, int lpc_passes,
175
                      int omethod, int max_shift, int zero_shift)
176 177 178 179 180 181 182
{
    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;

Michael Niedermayer's avatar
Michael Niedermayer committed
183
    av_assert2(max_order >= MIN_LPC_ORDER && max_order <= MAX_LPC_ORDER &&
184
           lpc_type > FF_LPC_TYPE_FIXED);
185

186 187 188 189 190 191 192
    /* reinit LPC context if parameters have changed */
    if (blocksize != s->blocksize || max_order != s->max_order ||
        lpc_type  != s->lpc_type) {
        ff_lpc_end(s);
        ff_lpc_init(s, blocksize, max_order, lpc_type);
    }

193
    if (lpc_type == FF_LPC_TYPE_LEVINSON) {
194
        s->lpc_apply_welch_window(samples, blocksize, s->windowed_samples);
195

196
        s->lpc_compute_autocorr(s->windowed_samples, blocksize, max_order, autoc);
197

198 199 200 201
        compute_lpc_coefs(autoc, max_order, &lpc[0][0], MAX_LPC_ORDER, 0, 1);

        for(i=0; i<max_order; i++)
            ref[i] = fabs(lpc[i][i]);
202
    } else if (lpc_type == FF_LPC_TYPE_CHOLESKY) {
203
        LLSModel m[2];
204
        double var[MAX_LPC_ORDER+1], av_uninit(weight);
205

206 207 208
        if(lpc_passes <= 0)
            lpc_passes = 2;

209
        for(pass=0; pass<lpc_passes; pass++){
210
            avpriv_init_lls(&m[pass&1], max_order);
211 212 213 214 215 216 217 218

            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;
219
                    eval= avpriv_evaluate_lls(&m[(pass-1)&1], var+1, max_order-1);
220 221 222 223 224 225 226 227 228
                    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++;

229
                avpriv_update_lls(&m[pass&1], var, 1.0);
230
            }
231
            avpriv_solve_lls(&m[pass&1], 0.001, 0);
232 233 234 235
        }

        for(i=0; i<max_order; i++){
            for(j=0; j<max_order; j++)
236
                lpc[i][j]=-m[(pass-1)&1].coeff[i][j];
237 238 239 240
            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];
241 242
    } else
        av_assert0(0);
243 244 245
    opt_order = max_order;

    if(omethod == ORDER_METHOD_EST) {
246
        opt_order = estimate_best_order(ref, min_order, max_order);
247 248 249
        i = opt_order-1;
        quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
    } else {
250
        for(i=min_order-1; i<max_order; i++) {
251 252 253 254 255 256
            quantize_lpc_coefs(lpc[i], i+1, precision, coefs[i], &shift[i], max_shift, zero_shift);
        }
    }

    return opt_order;
}
257

258
av_cold int ff_lpc_init(LPCContext *s, int blocksize, int max_order,
259
                        enum FFLPCType lpc_type)
260
{
261 262 263 264
    s->blocksize = blocksize;
    s->max_order = max_order;
    s->lpc_type  = lpc_type;

265
    if (lpc_type == FF_LPC_TYPE_LEVINSON) {
266 267 268
        s->windowed_buffer = av_mallocz((blocksize + 2 + FFALIGN(max_order, 4)) *
                                        sizeof(*s->windowed_samples));
        if (!s->windowed_buffer)
269
            return AVERROR(ENOMEM);
270
        s->windowed_samples = s->windowed_buffer + FFALIGN(max_order, 4);
271 272 273 274
    } else {
        s->windowed_samples = NULL;
    }

275 276
    s->lpc_apply_welch_window = lpc_apply_welch_window_c;
    s->lpc_compute_autocorr   = lpc_compute_autocorr_c;
277

278
    if (ARCH_X86)
279
        ff_lpc_init_x86(s);
280 281 282 283 284 285

    return 0;
}

av_cold void ff_lpc_end(LPCContext *s)
{
286
    av_freep(&s->windowed_buffer);
287
}