1/*
2 * (c) 2002 Fabrice Bellard
3 *
4 * This file is part of FFmpeg.
5 *
6 * FFmpeg 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 * FFmpeg 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 FFmpeg; 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/cpu.h"
27#include "libavutil/mathematics.h"
28#include "libavutil/lfg.h"
29#include "libavutil/log.h"
30#include "libavutil/time.h"
31#include "fft.h"
32#if FFT_FLOAT
33#include "dct.h"
34#include "rdft.h"
35#endif
36#include <math.h>
37#if HAVE_UNISTD_H
38#include <unistd.h>
39#endif
40#include <stdio.h>
41#include <stdlib.h>
42#include <string.h>
43
44/* reference fft */
45
46#define MUL16(a,b) ((a) * (b))
47
48#define CMAC(pre, pim, are, aim, bre, bim) \
49{\
50   pre += (MUL16(are, bre) - MUL16(aim, bim));\
51   pim += (MUL16(are, bim) + MUL16(bre, aim));\
52}
53
54#if FFT_FLOAT
55#   define RANGE 1.0
56#   define REF_SCALE(x, bits)  (x)
57#   define FMT "%10.6f"
58#elif FFT_FIXED_32
59#   define RANGE 8388608
60#   define REF_SCALE(x, bits) (x)
61#   define FMT "%6d"
62#else
63#   define RANGE 16384
64#   define REF_SCALE(x, bits) ((x) / (1<<(bits)))
65#   define FMT "%6d"
66#endif
67
68struct {
69    float re, im;
70} *exptab;
71
72static void fft_ref_init(int nbits, int inverse)
73{
74    int n, i;
75    double c1, s1, alpha;
76
77    n = 1 << nbits;
78    exptab = av_malloc_array((n / 2), sizeof(*exptab));
79
80    for (i = 0; i < (n/2); i++) {
81        alpha = 2 * M_PI * (float)i / (float)n;
82        c1 = cos(alpha);
83        s1 = sin(alpha);
84        if (!inverse)
85            s1 = -s1;
86        exptab[i].re = c1;
87        exptab[i].im = s1;
88    }
89}
90
91static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
92{
93    int n, i, j, k, n2;
94    double tmp_re, tmp_im, s, c;
95    FFTComplex *q;
96
97    n = 1 << nbits;
98    n2 = n >> 1;
99    for (i = 0; i < n; i++) {
100        tmp_re = 0;
101        tmp_im = 0;
102        q = tab;
103        for (j = 0; j < n; j++) {
104            k = (i * j) & (n - 1);
105            if (k >= n2) {
106                c = -exptab[k - n2].re;
107                s = -exptab[k - n2].im;
108            } else {
109                c = exptab[k].re;
110                s = exptab[k].im;
111            }
112            CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
113            q++;
114        }
115        tabr[i].re = REF_SCALE(tmp_re, nbits);
116        tabr[i].im = REF_SCALE(tmp_im, nbits);
117    }
118}
119
120#if CONFIG_MDCT
121static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
122{
123    int n = 1<<nbits;
124    int k, i, a;
125    double sum, f;
126
127    for (i = 0; i < n; i++) {
128        sum = 0;
129        for (k = 0; k < n/2; k++) {
130            a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
131            f = cos(M_PI * a / (double)(2 * n));
132            sum += f * in[k];
133        }
134        out[i] = REF_SCALE(-sum, nbits - 2);
135    }
136}
137
138/* NOTE: no normalisation by 1 / N is done */
139static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
140{
141    int n = 1<<nbits;
142    int k, i;
143    double a, s;
144
145    /* do it by hand */
146    for (k = 0; k < n/2; k++) {
147        s = 0;
148        for (i = 0; i < n; i++) {
149            a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
150            s += input[i] * cos(a);
151        }
152        output[k] = REF_SCALE(s, nbits - 1);
153    }
154}
155#endif /* CONFIG_MDCT */
156
157#if FFT_FLOAT
158#if CONFIG_DCT
159static void idct_ref(FFTSample *output, FFTSample *input, int nbits)
160{
161    int n = 1<<nbits;
162    int k, i;
163    double a, s;
164
165    /* do it by hand */
166    for (i = 0; i < n; i++) {
167        s = 0.5 * input[0];
168        for (k = 1; k < n; k++) {
169            a = M_PI*k*(i+0.5) / n;
170            s += input[k] * cos(a);
171        }
172        output[i] = 2 * s / n;
173    }
174}
175static void dct_ref(FFTSample *output, FFTSample *input, int nbits)
176{
177    int n = 1<<nbits;
178    int k, i;
179    double a, s;
180
181    /* do it by hand */
182    for (k = 0; k < n; k++) {
183        s = 0;
184        for (i = 0; i < n; i++) {
185            a = M_PI*k*(i+0.5) / n;
186            s += input[i] * cos(a);
187        }
188        output[k] = s;
189    }
190}
191#endif /* CONFIG_DCT */
192#endif
193
194
195static FFTSample frandom(AVLFG *prng)
196{
197    return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;
198}
199
200static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
201{
202    int i;
203    double max= 0;
204    double error= 0;
205    int err = 0;
206
207    for (i = 0; i < n; i++) {
208        double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;
209        if (e >= 1e-3) {
210            av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
211                   i, tab1[i], tab2[i]);
212            err = 1;
213        }
214        error+= e*e;
215        if(e>max) max= e;
216    }
217    av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error/n));
218    return err;
219}
220
221
222static void help(void)
223{
224    av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
225           "-h     print this help\n"
226           "-s     speed test\n"
227           "-m     (I)MDCT test\n"
228           "-d     (I)DCT test\n"
229           "-r     (I)RDFT test\n"
230           "-i     inverse transform test\n"
231           "-n b   set the transform size to 2^b\n"
232           "-f x   set scale factor for output data of (I)MDCT to x\n"
233           );
234}
235
236enum tf_transform {
237    TRANSFORM_FFT,
238    TRANSFORM_MDCT,
239    TRANSFORM_RDFT,
240    TRANSFORM_DCT,
241};
242
243#if !HAVE_GETOPT
244#include "compat/getopt.c"
245#endif
246
247int main(int argc, char **argv)
248{
249    FFTComplex *tab, *tab1, *tab_ref;
250    FFTSample *tab2;
251    int it, i, c;
252    int cpuflags;
253    int do_speed = 0;
254    int err = 1;
255    enum tf_transform transform = TRANSFORM_FFT;
256    int do_inverse = 0;
257    FFTContext s1, *s = &s1;
258    FFTContext m1, *m = &m1;
259#if FFT_FLOAT
260    RDFTContext r1, *r = &r1;
261    DCTContext d1, *d = &d1;
262    int fft_size_2;
263#endif
264    int fft_nbits, fft_size;
265    double scale = 1.0;
266    AVLFG prng;
267    av_lfg_init(&prng, 1);
268
269    fft_nbits = 9;
270    for(;;) {
271        c = getopt(argc, argv, "hsimrdn:f:c:");
272        if (c == -1)
273            break;
274        switch(c) {
275        case 'h':
276            help();
277            return 1;
278        case 's':
279            do_speed = 1;
280            break;
281        case 'i':
282            do_inverse = 1;
283            break;
284        case 'm':
285            transform = TRANSFORM_MDCT;
286            break;
287        case 'r':
288            transform = TRANSFORM_RDFT;
289            break;
290        case 'd':
291            transform = TRANSFORM_DCT;
292            break;
293        case 'n':
294            fft_nbits = atoi(optarg);
295            break;
296        case 'f':
297            scale = atof(optarg);
298            break;
299        case 'c':
300            cpuflags = av_get_cpu_flags();
301
302            if (av_parse_cpu_caps(&cpuflags, optarg) < 0)
303                return 1;
304
305            av_force_cpu_flags(cpuflags);
306            break;
307        }
308    }
309
310    fft_size = 1 << fft_nbits;
311    tab = av_malloc_array(fft_size, sizeof(FFTComplex));
312    tab1 = av_malloc_array(fft_size, sizeof(FFTComplex));
313    tab_ref = av_malloc_array(fft_size, sizeof(FFTComplex));
314    tab2 = av_malloc_array(fft_size, sizeof(FFTSample));
315
316    switch (transform) {
317#if CONFIG_MDCT
318    case TRANSFORM_MDCT:
319        av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale);
320        if (do_inverse)
321            av_log(NULL, AV_LOG_INFO,"IMDCT");
322        else
323            av_log(NULL, AV_LOG_INFO,"MDCT");
324        ff_mdct_init(m, fft_nbits, do_inverse, scale);
325        break;
326#endif /* CONFIG_MDCT */
327    case TRANSFORM_FFT:
328        if (do_inverse)
329            av_log(NULL, AV_LOG_INFO,"IFFT");
330        else
331            av_log(NULL, AV_LOG_INFO,"FFT");
332        ff_fft_init(s, fft_nbits, do_inverse);
333        fft_ref_init(fft_nbits, do_inverse);
334        break;
335#if FFT_FLOAT
336#    if CONFIG_RDFT
337    case TRANSFORM_RDFT:
338        if (do_inverse)
339            av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
340        else
341            av_log(NULL, AV_LOG_INFO,"DFT_R2C");
342        ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
343        fft_ref_init(fft_nbits, do_inverse);
344        break;
345#    endif /* CONFIG_RDFT */
346#    if CONFIG_DCT
347    case TRANSFORM_DCT:
348        if (do_inverse)
349            av_log(NULL, AV_LOG_INFO,"DCT_III");
350        else
351            av_log(NULL, AV_LOG_INFO,"DCT_II");
352        ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);
353        break;
354#    endif /* CONFIG_DCT */
355#endif
356    default:
357        av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
358        return 1;
359    }
360    av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
361
362    /* generate random data */
363
364    for (i = 0; i < fft_size; i++) {
365        tab1[i].re = frandom(&prng);
366        tab1[i].im = frandom(&prng);
367    }
368
369    /* checking result */
370    av_log(NULL, AV_LOG_INFO,"Checking...\n");
371
372    switch (transform) {
373#if CONFIG_MDCT
374    case TRANSFORM_MDCT:
375        if (do_inverse) {
376            imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
377            m->imdct_calc(m, tab2, (FFTSample *)tab1);
378            err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);
379        } else {
380            mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
381
382            m->mdct_calc(m, tab2, (FFTSample *)tab1);
383
384            err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);
385        }
386        break;
387#endif /* CONFIG_MDCT */
388    case TRANSFORM_FFT:
389        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
390        s->fft_permute(s, tab);
391        s->fft_calc(s, tab);
392
393        fft_ref(tab_ref, tab1, fft_nbits);
394        err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);
395        break;
396#if FFT_FLOAT
397#if CONFIG_RDFT
398    case TRANSFORM_RDFT:
399        fft_size_2 = fft_size >> 1;
400        if (do_inverse) {
401            tab1[         0].im = 0;
402            tab1[fft_size_2].im = 0;
403            for (i = 1; i < fft_size_2; i++) {
404                tab1[fft_size_2+i].re =  tab1[fft_size_2-i].re;
405                tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;
406            }
407
408            memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
409            tab2[1] = tab1[fft_size_2].re;
410
411            r->rdft_calc(r, tab2);
412            fft_ref(tab_ref, tab1, fft_nbits);
413            for (i = 0; i < fft_size; i++) {
414                tab[i].re = tab2[i];
415                tab[i].im = 0;
416            }
417            err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5);
418        } else {
419            for (i = 0; i < fft_size; i++) {
420                tab2[i]    = tab1[i].re;
421                tab1[i].im = 0;
422            }
423            r->rdft_calc(r, tab2);
424            fft_ref(tab_ref, tab1, fft_nbits);
425            tab_ref[0].im = tab_ref[fft_size_2].re;
426            err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);
427        }
428        break;
429#endif /* CONFIG_RDFT */
430#if CONFIG_DCT
431    case TRANSFORM_DCT:
432        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
433        d->dct_calc(d, (FFTSample *)tab);
434        if (do_inverse) {
435            idct_ref((FFTSample*)tab_ref, (FFTSample *)tab1, fft_nbits);
436        } else {
437            dct_ref((FFTSample*)tab_ref, (FFTSample *)tab1, fft_nbits);
438        }
439        err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);
440        break;
441#endif /* CONFIG_DCT */
442#endif
443    }
444
445    /* do a speed test */
446
447    if (do_speed) {
448        int64_t time_start, duration;
449        int nb_its;
450
451        av_log(NULL, AV_LOG_INFO,"Speed test...\n");
452        /* we measure during about 1 seconds */
453        nb_its = 1;
454        for(;;) {
455            time_start = av_gettime_relative();
456            for (it = 0; it < nb_its; it++) {
457                switch (transform) {
458                case TRANSFORM_MDCT:
459                    if (do_inverse) {
460                        m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
461                    } else {
462                        m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
463                    }
464                    break;
465                case TRANSFORM_FFT:
466                    memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
467                    s->fft_calc(s, tab);
468                    break;
469#if FFT_FLOAT
470                case TRANSFORM_RDFT:
471                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
472                    r->rdft_calc(r, tab2);
473                    break;
474                case TRANSFORM_DCT:
475                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
476                    d->dct_calc(d, tab2);
477                    break;
478#endif
479                }
480            }
481            duration = av_gettime_relative() - time_start;
482            if (duration >= 1000000)
483                break;
484            nb_its *= 2;
485        }
486        av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
487               (double)duration / nb_its,
488               (double)duration / 1000000.0,
489               nb_its);
490    }
491
492    switch (transform) {
493#if CONFIG_MDCT
494    case TRANSFORM_MDCT:
495        ff_mdct_end(m);
496        break;
497#endif /* CONFIG_MDCT */
498    case TRANSFORM_FFT:
499        ff_fft_end(s);
500        break;
501#if FFT_FLOAT
502#    if CONFIG_RDFT
503    case TRANSFORM_RDFT:
504        ff_rdft_end(r);
505        break;
506#    endif /* CONFIG_RDFT */
507#    if CONFIG_DCT
508    case TRANSFORM_DCT:
509        ff_dct_end(d);
510        break;
511#    endif /* CONFIG_DCT */
512#endif
513    }
514
515    av_free(tab);
516    av_free(tab1);
517    av_free(tab2);
518    av_free(tab_ref);
519    av_free(exptab);
520
521    if (err)
522        printf("Error: %d.\n", err);
523
524    return !!err;
525}
526