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