view fft-test.c @ 2497:69adfbbdcdeb libavcodec

- samples from mplayer ftp in the "adv" profile seem to have profile=2, which isn't the advanced one; and indeed, using adv. profile parser fails. Using normal parser works, and that's what is done - attempt at taking care of stride for NORM2 bitplane decoding - duplication of much code from msmpeg4.c; this code isn't yet used, but goes down as far as the block layer (mainly Transform Type stuff, the remains are wild editing without checking). Unusable yet, and lacks the AC decoding (but a step further in bitstream parsing) patch by anonymous

/**
 * @file fft-test.c
 * FFT and MDCT tests.
 */

#include "dsputil.h"
#include <math.h>
#include <unistd.h>
#include <sys/time.h>

int mm_flags;

/* reference fft */

#define MUL16(a,b) ((a) * (b))

#define CMAC(pre, pim, are, aim, bre, bim) \
{\
   pre += (MUL16(are, bre) - MUL16(aim, bim));\
   pim += (MUL16(are, bim) + MUL16(bre, aim));\
}

FFTComplex *exptab;

void fft_ref_init(int nbits, int inverse)
{
    int n, i;
    float c1, s1, alpha;

    n = 1 << nbits;
    exptab = av_malloc((n / 2) * sizeof(FFTComplex));

    for(i=0;i<(n/2);i++) {
        alpha = 2 * M_PI * (float)i / (float)n;
        c1 = cos(alpha);
        s1 = sin(alpha);
        if (!inverse)
            s1 = -s1;
        exptab[i].re = c1;
        exptab[i].im = s1;
    }
}

void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
{
    int n, i, j, k, n2;
    float tmp_re, tmp_im, s, c;
    FFTComplex *q;

    n = 1 << nbits;
    n2 = n >> 1;
    for(i=0;i<n;i++) {
        tmp_re = 0;
        tmp_im = 0;
        q = tab;
        for(j=0;j<n;j++) {
            k = (i * j) & (n - 1);
            if (k >= n2) {
                c = -exptab[k - n2].re;
                s = -exptab[k - n2].im;
            } else {
                c = exptab[k].re;
                s = exptab[k].im;
            }
            CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
            q++;
        }
        tabr[i].re = tmp_re;
        tabr[i].im = tmp_im;
    }
}

void imdct_ref(float *out, float *in, int n)
{
    int k, i, a;
    float sum, f;

    for(i=0;i<n;i++) {
        sum = 0;
        for(k=0;k<n/2;k++) {
            a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
            f = cos(M_PI * a / (double)(2 * n));
            sum += f * in[k];
        }
        out[i] = -sum;
    }
}

/* NOTE: no normalisation by 1 / N is done */
void mdct_ref(float *output, float *input, int n)
{
    int k, i;
    float a, s;

    /* do it by hand */
    for(k=0;k<n/2;k++) {
        s = 0;
        for(i=0;i<n;i++) {
            a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
            s += input[i] * cos(a);
        }
        output[k] = s;
    }
}


float frandom(void)
{
    return (float)((random() & 0xffff) - 32768) / 32768.0;
}

int64_t gettime(void)
{
    struct timeval tv;
    gettimeofday(&tv,NULL);
    return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
}

void check_diff(float *tab1, float *tab2, int n)
{
    int i;

    for(i=0;i<n;i++) {
        if (fabsf(tab1[i] - tab2[i]) >= 1e-3) {
            printf("ERROR %d: %f %f\n", 
                   i, tab1[i], tab2[i]);
        }
    }
}


void help(void)
{
    printf("usage: fft-test [-h] [-s] [-i] [-n b]\n"
           "-h     print this help\n"
           "-s     speed test\n"
           "-m     (I)MDCT test\n"
           "-i     inverse transform test\n"
           "-n b   set the transform size to 2^b\n"
           );
    exit(1);
}



int main(int argc, char **argv)
{
    FFTComplex *tab, *tab1, *tab_ref;
    FFTSample *tabtmp, *tab2;
    int it, i, c;
    int do_speed = 0;
    int do_mdct = 0;
    int do_inverse = 0;
    FFTContext s1, *s = &s1;
    MDCTContext m1, *m = &m1;
    int fft_nbits, fft_size;

    mm_flags = 0;
    fft_nbits = 9;
    for(;;) {
        c = getopt(argc, argv, "hsimn:");
        if (c == -1)
            break;
        switch(c) {
        case 'h':
            help();
            break;
        case 's':
            do_speed = 1;
            break;
        case 'i':
            do_inverse = 1;
            break;
        case 'm':
            do_mdct = 1;
            break;
        case 'n':
            fft_nbits = atoi(optarg);
            break;
        }
    }

    fft_size = 1 << fft_nbits;
    tab = av_malloc(fft_size * sizeof(FFTComplex));
    tab1 = av_malloc(fft_size * sizeof(FFTComplex));
    tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
    tabtmp = av_malloc(fft_size / 2 * sizeof(FFTSample));
    tab2 = av_malloc(fft_size * sizeof(FFTSample));

    if (do_mdct) {
        if (do_inverse)
            printf("IMDCT");
        else
            printf("MDCT");
        ff_mdct_init(m, fft_nbits, do_inverse);
    } else {
        if (do_inverse)
            printf("IFFT");
        else
            printf("FFT");
        ff_fft_init(s, fft_nbits, do_inverse);
        fft_ref_init(fft_nbits, do_inverse);
    }
    printf(" %d test\n", fft_size);

    /* generate random data */

    for(i=0;i<fft_size;i++) {
        tab1[i].re = frandom();
        tab1[i].im = frandom();
    }

    /* checking result */
    printf("Checking...\n");

    if (do_mdct) {
        if (do_inverse) {
            imdct_ref((float *)tab_ref, (float *)tab1, fft_size);
            ff_imdct_calc(m, tab2, (float *)tab1, tabtmp);
            check_diff((float *)tab_ref, tab2, fft_size);
        } else {
            mdct_ref((float *)tab_ref, (float *)tab1, fft_size);
            
            ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);

            check_diff((float *)tab_ref, tab2, fft_size / 2);
        }
    } else {
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
        ff_fft_permute(s, tab);
        ff_fft_calc(s, tab);
        
        fft_ref(tab_ref, tab1, fft_nbits);
        check_diff((float *)tab_ref, (float *)tab, fft_size * 2);
    }

    /* do a speed test */

    if (do_speed) {
        int64_t time_start, duration;
        int nb_its;

        printf("Speed test...\n");
        /* we measure during about 1 seconds */
        nb_its = 1;
        for(;;) {
            time_start = gettime();
            for(it=0;it<nb_its;it++) {
                if (do_mdct) {
                    if (do_inverse) {
                        ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp);
                    } else {
                        ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp);
                    }
                } else {
                    memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
                    ff_fft_calc(s, tab);
                }
            }
            duration = gettime() - time_start;
            if (duration >= 1000000)
                break;
            nb_its *= 2;
        }
        printf("time: %0.1f us/transform [total time=%0.2f s its=%d]\n", 
               (double)duration / nb_its, 
               (double)duration / 1000000.0,
               nb_its);
    }
    
    if (do_mdct) {
        ff_mdct_end(m);
    } else {
        ff_fft_end(s);
    }
    return 0;
}