fft.c 16.1 KB
Newer Older
1 2 3
/*
 * (c) 2002 Fabrice Bellard
 *
4 5 6
 * This file is part of FFmpeg.
 *
 * FFmpeg 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
 * FFmpeg 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 FFmpeg; 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 27
#include "config.h"

28 29 30 31
#ifndef AVFFT
#define AVFFT 0
#endif

32
#include <math.h>
33
#if HAVE_UNISTD_H
Fabrice Bellard's avatar
Fabrice Bellard committed
34
#include <unistd.h>
35
#endif
36
#include <stdio.h>
37 38
#include <stdlib.h>
#include <string.h>
39

40 41 42 43 44 45
#include "libavutil/cpu.h"
#include "libavutil/lfg.h"
#include "libavutil/log.h"
#include "libavutil/mathematics.h"
#include "libavutil/time.h"

46 47 48
#if AVFFT
#include "libavcodec/avfft.h"
#else
49
#include "libavcodec/fft.h"
50 51
#endif

52
#if FFT_FLOAT
53 54
#include "libavcodec/dct.h"
#include "libavcodec/rdft.h"
55
#endif
56

57 58
/* reference fft */

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

61 62 63 64 65
#define CMAC(pre, pim, are, aim, bre, bim)          \
    {                                               \
        pre += (MUL16(are, bre) - MUL16(aim, bim)); \
        pim += (MUL16(are, bim) + MUL16(bre, aim)); \
    }
66

67
#if FFT_FLOAT || AVFFT
68 69 70
#define RANGE 1.0
#define REF_SCALE(x, bits)  (x)
#define FMT "%10.6f"
71
#elif FFT_FIXED_32
72 73 74
#define RANGE 8388608
#define REF_SCALE(x, bits) (x)
#define FMT "%6d"
75
#else
76 77 78
#define RANGE 16384
#define REF_SCALE(x, bits) ((x) / (1 << (bits)))
#define FMT "%6d"
79 80
#endif

81
static struct {
82 83
    float re, im;
} *exptab;
84

85
static int fft_ref_init(int nbits, int inverse)
86
{
87
    int i, n = 1 << nbits;
88

89
    exptab = av_malloc_array((n / 2), sizeof(*exptab));
90 91
    if (!exptab)
        return AVERROR(ENOMEM);
92

93 94
    for (i = 0; i < (n / 2); i++) {
        double alpha = 2 * M_PI * (float) i / (float) n;
95
        double c1 = cos(alpha), s1 = sin(alpha);
96 97 98 99 100
        if (!inverse)
            s1 = -s1;
        exptab[i].re = c1;
        exptab[i].im = s1;
    }
101
    return 0;
102 103
}

104
static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
105
{
106 107 108
    int i, j;
    int n  = 1 << nbits;
    int n2 = n >> 1;
109

110
    for (i = 0; i < n; i++) {
111 112
        double tmp_re = 0, tmp_im = 0;
        FFTComplex *q = tab;
113
        for (j = 0; j < n; j++) {
114 115
            double s, c;
            int k = (i * j) & (n - 1);
116 117 118 119 120 121 122 123 124 125
            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++;
        }
126 127
        tabr[i].re = REF_SCALE(tmp_re, nbits);
        tabr[i].im = REF_SCALE(tmp_im, nbits);
128 129 130
    }
}

131
#if CONFIG_MDCT
132
static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
133
{
134
    int i, k, n = 1 << nbits;
135

136
    for (i = 0; i < n; i++) {
137
        double sum = 0;
138
        for (k = 0; k < n / 2; k++) {
139
            int a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
140
            double f = cos(M_PI * a / (double) (2 * n));
141 142
            sum += f * in[k];
        }
143
        out[i] = REF_SCALE(-sum, nbits - 2);
144 145 146 147
    }
}

/* NOTE: no normalisation by 1 / N is done */
148
static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
149
{
150
    int i, k, n = 1 << nbits;
151 152

    /* do it by hand */
153
    for (k = 0; k < n / 2; k++) {
154
        double s = 0;
155
        for (i = 0; i < n; i++) {
156
            double a = (2 * M_PI * (2 * i + 1 + n / 2) * (2 * k + 1) / (4 * n));
157 158
            s += input[i] * cos(a);
        }
159
        output[k] = REF_SCALE(s, nbits - 1);
160 161
    }
}
162
#endif /* CONFIG_MDCT */
163

164
#if FFT_FLOAT
165
#if CONFIG_DCT
166
static void idct_ref(FFTSample *output, FFTSample *input, int nbits)
167
{
168
    int i, k, n = 1 << nbits;
169 170 171

    /* do it by hand */
    for (i = 0; i < n; i++) {
172
        double s = 0.5 * input[0];
173
        for (k = 1; k < n; k++) {
174
            double a = M_PI * k * (i + 0.5) / n;
175 176 177 178 179
            s += input[k] * cos(a);
        }
        output[i] = 2 * s / n;
    }
}
180

181
static void dct_ref(FFTSample *output, FFTSample *input, int nbits)
182
{
183
    int i, k, n = 1 << nbits;
184 185 186

    /* do it by hand */
    for (k = 0; k < n; k++) {
187
        double s = 0;
188
        for (i = 0; i < n; i++) {
189
            double a = M_PI * k * (i + 0.5) / n;
190 191 192 193 194
            s += input[i] * cos(a);
        }
        output[k] = s;
    }
}
195
#endif /* CONFIG_DCT */
196
#endif /* FFT_FLOAT */
197

198
static FFTSample frandom(AVLFG *prng)
199
{
200
    return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE;
201 202
}

203
static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
204
{
205 206
    int i, err = 0;
    double error = 0, max = 0;
207

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

223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
static inline void fft_init(FFTContext **s, int nbits, int inverse)
{
#if AVFFT
    *s = av_fft_init(nbits, inverse);
#else
    ff_fft_init(*s, nbits, inverse);
#endif
}

static inline void mdct_init(FFTContext **s, int nbits, int inverse, double scale)
{
#if AVFFT
    *s = av_mdct_init(nbits, inverse, scale);
#else
    ff_mdct_init(*s, nbits, inverse, scale);
#endif
}

static inline void mdct_calc(FFTContext *s, FFTSample *output, const FFTSample *input)
{
#if AVFFT
    av_mdct_calc(s, output, input);
#else
    s->mdct_calc(s, output, input);
#endif
}

static inline void imdct_calc(struct FFTContext *s, FFTSample *output, const FFTSample *input)
{
#if AVFFT
    av_imdct_calc(s, output, input);
#else
    s->imdct_calc(s, output, input);
#endif
}

static inline void fft_permute(FFTContext *s, FFTComplex *z)
{
#if AVFFT
    av_fft_permute(s, z);
#else
    s->fft_permute(s, z);
#endif
}

static inline void fft_calc(FFTContext *s, FFTComplex *z)
{
#if AVFFT
    av_fft_calc(s, z);
#else
    s->fft_calc(s, z);
#endif
}

static inline void mdct_end(FFTContext *s)
{
#if AVFFT
    av_mdct_end(s);
#else
    ff_mdct_end(s);
#endif
}

static inline void fft_end(FFTContext *s)
{
#if AVFFT
    av_fft_end(s);
#else
    ff_fft_end(s);
#endif
}

#if FFT_FLOAT
static inline void rdft_init(RDFTContext **r, int nbits, enum RDFTransformType trans)
{
#if AVFFT
    *r = av_rdft_init(nbits, trans);
#else
    ff_rdft_init(*r, nbits, trans);
#endif
}

static inline void dct_init(DCTContext **d, int nbits, enum DCTTransformType trans)
{
#if AVFFT
    *d = av_dct_init(nbits, trans);
#else
    ff_dct_init(*d, nbits, trans);
#endif
}

static inline void rdft_calc(RDFTContext *r, FFTSample *tab)
{
#if AVFFT
    av_rdft_calc(r, tab);
#else
    r->rdft_calc(r, tab);
#endif
}

static inline void dct_calc(DCTContext *d, FFTSample *data)
{
#if AVFFT
    av_dct_calc(d, data);
#else
    d->dct_calc(d, data);
#endif
}

static inline void rdft_end(RDFTContext *r)
{
#if AVFFT
    av_rdft_end(r);
#else
    ff_rdft_end(r);
#endif
}

static inline void dct_end(DCTContext *d)
{
#if AVFFT
    av_dct_end(d);
#else
    ff_dct_end(d);
#endif
}
#endif /* FFT_FLOAT */

351
static void help(void)
352
{
353 354
    av_log(NULL, AV_LOG_INFO,
           "usage: fft-test [-h] [-s] [-i] [-n b]\n"
355 356 357
           "-h     print this help\n"
           "-s     speed test\n"
           "-m     (I)MDCT test\n"
358
           "-d     (I)DCT test\n"
359
           "-r     (I)RDFT test\n"
360 361
           "-i     inverse transform test\n"
           "-n b   set the transform size to 2^b\n"
362
           "-f x   set scale factor for output data of (I)MDCT to x\n");
363 364
}

365 366 367
enum tf_transform {
    TRANSFORM_FFT,
    TRANSFORM_MDCT,
368
    TRANSFORM_RDFT,
369
    TRANSFORM_DCT,
370
};
371

372 373 374 375
#if !HAVE_GETOPT
#include "compat/getopt.c"
#endif

376 377 378
int main(int argc, char **argv)
{
    FFTComplex *tab, *tab1, *tab_ref;
Loren Merritt's avatar
Loren Merritt committed
379
    FFTSample *tab2;
380
    enum tf_transform transform = TRANSFORM_FFT;
381
    FFTContext *m, *s;
382
#if FFT_FLOAT
383 384
    RDFTContext *r;
    DCTContext *d;
385
#endif /* FFT_FLOAT */
386 387 388
    int it, i, err = 1;
    int do_speed = 0, do_inverse = 0;
    int fft_nbits = 9, fft_size;
389
    double scale = 1.0;
390
    AVLFG prng;
391

392 393 394 395 396 397 398 399 400 401
#if !AVFFT
    s = av_mallocz(sizeof(*s));
    m = av_mallocz(sizeof(*m));
#endif

#if !AVFFT && FFT_FLOAT
    r = av_mallocz(sizeof(*r));
    d = av_mallocz(sizeof(*d));
#endif

402
    av_lfg_init(&prng, 1);
403

404
    for (;;) {
405
        int c = getopt(argc, argv, "hsimrdn:f:c:");
406 407
        if (c == -1)
            break;
408
        switch (c) {
409 410
        case 'h':
            help();
411
            return 1;
412 413 414 415 416 417 418
        case 's':
            do_speed = 1;
            break;
        case 'i':
            do_inverse = 1;
            break;
        case 'm':
419
            transform = TRANSFORM_MDCT;
420
            break;
421 422 423
        case 'r':
            transform = TRANSFORM_RDFT;
            break;
424 425 426
        case 'd':
            transform = TRANSFORM_DCT;
            break;
427 428 429
        case 'n':
            fft_nbits = atoi(optarg);
            break;
430 431 432
        case 'f':
            scale = atof(optarg);
            break;
433
        case 'c':
434
        {
435
            unsigned cpuflags = av_get_cpu_flags();
436 437

            if (av_parse_cpu_caps(&cpuflags, optarg) < 0)
438
                return 1;
439 440

            av_force_cpu_flags(cpuflags);
441
            break;
442
        }
443
        }
444 445 446
    }

    fft_size = 1 << fft_nbits;
447 448 449 450
    tab      = av_malloc_array(fft_size, sizeof(FFTComplex));
    tab1     = av_malloc_array(fft_size, sizeof(FFTComplex));
    tab_ref  = av_malloc_array(fft_size, sizeof(FFTComplex));
    tab2     = av_malloc_array(fft_size, sizeof(FFTSample));
451

452 453 454
    if (!(tab && tab1 && tab_ref && tab2))
        goto cleanup;

455
    switch (transform) {
456
#if CONFIG_MDCT
457
    case TRANSFORM_MDCT:
458
        av_log(NULL, AV_LOG_INFO, "Scale factor is set to %f\n", scale);
459
        if (do_inverse)
460
            av_log(NULL, AV_LOG_INFO, "IMDCT");
461
        else
462
            av_log(NULL, AV_LOG_INFO, "MDCT");
463
        mdct_init(&m, fft_nbits, do_inverse, scale);
464
        break;
465
#endif /* CONFIG_MDCT */
466
    case TRANSFORM_FFT:
467
        if (do_inverse)
468
            av_log(NULL, AV_LOG_INFO, "IFFT");
469
        else
470
            av_log(NULL, AV_LOG_INFO, "FFT");
471
        fft_init(&s, fft_nbits, do_inverse);
472
        if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
473
            goto cleanup;
474
        break;
475
#if FFT_FLOAT
476
#    if CONFIG_RDFT
477 478
    case TRANSFORM_RDFT:
        if (do_inverse)
479
            av_log(NULL, AV_LOG_INFO, "IDFT_C2R");
480
        else
481
            av_log(NULL, AV_LOG_INFO, "DFT_R2C");
482
        rdft_init(&r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
483
        if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
484
            goto cleanup;
485
        break;
486
#    endif /* CONFIG_RDFT */
487
#    if CONFIG_DCT
488 489
    case TRANSFORM_DCT:
        if (do_inverse)
490
            av_log(NULL, AV_LOG_INFO, "DCT_III");
491
        else
492
            av_log(NULL, AV_LOG_INFO, "DCT_II");
493
        dct_init(&d, fft_nbits, do_inverse ? DCT_III : DCT_II);
494
        break;
495
#    endif /* CONFIG_DCT */
496
#endif /* FFT_FLOAT */
497 498
    default:
        av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
499
        goto cleanup;
500
    }
501
    av_log(NULL, AV_LOG_INFO, " %d test\n", fft_size);
502 503 504

    /* generate random data */

505
    for (i = 0; i < fft_size; i++) {
506 507
        tab1[i].re = frandom(&prng);
        tab1[i].im = frandom(&prng);
508 509 510
    }

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

513
    switch (transform) {
514
#if CONFIG_MDCT
515
    case TRANSFORM_MDCT:
516
        if (do_inverse) {
517
            imdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
518
            imdct_calc(m, tab2, &tab1->re);
519
            err = check_diff(&tab_ref->re, tab2, fft_size, scale);
520
        } else {
521
            mdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
522
            mdct_calc(m, tab2, &tab1->re);
523
            err = check_diff(&tab_ref->re, tab2, fft_size / 2, scale);
524
        }
525
        break;
526
#endif /* CONFIG_MDCT */
527
    case TRANSFORM_FFT:
528
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
529 530
        fft_permute(s, tab);
        fft_calc(s, tab);
531

532
        fft_ref(tab_ref, tab1, fft_nbits);
533
        err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 1.0);
534
        break;
535
#if FFT_FLOAT
536
#if CONFIG_RDFT
537
    case TRANSFORM_RDFT:
538 539
    {
        int fft_size_2 = fft_size >> 1;
540
        if (do_inverse) {
541
            tab1[0].im          = 0;
542 543
            tab1[fft_size_2].im = 0;
            for (i = 1; i < fft_size_2; i++) {
544 545
                tab1[fft_size_2 + i].re =  tab1[fft_size_2 - i].re;
                tab1[fft_size_2 + i].im = -tab1[fft_size_2 - i].im;
546 547 548 549 550
            }

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

551
            rdft_calc(r, tab2);
552 553 554 555 556
            fft_ref(tab_ref, tab1, fft_nbits);
            for (i = 0; i < fft_size; i++) {
                tab[i].re = tab2[i];
                tab[i].im = 0;
            }
557
            err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 0.5);
558 559 560 561 562
        } else {
            for (i = 0; i < fft_size; i++) {
                tab2[i]    = tab1[i].re;
                tab1[i].im = 0;
            }
563
            rdft_calc(r, tab2);
564 565
            fft_ref(tab_ref, tab1, fft_nbits);
            tab_ref[0].im = tab_ref[fft_size_2].re;
566
            err = check_diff(&tab_ref->re, tab2, fft_size, 1.0);
567
        }
568
        break;
569
    }
570 571
#endif /* CONFIG_RDFT */
#if CONFIG_DCT
572 573
    case TRANSFORM_DCT:
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
574
        dct_calc(d, &tab->re);
575
        if (do_inverse)
576
            idct_ref(&tab_ref->re, &tab1->re, fft_nbits);
577
        else
578
            dct_ref(&tab_ref->re, &tab1->re, fft_nbits);
579
        err = check_diff(&tab_ref->re, &tab->re, fft_size, 1.0);
580
        break;
581
#endif /* CONFIG_DCT */
582
#endif /* FFT_FLOAT */
583 584 585 586 587
    }

    /* do a speed test */

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

591
        av_log(NULL, AV_LOG_INFO, "Speed test...\n");
592 593
        /* we measure during about 1 seconds */
        nb_its = 1;
594
        for (;;) {
595
            time_start = av_gettime_relative();
596
            for (it = 0; it < nb_its; it++) {
597 598
                switch (transform) {
                case TRANSFORM_MDCT:
599
                    if (do_inverse)
600
                        imdct_calc(m, &tab->re, &tab1->re);
601
                    else
602
                        mdct_calc(m, &tab->re, &tab1->re);
603 604
                    break;
                case TRANSFORM_FFT:
605
                    memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
606
                    fft_calc(s, tab);
607
                    break;
608
#if FFT_FLOAT
609 610
                case TRANSFORM_RDFT:
                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
611
                    rdft_calc(r, tab2);
612
                    break;
613 614
                case TRANSFORM_DCT:
                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
615
                    dct_calc(d, tab2);
616
                    break;
617
#endif /* FFT_FLOAT */
618 619
                }
            }
620
            duration = av_gettime_relative() - time_start;
621 622 623 624
            if (duration >= 1000000)
                break;
            nb_its *= 2;
        }
625 626 627 628
        av_log(NULL, AV_LOG_INFO,
               "time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
               (double) duration / nb_its,
               (double) duration / 1000000.0,
629 630
               nb_its);
    }
631

632
    switch (transform) {
633
#if CONFIG_MDCT
634
    case TRANSFORM_MDCT:
635
        mdct_end(m);
636
        break;
637
#endif /* CONFIG_MDCT */
638
    case TRANSFORM_FFT:
639
        fft_end(s);
640
        break;
641
#if FFT_FLOAT
642
#    if CONFIG_RDFT
643
    case TRANSFORM_RDFT:
644
        rdft_end(r);
645
        break;
646
#    endif /* CONFIG_RDFT */
647
#    if CONFIG_DCT
648
    case TRANSFORM_DCT:
649
        dct_end(d);
650
        break;
651
#    endif /* CONFIG_DCT */
652
#endif /* FFT_FLOAT */
653
    }
654

655
cleanup:
656 657 658 659 660 661
    av_free(tab);
    av_free(tab1);
    av_free(tab2);
    av_free(tab_ref);
    av_free(exptab);

662 663 664 665 666 667 668 669 670 671
#if !AVFFT
    av_free(s);
    av_free(m);
#endif

#if !AVFFT && FFT_FLOAT
    av_free(r);
    av_free(d);
#endif

672 673 674 675
    if (err)
        printf("Error: %d.\n", err);

    return !!err;
676
}