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
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"
50 #define MUL16(a, b) ((a) * (b))
52 #define CMAC(pre, pim, are, aim, bre, bim) \
54 pre += (MUL16(are, bre) - MUL16(aim, bim)); \
55 pim += (MUL16(are, bim) + MUL16(bre, aim)); \
60 #define REF_SCALE(x, bits) (x)
64 #define REF_SCALE(x, bits) (x)
68 #define REF_SCALE(x, bits) ((x) / (1 << (bits)))
76 static int fft_ref_init(int nbits
, int inverse
)
78 int i
, n
= 1 << nbits
;
80 exptab
= av_malloc_array((n
/ 2), sizeof(*exptab
));
82 return AVERROR(ENOMEM
);
84 for (i
= 0; i
< (n
/ 2); i
++) {
85 double alpha
= 2 * M_PI
* (float) i
/ (float) n
;
86 double c1
= cos(alpha
), s1
= sin(alpha
);
95 static void fft_ref(FFTComplex
*tabr
, FFTComplex
*tab
, int nbits
)
101 for (i
= 0; i
< n
; i
++) {
102 double tmp_re
= 0, tmp_im
= 0;
104 for (j
= 0; j
< n
; j
++) {
106 int k
= (i
* j
) & (n
- 1);
108 c
= -exptab
[k
- n2
].re
;
109 s
= -exptab
[k
- n2
].im
;
114 CMAC(tmp_re
, tmp_im
, c
, s
, q
->re
, q
->im
);
117 tabr
[i
].re
= REF_SCALE(tmp_re
, nbits
);
118 tabr
[i
].im
= REF_SCALE(tmp_im
, nbits
);
123 static void imdct_ref(FFTSample
*out
, FFTSample
*in
, int nbits
)
125 int i
, k
, n
= 1 << nbits
;
127 for (i
= 0; i
< n
; i
++) {
129 for (k
= 0; k
< n
/ 2; k
++) {
130 int a
= (2 * i
+ 1 + (n
/ 2)) * (2 * k
+ 1);
131 double f
= cos(M_PI
* a
/ (double) (2 * n
));
134 out
[i
] = REF_SCALE(-sum
, nbits
- 2);
138 /* NOTE: no normalisation by 1 / N is done */
139 static void mdct_ref(FFTSample
*output
, FFTSample
*input
, int nbits
)
141 int i
, k
, n
= 1 << nbits
;
144 for (k
= 0; k
< n
/ 2; k
++) {
146 for (i
= 0; i
< n
; i
++) {
147 double a
= (2 * M_PI
* (2 * i
+ 1 + n
/ 2) * (2 * k
+ 1) / (4 * n
));
148 s
+= input
[i
] * cos(a
);
150 output
[k
] = REF_SCALE(s
, nbits
- 1);
153 #endif /* CONFIG_MDCT */
157 static void idct_ref(FFTSample
*output
, FFTSample
*input
, int nbits
)
159 int i
, k
, n
= 1 << nbits
;
162 for (i
= 0; i
< n
; i
++) {
163 double s
= 0.5 * input
[0];
164 for (k
= 1; k
< n
; k
++) {
165 double a
= M_PI
* k
* (i
+ 0.5) / n
;
166 s
+= input
[k
] * cos(a
);
168 output
[i
] = 2 * s
/ n
;
172 static void dct_ref(FFTSample
*output
, FFTSample
*input
, int nbits
)
174 int i
, k
, n
= 1 << nbits
;
177 for (k
= 0; k
< n
; k
++) {
179 for (i
= 0; i
< n
; i
++) {
180 double a
= M_PI
* k
* (i
+ 0.5) / n
;
181 s
+= input
[i
] * cos(a
);
186 #endif /* CONFIG_DCT */
187 #endif /* FFT_FLOAT */
189 static FFTSample
frandom(AVLFG
*prng
)
191 return (int16_t) av_lfg_get(prng
) / 32768.0 * RANGE
;
194 static int check_diff(FFTSample
*tab1
, FFTSample
*tab2
, int n
, double scale
)
197 double error
= 0, max
= 0;
199 for (i
= 0; i
< n
; i
++) {
200 double e
= fabsf(tab1
[i
] - (tab2
[i
] / scale
)) / RANGE
;
202 av_log(NULL
, AV_LOG_ERROR
, "ERROR %5d: "FMT
" "FMT
"\n",
203 i
, tab1
[i
], tab2
[i
]);
210 av_log(NULL
, AV_LOG_INFO
, "max:%f e:%g\n", max
, sqrt(error
/ n
));
214 static void help(void)
216 av_log(NULL
, AV_LOG_INFO
,
217 "usage: fft-test [-h] [-s] [-i] [-n b]\n"
218 "-h print this help\n"
223 "-i inverse transform test\n"
224 "-n b set the transform size to 2^b\n"
225 "-f x set scale factor for output data of (I)MDCT to x\n");
236 #include "compat/getopt.c"
239 int main(int argc
, char **argv
)
241 FFTComplex
*tab
, *tab1
, *tab_ref
;
243 enum tf_transform transform
= TRANSFORM_FFT
;
248 #endif /* FFT_FLOAT */
250 int do_speed
= 0, do_inverse
= 0;
251 int fft_nbits
= 9, fft_size
;
255 av_lfg_init(&prng
, 1);
258 int c
= getopt(argc
, argv
, "hsimrdn:f:c:");
272 transform
= TRANSFORM_MDCT
;
275 transform
= TRANSFORM_RDFT
;
278 transform
= TRANSFORM_DCT
;
281 fft_nbits
= atoi(optarg
);
284 scale
= atof(optarg
);
288 int cpuflags
= av_get_cpu_flags();
290 if (av_parse_cpu_caps(&cpuflags
, optarg
) < 0)
293 av_force_cpu_flags(cpuflags
);
299 fft_size
= 1 << fft_nbits
;
300 tab
= av_malloc_array(fft_size
, sizeof(FFTComplex
));
301 tab1
= av_malloc_array(fft_size
, sizeof(FFTComplex
));
302 tab_ref
= av_malloc_array(fft_size
, sizeof(FFTComplex
));
303 tab2
= av_malloc_array(fft_size
, sizeof(FFTSample
));
305 if (!(tab
&& tab1
&& tab_ref
&& tab2
))
311 av_log(NULL
, AV_LOG_INFO
, "Scale factor is set to %f\n", scale
);
313 av_log(NULL
, AV_LOG_INFO
, "IMDCT");
315 av_log(NULL
, AV_LOG_INFO
, "MDCT");
316 ff_mdct_init(&m
, fft_nbits
, do_inverse
, scale
);
318 #endif /* CONFIG_MDCT */
321 av_log(NULL
, AV_LOG_INFO
, "IFFT");
323 av_log(NULL
, AV_LOG_INFO
, "FFT");
324 ff_fft_init(&s
, fft_nbits
, do_inverse
);
325 if ((err
= fft_ref_init(fft_nbits
, do_inverse
)) < 0)
332 av_log(NULL
, AV_LOG_INFO
, "IDFT_C2R");
334 av_log(NULL
, AV_LOG_INFO
, "DFT_R2C");
335 ff_rdft_init(&r
, fft_nbits
, do_inverse
? IDFT_C2R
: DFT_R2C
);
336 if ((err
= fft_ref_init(fft_nbits
, do_inverse
)) < 0)
339 # endif /* CONFIG_RDFT */
343 av_log(NULL
, AV_LOG_INFO
, "DCT_III");
345 av_log(NULL
, AV_LOG_INFO
, "DCT_II");
346 ff_dct_init(&d
, fft_nbits
, do_inverse
? DCT_III
: DCT_II
);
348 # endif /* CONFIG_DCT */
349 #endif /* FFT_FLOAT */
351 av_log(NULL
, AV_LOG_ERROR
, "Requested transform not supported\n");
354 av_log(NULL
, AV_LOG_INFO
, " %d test\n", fft_size
);
356 /* generate random data */
358 for (i
= 0; i
< fft_size
; i
++) {
359 tab1
[i
].re
= frandom(&prng
);
360 tab1
[i
].im
= frandom(&prng
);
363 /* checking result */
364 av_log(NULL
, AV_LOG_INFO
, "Checking...\n");
370 imdct_ref(&tab_ref
->re
, &tab1
->re
, fft_nbits
);
371 m
.imdct_calc(&m
, tab2
, &tab1
->re
);
372 err
= check_diff(&tab_ref
->re
, tab2
, fft_size
, scale
);
374 mdct_ref(&tab_ref
->re
, &tab1
->re
, fft_nbits
);
375 m
.mdct_calc(&m
, tab2
, &tab1
->re
);
376 err
= check_diff(&tab_ref
->re
, tab2
, fft_size
/ 2, scale
);
379 #endif /* CONFIG_MDCT */
381 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
382 s
.fft_permute(&s
, tab
);
385 fft_ref(tab_ref
, tab1
, fft_nbits
);
386 err
= check_diff(&tab_ref
->re
, &tab
->re
, fft_size
* 2, 1.0);
392 int fft_size_2
= fft_size
>> 1;
395 tab1
[fft_size_2
].im
= 0;
396 for (i
= 1; i
< fft_size_2
; i
++) {
397 tab1
[fft_size_2
+ i
].re
= tab1
[fft_size_2
- i
].re
;
398 tab1
[fft_size_2
+ i
].im
= -tab1
[fft_size_2
- i
].im
;
401 memcpy(tab2
, tab1
, fft_size
* sizeof(FFTSample
));
402 tab2
[1] = tab1
[fft_size_2
].re
;
404 r
.rdft_calc(&r
, tab2
);
405 fft_ref(tab_ref
, tab1
, fft_nbits
);
406 for (i
= 0; i
< fft_size
; i
++) {
410 err
= check_diff(&tab_ref
->re
, &tab
->re
, fft_size
* 2, 0.5);
412 for (i
= 0; i
< fft_size
; i
++) {
413 tab2
[i
] = tab1
[i
].re
;
416 r
.rdft_calc(&r
, tab2
);
417 fft_ref(tab_ref
, tab1
, fft_nbits
);
418 tab_ref
[0].im
= tab_ref
[fft_size_2
].re
;
419 err
= check_diff(&tab_ref
->re
, tab2
, fft_size
, 1.0);
423 #endif /* CONFIG_RDFT */
426 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
427 d
.dct_calc(&d
, &tab
->re
);
429 idct_ref(&tab_ref
->re
, &tab1
->re
, fft_nbits
);
431 dct_ref(&tab_ref
->re
, &tab1
->re
, fft_nbits
);
432 err
= check_diff(&tab_ref
->re
, &tab
->re
, fft_size
, 1.0);
434 #endif /* CONFIG_DCT */
435 #endif /* FFT_FLOAT */
438 /* do a speed test */
441 int64_t time_start
, duration
;
444 av_log(NULL
, AV_LOG_INFO
, "Speed test...\n");
445 /* we measure during about 1 seconds */
448 time_start
= av_gettime_relative();
449 for (it
= 0; it
< nb_its
; it
++) {
453 m
.imdct_calc(&m
, &tab
->re
, &tab1
->re
);
455 m
.mdct_calc(&m
, &tab
->re
, &tab1
->re
);
458 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
463 memcpy(tab2
, tab1
, fft_size
* sizeof(FFTSample
));
464 r
.rdft_calc(&r
, tab2
);
467 memcpy(tab2
, tab1
, fft_size
* sizeof(FFTSample
));
468 d
.dct_calc(&d
, tab2
);
470 #endif /* FFT_FLOAT */
473 duration
= av_gettime_relative() - time_start
;
474 if (duration
>= 1000000)
478 av_log(NULL
, AV_LOG_INFO
,
479 "time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
480 (double) duration
/ nb_its
,
481 (double) duration
/ 1000000.0,
490 #endif /* CONFIG_MDCT */
499 # endif /* CONFIG_RDFT */
504 # endif /* CONFIG_DCT */
505 #endif /* FFT_FLOAT */
516 printf("Error: %d.\n", err
);