FFmpeg coverage

GCC Code Coverage Report

Directory:	../../../ffmpeg/		Exec	Total	Coverage
File:	src/libavcodec/tests/fft.c	Lines:	287	339	84.7 %
Date:	2021-04-14 23:45:22	Branches:	83	112	74.1 %


/*
 * (c) 2002 Fabrice Bellard
 *
 * This file is part of FFmpeg.
 *
 * FFmpeg is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * FFmpeg is distributed in the hope that it will be useful,
 * 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
 * License along with FFmpeg; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
 */

/**
 * @file
 * FFT and MDCT tests.
 */

#include "config.h"

#ifndef AVFFT
#define AVFFT 0
#endif

#include <math.h>
#if HAVE_UNISTD_H
#include <unistd.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "libavutil/cpu.h"
#include "libavutil/lfg.h"
#include "libavutil/log.h"
#include "libavutil/mathematics.h"
#include "libavutil/time.h"

#if AVFFT
#include "libavcodec/avfft.h"
#else
#include "libavcodec/fft.h"
#endif

#if FFT_FLOAT
#include "libavcodec/dct.h"
#include "libavcodec/rdft.h"
#endif

/* 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)); \
    }

#if FFT_FLOAT || AVFFT
#define RANGE 1.0
#define REF_SCALE(x, bits)  (x)
#define FMT "%10.6f"
#elif FFT_FIXED_32
#define RANGE 8388608
#define REF_SCALE(x, bits) (x)
#define FMT "%6d"
#else
#define RANGE 16384
#define REF_SCALE(x, bits) ((x) / (1 << (bits)))
#define FMT "%6d"
#endif

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

static int fft_ref_init(int nbits, int inverse)
{
    int i, n = 1 << nbits;

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

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

static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
{
    int i, j;
    int n  = 1 << nbits;
    int n2 = n >> 1;

    for (i = 0; i < n; i++) {
        double tmp_re = 0, tmp_im = 0;
        FFTComplex *q = tab;
        for (j = 0; j < n; j++) {
            double s, c;
            int 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 = REF_SCALE(tmp_re, nbits);
        tabr[i].im = REF_SCALE(tmp_im, nbits);
    }
}

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

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

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

    /* do it by hand */
    for (k = 0; k < n / 2; k++) {
        double s = 0;
        for (i = 0; i < n; i++) {
            double a = (2 * M_PI * (2 * i + 1 + n / 2) * (2 * k + 1) / (4 * n));
            s += input[i] * cos(a);
        }
        output[k] = REF_SCALE(s, nbits - 1);
    }
}
#endif /* CONFIG_MDCT */

#if FFT_FLOAT
#if CONFIG_DCT
static void idct_ref(FFTSample *output, FFTSample *input, int nbits)
{
    int i, k, n = 1 << nbits;

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

static void dct_ref(FFTSample *output, FFTSample *input, int nbits)
{
    int i, k, n = 1 << nbits;

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

static FFTSample frandom(AVLFG *prng)
{
    return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE;
}

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

    for (i = 0; i < n; i++) {
        double e = fabs(tab1[i] - (tab2[i] / scale)) / RANGE;
        if (e >= 1e-3) {
            av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
                   i, tab1[i], tab2[i]);
            err = 1;
        }
        error += e * e;
        if (e > max)
            max = e;
    }
    av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error / n));
    return err;
}

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 */

static void help(void)
{
    av_log(NULL, AV_LOG_INFO,
           "usage: fft-test [-h] [-s] [-i] [-n b]\n"
           "-h     print this help\n"
           "-s     speed test\n"
           "-m     (I)MDCT test\n"
           "-d     (I)DCT test\n"
           "-r     (I)RDFT test\n"
           "-i     inverse transform test\n"
           "-n b   set the transform size to 2^b\n"
           "-f x   set scale factor for output data of (I)MDCT to x\n");
}

enum tf_transform {
    TRANSFORM_FFT,
    TRANSFORM_MDCT,
    TRANSFORM_RDFT,
    TRANSFORM_DCT,
};

#if !HAVE_GETOPT
#include "compat/getopt.c"
#endif

int main(int argc, char **argv)
{
    FFTComplex *tab, *tab1, *tab_ref;
    FFTSample *tab2;
    enum tf_transform transform = TRANSFORM_FFT;
    FFTContext *m, *s;
#if FFT_FLOAT
    RDFTContext *r;
    DCTContext *d;
#endif /* FFT_FLOAT */
    int it, i, err = 1;
    int do_speed = 0, do_inverse = 0;
    int fft_nbits = 9, fft_size;
    double scale = 1.0;
    AVLFG prng;

#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

    av_lfg_init(&prng, 1);

    for (;;) {
        int c = getopt(argc, argv, "hsimrdn:f:c:");
        if (c == -1)
            break;
        switch (c) {
        case 'h':
            help();
            return 1;
        case 's':
            do_speed = 1;
            break;
        case 'i':
            do_inverse = 1;
            break;
        case 'm':
            transform = TRANSFORM_MDCT;
            break;
        case 'r':
            transform = TRANSFORM_RDFT;
            break;
        case 'd':
            transform = TRANSFORM_DCT;
            break;
        case 'n':
            fft_nbits = atoi(optarg);
            break;
        case 'f':
            scale = atof(optarg);
            break;
        case 'c':
        {
            unsigned cpuflags = av_get_cpu_flags();

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

            av_force_cpu_flags(cpuflags);
            break;
        }
        }
    }

    fft_size = 1 << fft_nbits;
    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));

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

    switch (transform) {
#if CONFIG_MDCT
    case TRANSFORM_MDCT:
        av_log(NULL, AV_LOG_INFO, "Scale factor is set to %f\n", scale);
        if (do_inverse)
            av_log(NULL, AV_LOG_INFO, "IMDCT");
        else
            av_log(NULL, AV_LOG_INFO, "MDCT");
        mdct_init(&m, fft_nbits, do_inverse, scale);
        break;
#endif /* CONFIG_MDCT */
    case TRANSFORM_FFT:
        if (do_inverse)
            av_log(NULL, AV_LOG_INFO, "IFFT");
        else
            av_log(NULL, AV_LOG_INFO, "FFT");
        fft_init(&s, fft_nbits, do_inverse);
        if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
            goto cleanup;
        break;
#if FFT_FLOAT
#    if CONFIG_RDFT
    case TRANSFORM_RDFT:
        if (do_inverse)
            av_log(NULL, AV_LOG_INFO, "IDFT_C2R");
        else
            av_log(NULL, AV_LOG_INFO, "DFT_R2C");
        rdft_init(&r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
        if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
            goto cleanup;
        break;
#    endif /* CONFIG_RDFT */
#    if CONFIG_DCT
    case TRANSFORM_DCT:
        if (do_inverse)
            av_log(NULL, AV_LOG_INFO, "DCT_III");
        else
            av_log(NULL, AV_LOG_INFO, "DCT_II");
        dct_init(&d, fft_nbits, do_inverse ? DCT_III : DCT_II);
        break;
#    endif /* CONFIG_DCT */
#endif /* FFT_FLOAT */
    default:
        av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
        goto cleanup;
    }
    av_log(NULL, AV_LOG_INFO, " %d test\n", fft_size);

    /* generate random data */

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

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

    switch (transform) {
#if CONFIG_MDCT
    case TRANSFORM_MDCT:
        if (do_inverse) {
            imdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
            imdct_calc(m, tab2, &tab1->re);
            err = check_diff(&tab_ref->re, tab2, fft_size, scale);
        } else {
            mdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
            mdct_calc(m, tab2, &tab1->re);
            err = check_diff(&tab_ref->re, tab2, fft_size / 2, scale);
        }
        break;
#endif /* CONFIG_MDCT */
    case TRANSFORM_FFT:
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
        fft_permute(s, tab);
        fft_calc(s, tab);

        fft_ref(tab_ref, tab1, fft_nbits);
        err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 1.0);
        break;
#if FFT_FLOAT
#if CONFIG_RDFT
    case TRANSFORM_RDFT:
    {
        int fft_size_2 = fft_size >> 1;
        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;

            rdft_calc(r, tab2);
            fft_ref(tab_ref, tab1, fft_nbits);
            for (i = 0; i < fft_size; i++) {
                tab[i].re = tab2[i];
                tab[i].im = 0;
            }
            err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 0.5);
        } else {
            for (i = 0; i < fft_size; i++) {
                tab2[i]    = tab1[i].re;
                tab1[i].im = 0;
            }
            rdft_calc(r, tab2);
            fft_ref(tab_ref, tab1, fft_nbits);
            tab_ref[0].im = tab_ref[fft_size_2].re;
            err = check_diff(&tab_ref->re, tab2, fft_size, 1.0);
        }
        break;
    }
#endif /* CONFIG_RDFT */
#if CONFIG_DCT
    case TRANSFORM_DCT:
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
        dct_calc(d, &tab->re);
        if (do_inverse)
            idct_ref(&tab_ref->re, &tab1->re, fft_nbits);
        else
            dct_ref(&tab_ref->re, &tab1->re, fft_nbits);
        err = check_diff(&tab_ref->re, &tab->re, fft_size, 1.0);
        break;
#endif /* CONFIG_DCT */
#endif /* FFT_FLOAT */
    }

    /* do a speed test */

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

        av_log(NULL, AV_LOG_INFO, "Speed test...\n");
        /* we measure during about 1 seconds */
        nb_its = 1;
        for (;;) {
            time_start = av_gettime_relative();
            for (it = 0; it < nb_its; it++) {
                switch (transform) {
                case TRANSFORM_MDCT:
                    if (do_inverse)
                        imdct_calc(m, &tab->re, &tab1->re);
                    else
                        mdct_calc(m, &tab->re, &tab1->re);
                    break;
                case TRANSFORM_FFT:
                    memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
                    fft_calc(s, tab);
                    break;
#if FFT_FLOAT
                case TRANSFORM_RDFT:
                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
                    rdft_calc(r, tab2);
                    break;
                case TRANSFORM_DCT:
                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
                    dct_calc(d, tab2);
                    break;
#endif /* FFT_FLOAT */
                }
            }
            duration = av_gettime_relative() - time_start;
            if (duration >= 1000000)
                break;
            nb_its *= 2;
        }
        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,
               nb_its);
    }

    switch (transform) {
#if CONFIG_MDCT
    case TRANSFORM_MDCT:
        mdct_end(m);
        break;
#endif /* CONFIG_MDCT */
    case TRANSFORM_FFT:
        fft_end(s);
        break;
#if FFT_FLOAT
#    if CONFIG_RDFT
    case TRANSFORM_RDFT:
        rdft_end(r);
        break;
#    endif /* CONFIG_RDFT */
#    if CONFIG_DCT
    case TRANSFORM_DCT:
        dct_end(d);
        break;
#    endif /* CONFIG_DCT */
#endif /* FFT_FLOAT */
    }

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

#if !AVFFT
    av_free(s);
    av_free(m);
#endif

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

    if (err)
        printf("Error: %d.\n", err);

    return !!err;
}


Generated by: GCOVR (Version 4.2)

Line	Branch	Exec	Source
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 "config.h"
27
28			#ifndef AVFFT
29			#define AVFFT 0
30			#endif
31
32			#include <math.h>
33			#if HAVE_UNISTD_H
34			#include <unistd.h>
35			#endif
36			#include <stdio.h>
37			#include <stdlib.h>
38			#include <string.h>
39
40			#include "libavutil/cpu.h"
41			#include "libavutil/lfg.h"
42			#include "libavutil/log.h"
43			#include "libavutil/mathematics.h"
44			#include "libavutil/time.h"
45
46			#if AVFFT
47			#include "libavcodec/avfft.h"
48			#else
49			#include "libavcodec/fft.h"
50			#endif
51
52			#if FFT_FLOAT
53			#include "libavcodec/dct.h"
54			#include "libavcodec/rdft.h"
55			#endif
56
57			/* reference fft */
58
59			#define MUL16(a, b) ((a) * (b))
60
61			#define CMAC(pre, pim, are, aim, bre, bim) \
62			{ \
63			pre += (MUL16(are, bre) - MUL16(aim, bim)); \
64			pim += (MUL16(are, bim) + MUL16(bre, aim)); \
65			}
66
67			#if FFT_FLOAT \|\| AVFFT
68			#define RANGE 1.0
69			#define REF_SCALE(x, bits) (x)
70			#define FMT "%10.6f"
71			#elif FFT_FIXED_32
72			#define RANGE 8388608
73			#define REF_SCALE(x, bits) (x)
74			#define FMT "%6d"
75			#else
76			#define RANGE 16384
77			#define REF_SCALE(x, bits) ((x) / (1 << (bits)))
78			#define FMT "%6d"
79			#endif
80
81			static struct {
82			float re, im;
83			} *exptab;
84
85		99	static int fft_ref_init(int nbits, int inverse)
86			{
87		99	int i, n = 1 << nbits;
88
89		99	exptab = av_malloc_array((n / 2), sizeof(*exptab));
90	✗✓	99	if (!exptab)
91			return AVERROR(ENOMEM);
92
93	✓✓	43283	for (i = 0; i < (n / 2); i++) {
94		43184	double alpha = 2 * M_PI * (float) i / (float) n;
95		43184	double c1 = cos(alpha), s1 = sin(alpha);
96	✓✓	43184	if (!inverse)
97		22744	s1 = -s1;
98		43184	exptab[i].re = c1;
99		43184	exptab[i].im = s1;
100			}
101		99	return 0;
102			}
103
104		99	static void fft_ref(FFTComplex tabr, FFTComplex tab, int nbits)
105			{
106			int i, j;
107		99	int n = 1 << nbits;
108		99	int n2 = n >> 1;
109
110	✓✓	86467	for (i = 0; i < n; i++) {
111		86368	double tmp_re = 0, tmp_im = 0;
112		86368	FFTComplex *q = tab;
113	✓✓	226141024	for (j = 0; j < n; j++) {
114			double s, c;
115		226054656	int k = (i * j) & (n - 1);
116	✓✓	226054656	if (k >= n2) {
117		112984144	c = -exptab[k - n2].re;
118		112984144	s = -exptab[k - n2].im;
119			} else {
120		113070512	c = exptab[k].re;
121		113070512	s = exptab[k].im;
122			}
123		226054656	CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
124		226054656	q++;
125			}
126		86368	tabr[i].re = REF_SCALE(tmp_re, nbits);
127		86368	tabr[i].im = REF_SCALE(tmp_im, nbits);
128			}
129		99	}
130
131			#if CONFIG_MDCT
132		27	static void imdct_ref(FFTSample out, FFTSample in, int nbits)
133			{
134		27	int i, k, n = 1 << nbits;
135
136	✓✓	24555	for (i = 0; i < n; i++) {
137		24528	double sum = 0;
138	✓✓	33578832	for (k = 0; k < n / 2; k++) {
139		33554304	int a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
140		33554304	double f = cos(M_PI * a / (double) (2 * n));
141		33554304	sum += f * in[k];
142			}
143		24528	out[i] = REF_SCALE(-sum, nbits - 2);
144			}
145		27	}
146
147			/* NOTE: no normalisation by 1 / N is done */
148		18	static void mdct_ref(FFTSample output, FFTSample input, int nbits)
149			{
150		18	int i, k, n = 1 << nbits;
151
152			/* do it by hand */
153	✓✓	8194	for (k = 0; k < n / 2; k++) {
154		8176	double s = 0;
155	✓✓	22377712	for (i = 0; i < n; i++) {
156		22369536	double a = (2 * M_PI * (2 * i + 1 + n / 2) * (2 * k + 1) / (4 * n));
157		22369536	s += input[i] * cos(a);
158			}
159		8176	output[k] = REF_SCALE(s, nbits - 1);
160			}
161		18	}
162			#endif /* CONFIG_MDCT */
163
164			#if FFT_FLOAT
165			#if CONFIG_DCT
166		18	static void idct_ref(FFTSample output, FFTSample input, int nbits)
167			{
168		18	int i, k, n = 1 << nbits;
169
170			/* do it by hand */
171	✓✓	16370	for (i = 0; i < n; i++) {
172		16352	double s = 0.5 * input[0];
173	✓✓	44739072	for (k = 1; k < n; k++) {
174		44722720	double a = M_PI * k * (i + 0.5) / n;
175		44722720	s += input[k] * cos(a);
176			}
177		16352	output[i] = 2 * s / n;
178			}
179		18	}
180
181		18	static void dct_ref(FFTSample output, FFTSample input, int nbits)
182			{
183		18	int i, k, n = 1 << nbits;
184
185			/* do it by hand */
186	✓✓	16370	for (k = 0; k < n; k++) {
187		16352	double s = 0;
188	✓✓	44755424	for (i = 0; i < n; i++) {
189		44739072	double a = M_PI * k * (i + 0.5) / n;
190		44739072	s += input[i] * cos(a);
191			}
192		16352	output[k] = s;
193			}
194		18	}
195			#endif /* CONFIG_DCT */
196			#endif /* FFT_FLOAT */
197
198		319904	static FFTSample frandom(AVLFG *prng)
199			{
200		319904	return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE;
201			}
202
203		180	static int check_diff(FFTSample tab1, FFTSample tab2, int n, double scale)
204			{
205		180	int i, err = 0;
206		180	double error = 0, max = 0;
207
208	✓✓	221972	for (i = 0; i < n; i++) {
209		221792	double e = fabs(tab1[i] - (tab2[i] / scale)) / RANGE;
210	✗✓	221792	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		221792	error += e * e;
216	✓✓	221792	if (e > max)
217		834	max = e;
218			}
219		180	av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error / n));
220		180	return err;
221			}
222
223		63	static inline void fft_init(FFTContext **s, int nbits, int inverse)
224			{
225			#if AVFFT
226		18	*s = av_fft_init(nbits, inverse);
227			#else
228		45	ff_fft_init(*s, nbits, inverse);
229			#endif
230		63	}
231
232		45	static inline void mdct_init(FFTContext **s, int nbits, int inverse, double scale)
233			{
234			#if AVFFT
235		18	*s = av_mdct_init(nbits, inverse, scale);
236			#else
237		27	ff_mdct_init(*s, nbits, inverse, scale);
238			#endif
239		45	}
240
241		18	static inline void mdct_calc(FFTContext s, FFTSample output, const FFTSample *input)
242			{
243			#if AVFFT
244		9	av_mdct_calc(s, output, input);
245			#else
246		9	s->mdct_calc(s, output, input);
247			#endif
248		18	}
249
250		27	static inline void imdct_calc(struct FFTContext s, FFTSample output, const FFTSample *input)
251			{
252			#if AVFFT
253		9	av_imdct_calc(s, output, input);
254			#else
255		18	s->imdct_calc(s, output, input);
256			#endif
257		27	}
258
259		63	static inline void fft_permute(FFTContext s, FFTComplex z)
260			{
261			#if AVFFT
262		18	av_fft_permute(s, z);
263			#else
264		45	s->fft_permute(s, z);
265			#endif
266		63	}
267
268		63	static inline void fft_calc(FFTContext s, FFTComplex z)
269			{
270			#if AVFFT
271		18	av_fft_calc(s, z);
272			#else
273		45	s->fft_calc(s, z);
274			#endif
275		63	}
276
277		45	static inline void mdct_end(FFTContext *s)
278			{
279			#if AVFFT
280		18	av_mdct_end(s);
281			#else
282		27	ff_mdct_end(s);
283			#endif
284		45	}
285
286		63	static inline void fft_end(FFTContext *s)
287			{
288			#if AVFFT
289		18	av_fft_end(s);
290			#else
291		45	ff_fft_end(s);
292			#endif
293		63	}
294
295			#if FFT_FLOAT
296		36	static inline void rdft_init(RDFTContext **r, int nbits, enum RDFTransformType trans)
297			{
298			#if AVFFT
299		18	*r = av_rdft_init(nbits, trans);
300			#else
301		18	ff_rdft_init(*r, nbits, trans);
302			#endif
303		36	}
304
305		36	static inline void dct_init(DCTContext **d, int nbits, enum DCTTransformType trans)
306			{
307			#if AVFFT
308		18	*d = av_dct_init(nbits, trans);
309			#else
310		18	ff_dct_init(*d, nbits, trans);
311			#endif
312		36	}
313
314		36	static inline void rdft_calc(RDFTContext r, FFTSample tab)
315			{
316			#if AVFFT
317		18	av_rdft_calc(r, tab);
318			#else
319		18	r->rdft_calc(r, tab);
320			#endif
321		36	}
322
323		36	static inline void dct_calc(DCTContext d, FFTSample data)
324			{
325			#if AVFFT
326		18	av_dct_calc(d, data);
327			#else
328		18	d->dct_calc(d, data);
329			#endif
330		36	}
331
332		36	static inline void rdft_end(RDFTContext *r)
333			{
334			#if AVFFT
335		18	av_rdft_end(r);
336			#else
337		18	ff_rdft_end(r);
338			#endif
339		36	}
340
341		36	static inline void dct_end(DCTContext *d)
342			{
343			#if AVFFT
344		18	av_dct_end(d);
345			#else
346		18	ff_dct_end(d);
347			#endif
348		36	}
349			#endif /* FFT_FLOAT */
350
351			static void help(void)
352			{
353			av_log(NULL, AV_LOG_INFO,
354			"usage: fft-test [-h] [-s] [-i] [-n b]\n"
355			"-h print this help\n"
356			"-s speed test\n"
357			"-m (I)MDCT test\n"
358			"-d (I)DCT test\n"
359			"-r (I)RDFT test\n"
360			"-i inverse transform test\n"
361			"-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			enum tf_transform {
366			TRANSFORM_FFT,
367			TRANSFORM_MDCT,
368			TRANSFORM_RDFT,
369			TRANSFORM_DCT,
370			};
371
372			#if !HAVE_GETOPT
373			#include "compat/getopt.c"
374			#endif
375
376		180	int main(int argc, char **argv)
377			{
378			FFTComplex tab, tab1, *tab_ref;
379			FFTSample *tab2;
380		180	enum tf_transform transform = TRANSFORM_FFT;
381			FFTContext m, s;
382			#if FFT_FLOAT
383			RDFTContext *r;
384			DCTContext *d;
385			#endif /* FFT_FLOAT */
386		180	int it, i, err = 1;
387		180	int do_speed = 0, do_inverse = 0;
388		180	int fft_nbits = 9, fft_size;
389		180	double scale = 1.0;
390			AVLFG prng;
391
392			#if !AVFFT
393		108	s = av_mallocz(sizeof(*s));
394		108	m = av_mallocz(sizeof(*m));
395			#endif
396
397			#if !AVFFT && FFT_FLOAT
398		72	r = av_mallocz(sizeof(*r));
399		72	d = av_mallocz(sizeof(*d));
400			#endif
401
402		180	av_lfg_init(&prng, 1);
403
404		558	for (;;) {
405		738	int c = getopt(argc, argv, "hsimrdn:f:c:");
406	✓✓	738	if (c == -1)
407		180	break;
408	✗✗✓✓ ✓✓✓✗ ✓✗	558	switch (c) {
409			case 'h':
410			help();
411			return 1;
412			case 's':
413			do_speed = 1;
414			break;
415		90	case 'i':
416		90	do_inverse = 1;
417		90	break;
418		45	case 'm':
419		45	transform = TRANSFORM_MDCT;
420		45	break;
421		36	case 'r':
422		36	transform = TRANSFORM_RDFT;
423		36	break;
424		36	case 'd':
425		36	transform = TRANSFORM_DCT;
426		36	break;
427		171	case 'n':
428		171	fft_nbits = atoi(optarg);
429		171	break;
430			case 'f':
431			scale = atof(optarg);
432			break;
433		180	case 'c':
434			{
435		180	unsigned cpuflags = av_get_cpu_flags();
436
437	✗✓	180	if (av_parse_cpu_caps(&cpuflags, optarg) < 0)
438			return 1;
439
440		180	av_force_cpu_flags(cpuflags);
441		180	break;
442			}
443			}
444		558	}
445
446		180	fft_size = 1 << fft_nbits;
447		180	tab = av_malloc_array(fft_size, sizeof(FFTComplex));
448		180	tab1 = av_malloc_array(fft_size, sizeof(FFTComplex));
449		180	tab_ref = av_malloc_array(fft_size, sizeof(FFTComplex));
450		180	tab2 = av_malloc_array(fft_size, sizeof(FFTSample));
451
452	✓✗✓✗ ✓✗✗✓	180	if (!(tab && tab1 && tab_ref && tab2))
453			goto cleanup;
454
455	✓✓✓✓ ✗	180	switch (transform) {
456			#if CONFIG_MDCT
457		45	case TRANSFORM_MDCT:
458		45	av_log(NULL, AV_LOG_INFO, "Scale factor is set to %f\n", scale);
459	✓✓	45	if (do_inverse)
460		27	av_log(NULL, AV_LOG_INFO, "IMDCT");
461			else
462		18	av_log(NULL, AV_LOG_INFO, "MDCT");
463		45	mdct_init(&m, fft_nbits, do_inverse, scale);
464		45	break;
465			#endif /* CONFIG_MDCT */
466		63	case TRANSFORM_FFT:
467	✓✓	63	if (do_inverse)
468		27	av_log(NULL, AV_LOG_INFO, "IFFT");
469			else
470		36	av_log(NULL, AV_LOG_INFO, "FFT");
471		63	fft_init(&s, fft_nbits, do_inverse);
472	✗✓	63	if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
473			goto cleanup;
474		63	break;
475			#if FFT_FLOAT
476			# if CONFIG_RDFT
477		36	case TRANSFORM_RDFT:
478	✓✓	36	if (do_inverse)
479		18	av_log(NULL, AV_LOG_INFO, "IDFT_C2R");
480			else
481		18	av_log(NULL, AV_LOG_INFO, "DFT_R2C");
482		36	rdft_init(&r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
483	✗✓	36	if ((err = fft_ref_init(fft_nbits, do_inverse)) < 0)
484			goto cleanup;
485		36	break;
486			# endif /* CONFIG_RDFT */
487			# if CONFIG_DCT
488		36	case TRANSFORM_DCT:
489	✓✓	36	if (do_inverse)
490		18	av_log(NULL, AV_LOG_INFO, "DCT_III");
491			else
492		18	av_log(NULL, AV_LOG_INFO, "DCT_II");
493		36	dct_init(&d, fft_nbits, do_inverse ? DCT_III : DCT_II);
494		36	break;
495			# endif /* CONFIG_DCT */
496			#endif /* FFT_FLOAT */
497			default:
498			av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
499			goto cleanup;
500			}
501		180	av_log(NULL, AV_LOG_INFO, " %d test\n", fft_size);
502
503			/* generate random data */
504
505	✓✓	160132	for (i = 0; i < fft_size; i++) {
506		159952	tab1[i].re = frandom(&prng);
507		159952	tab1[i].im = frandom(&prng);
508			}
509
510			/* checking result */
511		180	av_log(NULL, AV_LOG_INFO, "Checking...\n");
512
513	✓✓✓✓ ✗	180	switch (transform) {
514			#if CONFIG_MDCT
515		45	case TRANSFORM_MDCT:
516	✓✓	45	if (do_inverse) {
517		27	imdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
518		27	imdct_calc(m, tab2, &tab1->re);
519		27	err = check_diff(&tab_ref->re, tab2, fft_size, scale);
520			} else {
521		18	mdct_ref(&tab_ref->re, &tab1->re, fft_nbits);
522		18	mdct_calc(m, tab2, &tab1->re);
523		18	err = check_diff(&tab_ref->re, tab2, fft_size / 2, scale);
524			}
525		45	break;
526			#endif /* CONFIG_MDCT */
527		63	case TRANSFORM_FFT:
528		63	memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
529		63	fft_permute(s, tab);
530		63	fft_calc(s, tab);
531
532		63	fft_ref(tab_ref, tab1, fft_nbits);
533		63	err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 1.0);
534		63	break;
535			#if FFT_FLOAT
536			#if CONFIG_RDFT
537		36	case TRANSFORM_RDFT:
538			{
539		36	int fft_size_2 = fft_size >> 1;
540	✓✓	36	if (do_inverse) {
541		18	tab1[0].im = 0;
542		18	tab1[fft_size_2].im = 0;
543	✓✓	8176	for (i = 1; i < fft_size_2; i++) {
544		8158	tab1[fft_size_2 + i].re = tab1[fft_size_2 - i].re;
545		8158	tab1[fft_size_2 + i].im = -tab1[fft_size_2 - i].im;
546			}
547
548		18	memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
549		18	tab2[1] = tab1[fft_size_2].re;
550
551		18	rdft_calc(r, tab2);
552		18	fft_ref(tab_ref, tab1, fft_nbits);
553	✓✓	16370	for (i = 0; i < fft_size; i++) {
554		16352	tab[i].re = tab2[i];
555		16352	tab[i].im = 0;
556			}
557		18	err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 0.5);
558			} else {
559	✓✓	16370	for (i = 0; i < fft_size; i++) {
560		16352	tab2[i] = tab1[i].re;
561		16352	tab1[i].im = 0;
562			}
563		18	rdft_calc(r, tab2);
564		18	fft_ref(tab_ref, tab1, fft_nbits);
565		18	tab_ref[0].im = tab_ref[fft_size_2].re;
566		18	err = check_diff(&tab_ref->re, tab2, fft_size, 1.0);
567			}
568		36	break;
569			}
570			#endif /* CONFIG_RDFT */
571			#if CONFIG_DCT
572		36	case TRANSFORM_DCT:
573		36	memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
574		36	dct_calc(d, &tab->re);
575	✓✓	36	if (do_inverse)
576		18	idct_ref(&tab_ref->re, &tab1->re, fft_nbits);
577			else
578		18	dct_ref(&tab_ref->re, &tab1->re, fft_nbits);
579		36	err = check_diff(&tab_ref->re, &tab->re, fft_size, 1.0);
580		36	break;
581			#endif /* CONFIG_DCT */
582			#endif /* FFT_FLOAT */
583			}
584
585			/* do a speed test */
586
587	✗✓	180	if (do_speed) {
588			int64_t time_start, duration;
589			int nb_its;
590
591			av_log(NULL, AV_LOG_INFO, "Speed test...\n");
592			/* we measure during about 1 seconds */
593			nb_its = 1;
594			for (;;) {
595			time_start = av_gettime_relative();
596			for (it = 0; it < nb_its; it++) {
597			switch (transform) {
598			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			break;
604			case TRANSFORM_FFT:
605			memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
606			fft_calc(s, tab);
607			break;
608			#if FFT_FLOAT
609			case TRANSFORM_RDFT:
610			memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
611			rdft_calc(r, tab2);
612			break;
613			case TRANSFORM_DCT:
614			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			if (duration >= 1000000)
622			break;
623			nb_its *= 2;
624			}
625			av_log(NULL, AV_LOG_INFO,
626			"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
627			(double) duration / nb_its,
628			(double) duration / 1000000.0,
629			nb_its);
630			}
631
632	✓✓✓✓ ✗	180	switch (transform) {
633			#if CONFIG_MDCT
634		45	case TRANSFORM_MDCT:
635		45	mdct_end(m);
636		45	break;
637			#endif /* CONFIG_MDCT */
638		63	case TRANSFORM_FFT:
639		63	fft_end(s);
640		63	break;
641			#if FFT_FLOAT
642			# if CONFIG_RDFT
643		36	case TRANSFORM_RDFT:
644		36	rdft_end(r);
645		36	break;
646			# endif /* CONFIG_RDFT */
647			# if CONFIG_DCT
648		36	case TRANSFORM_DCT:
649		36	dct_end(d);
650		36	break;
651			# endif /* CONFIG_DCT */
652			#endif /* FFT_FLOAT */
653			}
654
655		180	cleanup:
656		180	av_free(tab);
657		180	av_free(tab1);
658		180	av_free(tab2);
659		180	av_free(tab_ref);
660		180	av_free(exptab);
661
662			#if !AVFFT
663		108	av_free(s);
664		108	av_free(m);
665			#endif
666
667			#if !AVFFT && FFT_FLOAT
668		72	av_free(r);
669		72	av_free(d);
670			#endif
671
672	✗✓	180	if (err)
673			printf("Error: %d.\n", err);
674
675		180	return !!err;
676			}