Statistics
| Branch: | Revision:

ffmpeg / libavcodec / fft-test.c @ 6efe6028

History | View | Annotate | Download (12.6 KB)

1
/*
2
 * (c) 2002 Fabrice Bellard
3
 *
4
 * This file is part of Libav.
5
 *
6
 * Libav is free software; you can redistribute it and/or
7
 * modify it under the terms of the GNU Lesser General Public
8
 * 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
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14
 * Lesser General Public License for more details.
15
 *
16
 * 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
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19
 */
20

    
21
/**
22
 * @file
23
 * FFT and MDCT tests.
24
 */
25

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

    
40
#undef exit
41

    
42
/* reference fft */
43

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

    
46
#define CMAC(pre, pim, are, aim, bre, bim) \
47
{\
48
   pre += (MUL16(are, bre) - MUL16(aim, bim));\
49
   pim += (MUL16(are, bim) + MUL16(bre, aim));\
50
}
51

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

    
62
struct {
63
    float re, im;
64
} *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
        alpha = 2 * M_PI * (float)i / (float)n;
76
        c1 = cos(alpha);
77
        s1 = sin(alpha);
78
        if (!inverse)
79
            s1 = -s1;
80
        exptab[i].re = c1;
81
        exptab[i].im = s1;
82
    }
83
}
84

    
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
    FFTComplex *q;
90

    
91
    n = 1 << nbits;
92
    n2 = n >> 1;
93
    for (i = 0; i < n; i++) {
94
        tmp_re = 0;
95
        tmp_im = 0;
96
        q = tab;
97
        for (j = 0; j < n; j++) {
98
            k = (i * j) & (n - 1);
99
            if (k >= n2) {
100
                c = -exptab[k - n2].re;
101
                s = -exptab[k - n2].im;
102
            } else {
103
                c = exptab[k].re;
104
                s = exptab[k].im;
105
            }
106
            CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
107
            q++;
108
        }
109
        tabr[i].re = REF_SCALE(tmp_re, nbits);
110
        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
            a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
124
            f = cos(M_PI * a / (double)(2 * n));
125
            sum += f * in[k];
126
        }
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
            a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
143
            s += input[i] * cos(a);
144
        }
145
        output[k] = REF_SCALE(s, nbits - 1);
146
    }
147
}
148

    
149
#if CONFIG_FFT_FLOAT
150
static void idct_ref(float *output, float *input, int nbits)
151
{
152
    int n = 1<<nbits;
153
    int k, i;
154
    double a, s;
155

    
156
    /* do it by hand */
157
    for (i = 0; i < n; i++) {
158
        s = 0.5 * input[0];
159
        for (k = 1; k < n; k++) {
160
            a = M_PI*k*(i+0.5) / n;
161
            s += input[k] * cos(a);
162
        }
163
        output[i] = 2 * s / n;
164
    }
165
}
166
static void dct_ref(float *output, float *input, int nbits)
167
{
168
    int n = 1<<nbits;
169
    int k, i;
170
    double a, s;
171

    
172
    /* do it by hand */
173
    for (k = 0; k < n; k++) {
174
        s = 0;
175
        for (i = 0; i < n; i++) {
176
            a = M_PI*k*(i+0.5) / n;
177
            s += input[i] * cos(a);
178
        }
179
        output[k] = s;
180
    }
181
}
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
    struct timeval tv;
193
    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;
200
    double max= 0;
201
    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;
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
        }
211
        error+= e*e;
212
        if(e>max) max= e;
213
    }
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
           "-h     print this help\n"
223
           "-s     speed test\n"
224
           "-m     (I)MDCT test\n"
225
           "-d     (I)DCT test\n"
226
           "-r     (I)RDFT test\n"
227
           "-i     inverse transform test\n"
228
           "-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
    exit(1);
232
}
233

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

    
241
int main(int argc, char **argv)
242
{
243
    FFTComplex *tab, *tab1, *tab_ref;
244
    FFTSample *tab2;
245
    int it, i, c;
246
    int do_speed = 0;
247
    int err = 1;
248
    enum tf_transform transform = TRANSFORM_FFT;
249
    int do_inverse = 0;
250
    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
#endif
256
    int fft_nbits, fft_size, fft_size_2;
257
    double scale = 1.0;
258
    AVLFG prng;
259
    av_lfg_init(&prng, 1);
260

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

    
294
    fft_size = 1 << fft_nbits;
295
    fft_size_2 = fft_size >> 1;
296
    tab = av_malloc(fft_size * sizeof(FFTComplex));
297
    tab1 = av_malloc(fft_size * sizeof(FFTComplex));
298
    tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
299
    tab2 = av_malloc(fft_size * sizeof(FFTSample));
300

    
301
    switch (transform) {
302
    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
        break;
310
    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
    case TRANSFORM_RDFT:
320
        if (do_inverse)
321
            av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
322
        else
323
            av_log(NULL, AV_LOG_INFO,"DFT_R2C");
324
        ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
325
        fft_ref_init(fft_nbits, do_inverse);
326
        break;
327
    case TRANSFORM_DCT:
328
        if (do_inverse)
329
            av_log(NULL, AV_LOG_INFO,"DCT_III");
330
        else
331
            av_log(NULL, AV_LOG_INFO,"DCT_II");
332
        ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);
333
        break;
334
#endif
335
    default:
336
        av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
337
        return 1;
338
    }
339
    av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
340

    
341
    /* generate random data */
342

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

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

    
351
    switch (transform) {
352
    case TRANSFORM_MDCT:
353
        if (do_inverse) {
354
            imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
355
            m->imdct_calc(m, tab2, (FFTSample *)tab1);
356
            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
        break;
365
    case TRANSFORM_FFT:
366
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
367
        s->fft_permute(s, tab);
368
        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
        if (do_inverse) {
376
            tab1[         0].im = 0;
377
            tab1[fft_size_2].im = 0;
378
            for (i = 1; i < fft_size_2; i++) {
379
                tab1[fft_size_2+i].re =  tab1[fft_size_2-i].re;
380
                tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;
381
            }
382

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

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

    
417
    /* do a speed test */
418

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

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

    
464
    switch (transform) {
465
    case TRANSFORM_MDCT:
466
        ff_mdct_end(m);
467
        break;
468
    case TRANSFORM_FFT:
469
        ff_fft_end(s);
470
        break;
471
#if CONFIG_FFT_FLOAT
472
    case TRANSFORM_RDFT:
473
        ff_rdft_end(r);
474
        break;
475
    case TRANSFORM_DCT:
476
        ff_dct_end(d);
477
        break;
478
#endif
479
    }
480

    
481
    av_free(tab);
482
    av_free(tab1);
483
    av_free(tab2);
484
    av_free(tab_ref);
485
    av_free(exptab);
486

    
487
    return err;
488
}