2 * (c) 2002 Fabrice Bellard
4 * This file is part of FFmpeg.
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.
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.
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
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"
47 #include "libavcodec/avfft.h"
49 #include "libavcodec/fft.h"
53 #include "libavcodec/dct.h"
54 #include "libavcodec/rdft.h"
59 #define MUL16(a, b) ((a) * (b))
61 #define CMAC(pre, pim, are, aim, bre, bim) \
63 pre += (MUL16(are, bre) - MUL16(aim, bim)); \
64 pim += (MUL16(are, bim) + MUL16(bre, aim)); \
67 #if FFT_FLOAT || AVFFT
69 #define REF_SCALE(x, bits) (x)
73 #define REF_SCALE(x, bits) (x)
77 #define REF_SCALE(x, bits) ((x) / (1 << (bits)))
85 static int fft_ref_init(int nbits, int inverse)
87 int i, n = 1 << nbits;
89 exptab = av_malloc_array((n / 2), sizeof(*exptab));
91 return AVERROR(ENOMEM);
93 for (i = 0; i < (n / 2); i++) {
94 double alpha = 2 * M_PI * (float) i / (float) n;
95 double c1 = cos(alpha), s1 = sin(alpha);
104 static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
110 for (i = 0; i < n; i++) {
111 double tmp_re = 0, tmp_im = 0;
113 for (j = 0; j < n; j++) {
115 int k = (i * j) & (n - 1);
117 c = -exptab[k - n2].re;
118 s = -exptab[k - n2].im;
123 CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
126 tabr[i].re = REF_SCALE(tmp_re, nbits);
127 tabr[i].im = REF_SCALE(tmp_im, nbits);
132 static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
134 int i, k, n = 1 << nbits;
136 for (i = 0; i < n; i++) {
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));
143 out[i] = REF_SCALE(-sum, nbits - 2);
147 /* NOTE: no normalisation by 1 / N is done */
148 static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
150 int i, k, n = 1 << nbits;
153 for (k = 0; k < n / 2; k++) {
155 for (i = 0; i < n; i++) {
156 double a = (2 * M_PI * (2 * i + 1 + n / 2) * (2 * k + 1) / (4 * n));
157 s += input[i] * cos(a);
159 output[k] = REF_SCALE(s, nbits - 1);
162 #endif /* CONFIG_MDCT */
166 static void idct_ref(FFTSample *output, FFTSample *input, int nbits)
168 int i, k, n = 1 << nbits;
171 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 s += input[k] * cos(a);
177 output[i] = 2 * s / n;
181 static void dct_ref(FFTSample *output, FFTSample *input, int nbits)
183 int i, k, n = 1 << nbits;
186 for (k = 0; k < n; k++) {
188 for (i = 0; i < n; i++) {
189 double a = M_PI * k * (i + 0.5) / n;
190 s += input[i] * cos(a);
195 #endif /* CONFIG_DCT */
196 #endif /* FFT_FLOAT */
198 static FFTSample frandom(AVLFG *prng)
200 return (int16_t) av_lfg_get(prng) / 32768.0 * RANGE;
203 static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
206 double error = 0, max = 0;
208 for (i = 0; i < n; i++) {
209 double e = fabs(tab1[i] - (tab2[i] / scale)) / RANGE;
211 av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
212 i, tab1[i], tab2[i]);
219 av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error / n));
223 static inline void fft_init(FFTContext **s, int nbits, int inverse)
226 *s = av_fft_init(nbits, inverse);
228 ff_fft_init(*s, nbits, inverse);
232 static inline void mdct_init(FFTContext **s, int nbits, int inverse, double scale)
235 *s = av_mdct_init(nbits, inverse, scale);
237 ff_mdct_init(*s, nbits, inverse, scale);
241 static inline void mdct_calc(FFTContext *s, FFTSample *output, const FFTSample *input)
244 av_mdct_calc(s, output, input);
246 s->mdct_calc(s, output, input);
250 static inline void imdct_calc(struct FFTContext *s, FFTSample *output, const FFTSample *input)
253 av_imdct_calc(s, output, input);
255 s->imdct_calc(s, output, input);
259 static inline void fft_permute(FFTContext *s, FFTComplex *z)
262 av_fft_permute(s, z);
264 s->fft_permute(s, z);
268 static inline void fft_calc(FFTContext *s, FFTComplex *z)
277 static inline void mdct_end(FFTContext *s)
286 static inline void fft_end(FFTContext *s)
296 static inline void rdft_init(RDFTContext **r, int nbits, enum RDFTransformType trans)
299 *r = av_rdft_init(nbits, trans);
301 ff_rdft_init(*r, nbits, trans);
305 static inline void dct_init(DCTContext **d, int nbits, enum DCTTransformType trans)
308 *d = av_dct_init(nbits, trans);
310 ff_dct_init(*d, nbits, trans);
314 static inline void rdft_calc(RDFTContext *r, FFTSample *tab)
317 av_rdft_calc(r, tab);
319 r->rdft_calc(r, tab);
323 static inline void dct_calc(DCTContext *d, FFTSample *data)
326 av_dct_calc(d, data);
328 d->dct_calc(d, data);
332 static inline void rdft_end(RDFTContext *r)
341 static inline void dct_end(DCTContext *d)
349 #endif /* FFT_FLOAT */
351 static void help(void)
353 av_log(NULL, AV_LOG_INFO,
354 "usage: fft-test [-h] [-s] [-i] [-n b]\n"
355 "-h print this help\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");
373 #include "compat/getopt.c"
376 int main(int argc, char **argv)
378 FFTComplex *tab, *tab1, *tab_ref;
380 enum tf_transform transform = TRANSFORM_FFT;
385 #endif /* FFT_FLOAT */
387 int do_speed = 0, do_inverse = 0;
388 int fft_nbits = 9, fft_size;
393 s = av_mallocz(sizeof(*s));
394 m = av_mallocz(sizeof(*m));
397 #if !AVFFT && FFT_FLOAT
398 r = av_mallocz(sizeof(*r));
399 d = av_mallocz(sizeof(*d));
402 av_lfg_init(&prng, 1);
405 int c = getopt(argc, argv, "hsimrdn:f:c:");
419 transform = TRANSFORM_MDCT;
422 transform = TRANSFORM_RDFT;
425 transform = TRANSFORM_DCT;
428 fft_nbits = atoi(optarg);
431 scale = atof(optarg);
435 unsigned cpuflags = av_get_cpu_flags();
437 if (av_parse_cpu_caps(&cpuflags, optarg) < 0)
440 av_force_cpu_flags(cpuflags);
446 fft_size = 1 << fft_nbits;
447 tab = av_malloc_array(fft_size, sizeof(FFTComplex));
448 tab1 = av_malloc_array(fft_size, sizeof(FFTComplex));
449 tab_ref = av_malloc_array(fft_size, sizeof(FFTComplex));
450 tab2 = av_malloc_array(fft_size, sizeof(FFTSample));
452 if (!(tab && tab1 && tab_ref && tab2))
458 av_log(NULL, AV_LOG_INFO, "Scale factor is set to %f\n", scale);
460 av_log(NULL, AV_LOG_INFO, "IMDCT");
462 av_log(NULL, AV_LOG_INFO, "MDCT");
463 mdct_init(&m, fft_nbits, do_inverse, scale);
465 #endif /* CONFIG_MDCT */
468 av_log(NULL, AV_LOG_INFO, "IFFT");
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)
479 av_log(NULL, AV_LOG_INFO, "IDFT_C2R");
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)
486 # endif /* CONFIG_RDFT */
490 av_log(NULL, AV_LOG_INFO, "DCT_III");
492 av_log(NULL, AV_LOG_INFO, "DCT_II");
493 dct_init(&d, fft_nbits, do_inverse ? DCT_III : DCT_II);
495 # endif /* CONFIG_DCT */
496 #endif /* FFT_FLOAT */
498 av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
501 av_log(NULL, AV_LOG_INFO, " %d test\n", fft_size);
503 /* generate random data */
505 for (i = 0; i < fft_size; i++) {
506 tab1[i].re = frandom(&prng);
507 tab1[i].im = frandom(&prng);
510 /* checking result */
511 av_log(NULL, AV_LOG_INFO, "Checking...\n");
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);
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);
526 #endif /* CONFIG_MDCT */
528 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
532 fft_ref(tab_ref, tab1, fft_nbits);
533 err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 1.0);
539 int fft_size_2 = fft_size >> 1;
542 tab1[fft_size_2].im = 0;
543 for (i = 1; i < fft_size_2; i++) {
544 tab1[fft_size_2 + i].re = tab1[fft_size_2 - i].re;
545 tab1[fft_size_2 + i].im = -tab1[fft_size_2 - i].im;
548 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
549 tab2[1] = tab1[fft_size_2].re;
552 fft_ref(tab_ref, tab1, fft_nbits);
553 for (i = 0; i < fft_size; i++) {
557 err = check_diff(&tab_ref->re, &tab->re, fft_size * 2, 0.5);
559 for (i = 0; i < fft_size; i++) {
560 tab2[i] = tab1[i].re;
564 fft_ref(tab_ref, tab1, fft_nbits);
565 tab_ref[0].im = tab_ref[fft_size_2].re;
566 err = check_diff(&tab_ref->re, tab2, fft_size, 1.0);
570 #endif /* CONFIG_RDFT */
573 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
574 dct_calc(d, &tab->re);
576 idct_ref(&tab_ref->re, &tab1->re, fft_nbits);
578 dct_ref(&tab_ref->re, &tab1->re, fft_nbits);
579 err = check_diff(&tab_ref->re, &tab->re, fft_size, 1.0);
581 #endif /* CONFIG_DCT */
582 #endif /* FFT_FLOAT */
585 /* do a speed test */
588 int64_t time_start, duration;
591 av_log(NULL, AV_LOG_INFO, "Speed test...\n");
592 /* we measure during about 1 seconds */
595 time_start = av_gettime_relative();
596 for (it = 0; it < nb_its; it++) {
600 imdct_calc(m, &tab->re, &tab1->re);
602 mdct_calc(m, &tab->re, &tab1->re);
605 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
610 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
614 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
617 #endif /* FFT_FLOAT */
620 duration = av_gettime_relative() - time_start;
621 if (duration >= 1000000)
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,
637 #endif /* CONFIG_MDCT */
646 # endif /* CONFIG_RDFT */
651 # endif /* CONFIG_DCT */
652 #endif /* FFT_FLOAT */
667 #if !AVFFT && FFT_FLOAT
673 printf("Error: %d.\n", err);