fft-test.c 12.7 KB
Newer Older
1 2 3
/*
 * (c) 2002 Fabrice Bellard
 *
4
 * This file is part of Libav.
5
 *
6
 * Libav is free software; you can redistribute it and/or
7 8
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
9
 * version 2.1 of the License, or (at your option) any later version.
10
 *
11
 * Libav is distributed in the hope that it will be useful,
12 13 14 15 16
 * 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
17
 * License along with Libav; if not, write to the Free Software
18 19 20
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

Michael Niedermayer's avatar
Michael Niedermayer committed
21
/**
22
 * @file
Michael Niedermayer's avatar
Michael Niedermayer committed
23 24 25
 * FFT and MDCT tests.
 */

26
#include "libavutil/mathematics.h"
27
#include "libavutil/lfg.h"
28 29
#include "libavutil/log.h"
#include "fft.h"
30
#if CONFIG_FFT_FLOAT
31 32
#include "dct.h"
#include "rdft.h"
33
#endif
34
#include <math.h>
Fabrice Bellard's avatar
Fabrice Bellard committed
35
#include <unistd.h>
36
#include <sys/time.h>
37 38
#include <stdlib.h>
#include <string.h>
39

40 41
#undef exit

42 43 44 45 46 47 48 49 50 51
/* 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));\
}

52 53 54 55 56 57 58 59 60 61 62 63 64
#if CONFIG_FFT_FLOAT
#   define RANGE 1.0
#   define REF_SCALE(x, bits)  (x)
#   define FMT "%10.6f"
#else
#   define RANGE 16384
#   define REF_SCALE(x, bits) ((x) / (1<<(bits)))
#   define FMT "%6d"
#endif

struct {
    float re, im;
} *exptab;
65

66
static void fft_ref_init(int nbits, int inverse)
67 68
{
    int n, i;
69
    double c1, s1, alpha;
70 71

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

74
    for (i = 0; i < (n/2); i++) {
75 76 77 78 79 80 81 82 83 84
        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;
    }
}

85
static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
86 87
{
    int n, i, j, k, n2;
88
    double tmp_re, tmp_im, s, c;
89 90 91 92
    FFTComplex *q;

    n = 1 << nbits;
    n2 = n >> 1;
93
    for (i = 0; i < n; i++) {
94 95 96
        tmp_re = 0;
        tmp_im = 0;
        q = tab;
97
        for (j = 0; j < n; j++) {
98 99 100 101 102 103 104 105 106 107 108
            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++;
        }
109 110
        tabr[i].re = REF_SCALE(tmp_re, nbits);
        tabr[i].im = REF_SCALE(tmp_im, nbits);
111 112 113
    }
}

114
static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
115
{
116
    int n = 1<<nbits;
117
    int k, i, a;
118
    double sum, f;
119

120
    for (i = 0; i < n; i++) {
121
        sum = 0;
122
        for (k = 0; k < n/2; k++) {
123 124 125 126
            a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
            f = cos(M_PI * a / (double)(2 * n));
            sum += f * in[k];
        }
127
        out[i] = REF_SCALE(-sum, nbits - 2);
128 129 130 131
    }
}

/* NOTE: no normalisation by 1 / N is done */
132
static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
133
{
134
    int n = 1<<nbits;
135
    int k, i;
136
    double a, s;
137 138

    /* do it by hand */
139
    for (k = 0; k < n/2; k++) {
140
        s = 0;
141
        for (i = 0; i < n; i++) {
142 143 144
            a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
            s += input[i] * cos(a);
        }
145
        output[k] = REF_SCALE(s, nbits - 1);
146 147 148
    }
}

149
#if CONFIG_FFT_FLOAT
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181
static void idct_ref(float *output, float *input, int nbits)
{
    int n = 1<<nbits;
    int k, i;
    double a, s;

    /* do it by hand */
    for (i = 0; i < n; i++) {
        s = 0.5 * input[0];
        for (k = 1; k < n; k++) {
            a = M_PI*k*(i+0.5) / n;
            s += input[k] * cos(a);
        }
        output[i] = 2 * s / n;
    }
}
static void dct_ref(float *output, float *input, int nbits)
{
    int n = 1<<nbits;
    int k, i;
    double a, s;

    /* do it by hand */
    for (k = 0; k < n; k++) {
        s = 0;
        for (i = 0; i < n; i++) {
            a = M_PI*k*(i+0.5) / n;
            s += input[i] * cos(a);
        }
        output[k] = s;
    }
}
182
#endif
183

184

185
static FFTSample frandom(AVLFG *prng)
186
{
187
    return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;
188 189
}

190
static int64_t gettime(void)
191 192 193
{
    struct timeval tv;
    gettimeofday(&tv,NULL);
194
    return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
195 196
}

197
static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
198 199
{
    int i;
Michael Niedermayer's avatar
Michael Niedermayer committed
200 201
    double max= 0;
    double error= 0;
202
    int err = 0;
203

204
    for (i = 0; i < n; i++) {
205
        double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;
Michael Niedermayer's avatar
Michael Niedermayer committed
206
        if (e >= 1e-3) {
207
            av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
208
                   i, tab1[i], tab2[i]);
209
            err = 1;
210
        }
Michael Niedermayer's avatar
Michael Niedermayer committed
211 212
        error+= e*e;
        if(e>max) max= e;
213
    }
Michael Niedermayer's avatar
Michael Niedermayer committed
214
    av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n);
215
    return err;
216 217 218
}


219
static void help(void)
220
{
221
    av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
222 223 224
           "-h     print this help\n"
           "-s     speed test\n"
           "-m     (I)MDCT test\n"
225
           "-d     (I)DCT test\n"
226
           "-r     (I)RDFT test\n"
227 228
           "-i     inverse transform test\n"
           "-n b   set the transform size to 2^b\n"
229
           "-f x   set scale factor for output data of (I)MDCT to x\n"
230 231 232 233
           );
    exit(1);
}

234 235 236
enum tf_transform {
    TRANSFORM_FFT,
    TRANSFORM_MDCT,
237
    TRANSFORM_RDFT,
238
    TRANSFORM_DCT,
239
};
240 241 242 243

int main(int argc, char **argv)
{
    FFTComplex *tab, *tab1, *tab_ref;
Loren Merritt's avatar
Loren Merritt committed
244
    FFTSample *tab2;
245 246
    int it, i, c;
    int do_speed = 0;
247
    int err = 1;
248
    enum tf_transform transform = TRANSFORM_FFT;
249 250
    int do_inverse = 0;
    FFTContext s1, *s = &s1;
251
    FFTContext m1, *m = &m1;
252
#if CONFIG_FFT_FLOAT
253
    RDFTContext r1, *r = &r1;
254
    DCTContext d1, *d = &d1;
255
    int fft_size_2;
256
#endif
257
    int fft_nbits, fft_size;
258
    double scale = 1.0;
259 260
    AVLFG prng;
    av_lfg_init(&prng, 1);
261 262 263

    fft_nbits = 9;
    for(;;) {
264
        c = getopt(argc, argv, "hsimrdn:f:");
265 266 267 268 269 270 271 272 273 274 275 276 277
        if (c == -1)
            break;
        switch(c) {
        case 'h':
            help();
            break;
        case 's':
            do_speed = 1;
            break;
        case 'i':
            do_inverse = 1;
            break;
        case 'm':
278
            transform = TRANSFORM_MDCT;
279
            break;
280 281 282
        case 'r':
            transform = TRANSFORM_RDFT;
            break;
283 284 285
        case 'd':
            transform = TRANSFORM_DCT;
            break;
286 287 288
        case 'n':
            fft_nbits = atoi(optarg);
            break;
289 290 291
        case 'f':
            scale = atof(optarg);
            break;
292 293 294 295 296 297 298 299 300
        }
    }

    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));
    tab2 = av_malloc(fft_size * sizeof(FFTSample));

301 302
    switch (transform) {
    case TRANSFORM_MDCT:
303
        av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale);
304
        if (do_inverse)
305
            av_log(NULL, AV_LOG_INFO,"IMDCT");
306
        else
307
            av_log(NULL, AV_LOG_INFO,"MDCT");
308
        ff_mdct_init(m, fft_nbits, do_inverse, scale);
309 310
        break;
    case TRANSFORM_FFT:
311
        if (do_inverse)
312
            av_log(NULL, AV_LOG_INFO,"IFFT");
313
        else
314
            av_log(NULL, AV_LOG_INFO,"FFT");
315
        ff_fft_init(s, fft_nbits, do_inverse);
316
        fft_ref_init(fft_nbits, do_inverse);
317
        break;
318
#if CONFIG_FFT_FLOAT
319 320
    case TRANSFORM_RDFT:
        if (do_inverse)
321
            av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
322
        else
323 324
            av_log(NULL, AV_LOG_INFO,"DFT_R2C");
        ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
325 326
        fft_ref_init(fft_nbits, do_inverse);
        break;
327 328
    case TRANSFORM_DCT:
        if (do_inverse)
329
            av_log(NULL, AV_LOG_INFO,"DCT_III");
330
        else
331 332
            av_log(NULL, AV_LOG_INFO,"DCT_II");
        ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);
333
        break;
334 335 336 337
#endif
    default:
        av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
        return 1;
338
    }
339
    av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
340 341 342

    /* generate random data */

343
    for (i = 0; i < fft_size; i++) {
344 345
        tab1[i].re = frandom(&prng);
        tab1[i].im = frandom(&prng);
346 347 348
    }

    /* checking result */
349
    av_log(NULL, AV_LOG_INFO,"Checking...\n");
350

351 352
    switch (transform) {
    case TRANSFORM_MDCT:
353
        if (do_inverse) {
354 355 356
            imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
            m->imdct_calc(m, tab2, (FFTSample *)tab1);
            err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);
357
        } else {
358
            mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
359

360
            m->mdct_calc(m, tab2, (FFTSample *)tab1);
361

362
            err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);
363
        }
364 365
        break;
    case TRANSFORM_FFT:
366
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
367 368
        s->fft_permute(s, tab);
        s->fft_calc(s, tab);
369

370
        fft_ref(tab_ref, tab1, fft_nbits);
371
        err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);
372
        break;
373
#if CONFIG_FFT_FLOAT
374
    case TRANSFORM_RDFT:
375
        fft_size_2 = fft_size >> 1;
376 377 378 379 380 381 382 383 384 385 386
        if (do_inverse) {
            tab1[         0].im = 0;
            tab1[fft_size_2].im = 0;
            for (i = 1; i < fft_size_2; i++) {
                tab1[fft_size_2+i].re =  tab1[fft_size_2-i].re;
                tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;
            }

            memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
            tab2[1] = tab1[fft_size_2].re;

387
            r->rdft_calc(r, tab2);
388 389 390 391 392
            fft_ref(tab_ref, tab1, fft_nbits);
            for (i = 0; i < fft_size; i++) {
                tab[i].re = tab2[i];
                tab[i].im = 0;
            }
393
            err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5);
394 395 396 397 398
        } else {
            for (i = 0; i < fft_size; i++) {
                tab2[i]    = tab1[i].re;
                tab1[i].im = 0;
            }
399
            r->rdft_calc(r, tab2);
400 401
            fft_ref(tab_ref, tab1, fft_nbits);
            tab_ref[0].im = tab_ref[fft_size_2].re;
402
            err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);
403
        }
404 405 406
        break;
    case TRANSFORM_DCT:
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
407
        d->dct_calc(d, tab);
408 409 410 411 412
        if (do_inverse) {
            idct_ref(tab_ref, tab1, fft_nbits);
        } else {
            dct_ref(tab_ref, tab1, fft_nbits);
        }
413
        err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);
414
        break;
415
#endif
416 417 418 419 420
    }

    /* do a speed test */

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

424
        av_log(NULL, AV_LOG_INFO,"Speed test...\n");
425 426 427 428
        /* we measure during about 1 seconds */
        nb_its = 1;
        for(;;) {
            time_start = gettime();
429
            for (it = 0; it < nb_its; it++) {
430 431
                switch (transform) {
                case TRANSFORM_MDCT:
432
                    if (do_inverse) {
433
                        m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
434
                    } else {
435
                        m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
436
                    }
437 438
                    break;
                case TRANSFORM_FFT:
439
                    memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
440
                    s->fft_calc(s, tab);
441
                    break;
442
#if CONFIG_FFT_FLOAT
443 444
                case TRANSFORM_RDFT:
                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
445
                    r->rdft_calc(r, tab2);
446
                    break;
447 448
                case TRANSFORM_DCT:
                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
449
                    d->dct_calc(d, tab2);
450
                    break;
451
#endif
452 453 454 455 456 457 458
                }
            }
            duration = gettime() - time_start;
            if (duration >= 1000000)
                break;
            nb_its *= 2;
        }
459 460
        av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
               (double)duration / nb_its,
461 462 463
               (double)duration / 1000000.0,
               nb_its);
    }
464

465 466
    switch (transform) {
    case TRANSFORM_MDCT:
Fabrice Bellard's avatar
Fabrice Bellard committed
467
        ff_mdct_end(m);
468 469
        break;
    case TRANSFORM_FFT:
470
        ff_fft_end(s);
471
        break;
472
#if CONFIG_FFT_FLOAT
473 474 475
    case TRANSFORM_RDFT:
        ff_rdft_end(r);
        break;
476 477 478
    case TRANSFORM_DCT:
        ff_dct_end(d);
        break;
479
#endif
480
    }
481 482 483 484 485 486 487

    av_free(tab);
    av_free(tab1);
    av_free(tab2);
    av_free(tab_ref);
    av_free(exptab);

488
    return err;
489
}