comparison mpegaudiodec.c @ 449:7b2f23b5dcdd libavcodec

fixed layer1/2 overflow if very loud sound - fixed broken free format decoding to pass all mpeg audio standard decoding tests (please avoid patching the parser without having all test streams available - contact me if necessary)

448:e8c8ca9106aa

 *  - test lsf / mpeg25 extensively.
 */

/* define USE_HIGHPRECISION to have a bit exact (but slower) mpeg
   audio decoder */
//#define USE_HIGHPRECISION

#ifdef USE_HIGHPRECISION
#define FRAC_BITS   23   /* fractional bits for sb_samples and dct */
#define WFRAC_BITS  16   /* fractional bits for window */
#else

#define SCALE_GEN(v) \
{ FIXR(1.0 * (v)), FIXR(0.7937005259 * (v)), FIXR(0.6299605249 * (v)) }

static INT32 scale_factor_mult2[3][3] = {
    SCALE_GEN(1.0 / 3.0), /* 3 steps */
    SCALE_GEN(1.0 / 5.0), /* 5 steps */
    SCALE_GEN(1.0 / 9.0), /* 9 steps */
};

/* 2^(n/4) */
static UINT32 scale_factor_mult3[4] = {
    FIXR(1.0),

    shift = scale_factor_modshift[scale_factor];
    mod = shift & 3;
    shift >>= 2;
    val = MUL64(mant + (-1 << n) + 1, scale_factor_mult[n-1][mod]);
    shift += n;
    return (int)((val + (1 << (shift - 1))) >> shift);
}

static inline int l2_unscale_group(int steps, int mant, int scale_factor)
{
    int shift, mod, val;

    shift = scale_factor_modshift[scale_factor];
    mod = shift & 3;
    shift >>= 2;
    /* XXX: store the result directly */
    val = (2 * (mant - (steps >> 1))) * scale_factor_mult2[steps >> 2][mod];
    return (val + (1 << (shift - 1))) >> shift;
}

/* compute value^(4/3) * 2^(exponent/4). It normalized to FRAC_BITS */
static inline int l3_unscale(int value, int exponent)
{

    while (a < POW_FRAC_ONE) {
        a = a << 1;
        eq--;
    }
    /* now POW_FRAC_ONE <= a < 2 * POW_FRAC_ONE */
#if (POW_FRAC_BITS - 1) > FRAC_BITS
    a = (a + (1 << (POW_FRAC_BITS - FRAC_BITS - 1))) >> (POW_FRAC_BITS - FRAC_BITS);
    /* correct overflow */
    if (a >= 2 * (1 << FRAC_BITS)) {
        a = a >> 1;
        eq++;

    if(!init) {
        /* scale factors table for layer 1/2 */
        for(i=0;i<64;i++) {
            int shift, mod;
            /* 1.0 (i = 3) is normalized to 2 ^ FRAC_BITS */
            shift = (i / 3) - 1;
            mod = i % 3;
#if FRAC_BITS <= 15
            if (shift > 31)
                shift = 31;
#endif
            scale_factor_modshift[i] = mod | (shift << 2);
        }

        /* scale factor multiply for layer 1 */
        for(i=0;i<15;i++) {
            int n, norm;
            n = i + 2;
            norm = ((INT64_C(1) << n) * FRAC_ONE) / ((1 << n) - 1);
            scale_factor_mult[i][0] = MULL(FIXR(1.0), norm);
            scale_factor_mult[i][1] = MULL(FIXR(0.7937005259), norm);
            scale_factor_mult[i][2] = MULL(FIXR(0.6299605249), norm);
            dprintf("%d: norm=%x s=%x %x %x\n",
                    i, norm, 
                    scale_factor_mult[i][0],
                    scale_factor_mult[i][1],
                    scale_factor_mult[i][2]);

    synth_buf = s1->synth_buf[ch] + offset;

    for(j=0;j<32;j++) {
        v = tmp[j];
#if FRAC_BITS <= 15
        if (v > 32767)
            v = 32767;
        else if (v < -32768)
            v = -32768;
#endif

    
    s->layer = 4 - ((header >> 17) & 3);
    /* extract frequency */
    sample_rate_index = (header >> 10) & 3;
    sample_rate = mpa_freq_tab[sample_rate_index] >> (s->lsf + mpeg25);
    if (sample_rate == 0)
        return 1;
    sample_rate_index += 3 * (s->lsf + mpeg25);
    s->sample_rate_index = sample_rate_index;
    s->error_protection = ((header >> 16) & 1) ^ 1;

    bitrate_index = (header >> 12) & 0xf;
    padding = (header >> 9) & 1;
    //extension = (header >> 8) & 1;
    s->mode = (header >> 6) & 3;

            s->frame_size += padding;
            s->bit_rate = (s->frame_size * (sample_rate << s->lsf)) / 144000;
            break;
        }
    }
    s->sample_rate = sample_rate;
    
#if defined(DEBUG)
    printf("layer%d, %d Hz, %d kbits/s, ",
           s->layer, s->sample_rate, s->bit_rate);
    if (s->nb_channels == 2) {

#if defined(DEBUG)
void sample_dump(int fnum, INT32 *tab, int n)
{
    static FILE *files[16], *f;
    char buf[512];

    f = files[fnum];
    if (!f) {
        sprintf(buf, "/tmp/out%d.pcm", fnum);
        f = fopen(buf, "w");
        if (!f)
            return;
        files[fnum] = f;
    }
    
    if (fnum == 0) {
        int i;
        static int pos = 0;
        printf("pos=%d\n", pos);
        for(i=0;i<n;i++) {
            printf(" %f", (double)tab[i] / 32768.0);
            if ((i % 18) == 17)
                printf("\n");
        }
        pos += n;
    }

    fwrite(tab, 1, n * sizeof(INT32), f);
}
#endif


/* main layer3 decoding function */

            reorder_block(s, g);
#if defined(DEBUG)
            sample_dump(0, g->sb_hybrid, 576);
#endif
            compute_antialias(s, g);
#ifdef DEBUG
            sample_dump(1, g->sb_hybrid, 576);
#endif
            compute_imdct(s, g, &s->sb_samples[ch][18 * gr][0], s->mdct_buf[ch]); 
#ifdef DEBUG
            sample_dump(2, &s->sb_samples[ch][18 * gr][0], 576);
#endif
        }
    } /* gr */
    return nb_granules * 18;

                    /* reset free format frame size to give a chance
                       to get a new bitrate */
                    s->free_format_frame_size = 0;
		} else {
		    if (decode_header(s, header) == 1) {
                        /* free format: compute frame size */
			s->frame_size = -1;
			memcpy(s->inbuf, s->inbuf + 1, s->inbuf_ptr - s->inbuf - 1);
			s->inbuf_ptr--;
                    } else {
                        /* update codec info */
                        avctx->sample_rate = s->sample_rate;
                        avctx->channels = s->nb_channels;
			avctx->bit_rate = s->bit_rate;
			avctx->frame_size = s->frame_size;
                    }
		}
	    }
        } else if (s->frame_size == -1) {
            /* free format : find next sync to compute frame size */
	    len = MPA_MAX_CODED_FRAME_SIZE - len;

449:7b2f23b5dcdd

 *  - test lsf / mpeg25 extensively.
 */

/* define USE_HIGHPRECISION to have a bit exact (but slower) mpeg
   audio decoder */
#ifdef CONFIG_MPEGAUDIO_HP
#define USE_HIGHPRECISION
#endif

#ifdef USE_HIGHPRECISION
#define FRAC_BITS   23   /* fractional bits for sb_samples and dct */
#define WFRAC_BITS  16   /* fractional bits for window */
#else

#define SCALE_GEN(v) \
{ FIXR(1.0 * (v)), FIXR(0.7937005259 * (v)), FIXR(0.6299605249 * (v)) }

static INT32 scale_factor_mult2[3][3] = {
    SCALE_GEN(4.0 / 3.0), /* 3 steps */
    SCALE_GEN(4.0 / 5.0), /* 5 steps */
    SCALE_GEN(4.0 / 9.0), /* 9 steps */
};

/* 2^(n/4) */
static UINT32 scale_factor_mult3[4] = {
    FIXR(1.0),

    shift = scale_factor_modshift[scale_factor];
    mod = shift & 3;
    shift >>= 2;
    val = MUL64(mant + (-1 << n) + 1, scale_factor_mult[n-1][mod]);
    shift += n;
    /* NOTE: at this point, 1 <= shift >= 21 + 15 */
    return (int)((val + (1LL << (shift - 1))) >> shift);
}

static inline int l2_unscale_group(int steps, int mant, int scale_factor)
{
    int shift, mod, val;

    shift = scale_factor_modshift[scale_factor];
    mod = shift & 3;
    shift >>= 2;

    val = (mant - (steps >> 1)) * scale_factor_mult2[steps >> 2][mod];
    /* NOTE: at this point, 0 <= shift <= 21 */
    if (shift > 0)
        val = (val + (1 << (shift - 1))) >> shift;
    return val;
}

/* compute value^(4/3) * 2^(exponent/4). It normalized to FRAC_BITS */
static inline int l3_unscale(int value, int exponent)
{

    while (a < POW_FRAC_ONE) {
        a = a << 1;
        eq--;
    }
    /* now POW_FRAC_ONE <= a < 2 * POW_FRAC_ONE */
#if POW_FRAC_BITS > FRAC_BITS
    a = (a + (1 << (POW_FRAC_BITS - FRAC_BITS - 1))) >> (POW_FRAC_BITS - FRAC_BITS);
    /* correct overflow */
    if (a >= 2 * (1 << FRAC_BITS)) {
        a = a >> 1;
        eq++;

    if(!init) {
        /* scale factors table for layer 1/2 */
        for(i=0;i<64;i++) {
            int shift, mod;
            /* 1.0 (i = 3) is normalized to 2 ^ FRAC_BITS */
            shift = (i / 3);
            mod = i % 3;
            scale_factor_modshift[i] = mod | (shift << 2);
        }

        /* scale factor multiply for layer 1 */
        for(i=0;i<15;i++) {
            int n, norm;
            n = i + 2;
            norm = ((INT64_C(1) << n) * FRAC_ONE) / ((1 << n) - 1);
            scale_factor_mult[i][0] = MULL(FIXR(1.0 * 2.0), norm);
            scale_factor_mult[i][1] = MULL(FIXR(0.7937005259 * 2.0), norm);
            scale_factor_mult[i][2] = MULL(FIXR(0.6299605249 * 2.0), norm);
            dprintf("%d: norm=%x s=%x %x %x\n",
                    i, norm, 
                    scale_factor_mult[i][0],
                    scale_factor_mult[i][1],
                    scale_factor_mult[i][2]);

    synth_buf = s1->synth_buf[ch] + offset;

    for(j=0;j<32;j++) {
        v = tmp[j];
#if FRAC_BITS <= 15
        /* NOTE: can cause a loss in precision if very high amplitude
           sound */
        if (v > 32767)
            v = 32767;
        else if (v < -32768)
            v = -32768;
#endif

    
    s->layer = 4 - ((header >> 17) & 3);
    /* extract frequency */
    sample_rate_index = (header >> 10) & 3;
    sample_rate = mpa_freq_tab[sample_rate_index] >> (s->lsf + mpeg25);
    sample_rate_index += 3 * (s->lsf + mpeg25);
    s->sample_rate_index = sample_rate_index;
    s->error_protection = ((header >> 16) & 1) ^ 1;
    s->sample_rate = sample_rate;

    bitrate_index = (header >> 12) & 0xf;
    padding = (header >> 9) & 1;
    //extension = (header >> 8) & 1;
    s->mode = (header >> 6) & 3;

            s->frame_size += padding;
            s->bit_rate = (s->frame_size * (sample_rate << s->lsf)) / 144000;
            break;
        }
    }
    
#if defined(DEBUG)
    printf("layer%d, %d Hz, %d kbits/s, ",
           s->layer, s->sample_rate, s->bit_rate);
    if (s->nb_channels == 2) {

#if defined(DEBUG)
void sample_dump(int fnum, INT32 *tab, int n)
{
    static FILE *files[16], *f;
    char buf[512];
    int i;
    INT32 v;
    
    f = files[fnum];
    if (!f) {
        sprintf(buf, "/tmp/out%d.%s.pcm", 
                fnum, 
#ifdef USE_HIGHPRECISION
                "hp"
#else
                "lp"
#endif
                );
        f = fopen(buf, "w");
        if (!f)
            return;
        files[fnum] = f;
    }
    
    if (fnum == 0) {
        static int pos = 0;
        printf("pos=%d\n", pos);
        for(i=0;i<n;i++) {
            printf(" %0.4f", (double)tab[i] / FRAC_ONE);
            if ((i % 18) == 17)
                printf("\n");
        }
        pos += n;
    }
    for(i=0;i<n;i++) {
        /* normalize to 23 frac bits */
        v = tab[i] << (23 - FRAC_BITS);
        fwrite(&v, 1, sizeof(INT32), f);
    }
}
#endif


/* main layer3 decoding function */

            reorder_block(s, g);
#if defined(DEBUG)
            sample_dump(0, g->sb_hybrid, 576);
#endif
            compute_antialias(s, g);
#if defined(DEBUG)
            sample_dump(1, g->sb_hybrid, 576);
#endif
            compute_imdct(s, g, &s->sb_samples[ch][18 * gr][0], s->mdct_buf[ch]); 
#if defined(DEBUG)
            sample_dump(2, &s->sb_samples[ch][18 * gr][0], 576);
#endif
        }
    } /* gr */
    return nb_granules * 18;

                    /* reset free format frame size to give a chance
                       to get a new bitrate */
                    s->free_format_frame_size = 0;
		} else {
		    if (decode_header(s, header) == 1) {
                        /* free format: prepare to compute frame size */
			s->frame_size = -1;
                    }
                    /* update codec info */
                    avctx->sample_rate = s->sample_rate;
                    avctx->channels = s->nb_channels;
                    avctx->bit_rate = s->bit_rate;
                    avctx->frame_size = s->frame_size;
		}
	    }
        } else if (s->frame_size == -1) {
            /* free format : find next sync to compute frame size */
	    len = MPA_MAX_CODED_FRAME_SIZE - len;