2 * (c) 2002 Fabrice Bellard
4 * This file is part of Libav.
6 * Libav 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 * Libav 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 Libav; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
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"
45 #define MUL16(a,b) ((a) * (b))
47 #define CMAC(pre, pim, are, aim, bre, bim) \
49 pre += (MUL16(are, bre) - MUL16(aim, bim));\
50 pim += (MUL16(are, bim) + MUL16(bre, aim));\
55 # define REF_SCALE(x, bits) (x)
59 # define REF_SCALE(x, bits) ((x) / (1<<(bits)))
67 static void fft_ref_init(int nbits
, int inverse
)
73 exptab
= av_malloc((n
/ 2) * sizeof(*exptab
));
75 for (i
= 0; i
< (n
/2); i
++) {
76 alpha
= 2 * M_PI
* (float)i
/ (float)n
;
86 static void fft_ref(FFTComplex
*tabr
, FFTComplex
*tab
, int nbits
)
89 double tmp_re
, tmp_im
, s
, c
;
94 for (i
= 0; i
< n
; i
++) {
98 for (j
= 0; j
< n
; j
++) {
99 k
= (i
* j
) & (n
- 1);
101 c
= -exptab
[k
- n2
].re
;
102 s
= -exptab
[k
- n2
].im
;
107 CMAC(tmp_re
, tmp_im
, c
, s
, q
->re
, q
->im
);
110 tabr
[i
].re
= REF_SCALE(tmp_re
, nbits
);
111 tabr
[i
].im
= REF_SCALE(tmp_im
, nbits
);
115 static void imdct_ref(FFTSample
*out
, FFTSample
*in
, int nbits
)
121 for (i
= 0; i
< n
; i
++) {
123 for (k
= 0; k
< n
/2; k
++) {
124 a
= (2 * i
+ 1 + (n
/ 2)) * (2 * k
+ 1);
125 f
= cos(M_PI
* a
/ (double)(2 * n
));
128 out
[i
] = REF_SCALE(-sum
, nbits
- 2);
132 /* NOTE: no normalisation by 1 / N is done */
133 static void mdct_ref(FFTSample
*output
, FFTSample
*input
, int nbits
)
140 for (k
= 0; k
< n
/2; k
++) {
142 for (i
= 0; i
< n
; i
++) {
143 a
= (2*M_PI
*(2*i
+1+n
/2)*(2*k
+1) / (4 * n
));
144 s
+= input
[i
] * cos(a
);
146 output
[k
] = REF_SCALE(s
, nbits
- 1);
151 static void idct_ref(float *output
, float *input
, int nbits
)
158 for (i
= 0; i
< n
; i
++) {
160 for (k
= 1; k
< n
; k
++) {
161 a
= M_PI
*k
*(i
+0.5) / n
;
162 s
+= input
[k
] * cos(a
);
164 output
[i
] = 2 * s
/ n
;
167 static void dct_ref(float *output
, float *input
, int nbits
)
174 for (k
= 0; k
< n
; k
++) {
176 for (i
= 0; i
< n
; i
++) {
177 a
= M_PI
*k
*(i
+0.5) / n
;
178 s
+= input
[i
] * cos(a
);
186 static FFTSample
frandom(AVLFG
*prng
)
188 return (int16_t)av_lfg_get(prng
) / 32768.0 * RANGE
;
191 static int check_diff(FFTSample
*tab1
, FFTSample
*tab2
, int n
, double scale
)
198 for (i
= 0; i
< n
; i
++) {
199 double e
= fabsf(tab1
[i
] - (tab2
[i
] / scale
)) / RANGE
;
201 av_log(NULL
, AV_LOG_ERROR
, "ERROR %5d: "FMT
" "FMT
"\n",
202 i
, tab1
[i
], tab2
[i
]);
208 av_log(NULL
, AV_LOG_INFO
, "max:%f e:%g\n", max
, sqrt(error
)/n
);
213 static void help(void)
215 av_log(NULL
, AV_LOG_INFO
,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
216 "-h print this help\n"
221 "-i inverse transform test\n"
222 "-n b set the transform size to 2^b\n"
223 "-f x set scale factor for output data of (I)MDCT to x\n"
235 #include "compat/getopt.c"
238 int main(int argc
, char **argv
)
240 FFTComplex
*tab
, *tab1
, *tab_ref
;
246 enum tf_transform transform
= TRANSFORM_FFT
;
248 FFTContext s1
, *s
= &s1
;
249 FFTContext m1
, *m
= &m1
;
251 RDFTContext r1
, *r
= &r1
;
252 DCTContext d1
, *d
= &d1
;
255 int fft_nbits
, fft_size
;
258 av_lfg_init(&prng
, 1);
262 c
= getopt(argc
, argv
, "hsimrdn:f:c:");
276 transform
= TRANSFORM_MDCT
;
279 transform
= TRANSFORM_RDFT
;
282 transform
= TRANSFORM_DCT
;
285 fft_nbits
= atoi(optarg
);
288 scale
= atof(optarg
);
291 cpuflags
= av_parse_cpu_flags(optarg
);
294 av_set_cpu_flags_mask(cpuflags
);
299 fft_size
= 1 << fft_nbits
;
300 tab
= av_malloc(fft_size
* sizeof(FFTComplex
));
301 tab1
= av_malloc(fft_size
* sizeof(FFTComplex
));
302 tab_ref
= av_malloc(fft_size
* sizeof(FFTComplex
));
303 tab2
= av_malloc(fft_size
* sizeof(FFTSample
));
307 av_log(NULL
, AV_LOG_INFO
,"Scale factor is set to %f\n", scale
);
309 av_log(NULL
, AV_LOG_INFO
,"IMDCT");
311 av_log(NULL
, AV_LOG_INFO
,"MDCT");
312 ff_mdct_init(m
, fft_nbits
, do_inverse
, scale
);
316 av_log(NULL
, AV_LOG_INFO
,"IFFT");
318 av_log(NULL
, AV_LOG_INFO
,"FFT");
319 ff_fft_init(s
, fft_nbits
, do_inverse
);
320 fft_ref_init(fft_nbits
, do_inverse
);
325 av_log(NULL
, AV_LOG_INFO
,"IDFT_C2R");
327 av_log(NULL
, AV_LOG_INFO
,"DFT_R2C");
328 ff_rdft_init(r
, fft_nbits
, do_inverse
? IDFT_C2R
: DFT_R2C
);
329 fft_ref_init(fft_nbits
, do_inverse
);
333 av_log(NULL
, AV_LOG_INFO
,"DCT_III");
335 av_log(NULL
, AV_LOG_INFO
,"DCT_II");
336 ff_dct_init(d
, fft_nbits
, do_inverse
? DCT_III
: DCT_II
);
340 av_log(NULL
, AV_LOG_ERROR
, "Requested transform not supported\n");
343 av_log(NULL
, AV_LOG_INFO
," %d test\n", fft_size
);
345 /* generate random data */
347 for (i
= 0; i
< fft_size
; i
++) {
348 tab1
[i
].re
= frandom(&prng
);
349 tab1
[i
].im
= frandom(&prng
);
352 /* checking result */
353 av_log(NULL
, AV_LOG_INFO
,"Checking...\n");
358 imdct_ref((FFTSample
*)tab_ref
, (FFTSample
*)tab1
, fft_nbits
);
359 m
->imdct_calc(m
, tab2
, (FFTSample
*)tab1
);
360 err
= check_diff((FFTSample
*)tab_ref
, tab2
, fft_size
, scale
);
362 mdct_ref((FFTSample
*)tab_ref
, (FFTSample
*)tab1
, fft_nbits
);
364 m
->mdct_calc(m
, tab2
, (FFTSample
*)tab1
);
366 err
= check_diff((FFTSample
*)tab_ref
, tab2
, fft_size
/ 2, scale
);
370 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
371 s
->fft_permute(s
, tab
);
374 fft_ref(tab_ref
, tab1
, fft_nbits
);
375 err
= check_diff((FFTSample
*)tab_ref
, (FFTSample
*)tab
, fft_size
* 2, 1.0);
379 fft_size_2
= fft_size
>> 1;
382 tab1
[fft_size_2
].im
= 0;
383 for (i
= 1; i
< fft_size_2
; i
++) {
384 tab1
[fft_size_2
+i
].re
= tab1
[fft_size_2
-i
].re
;
385 tab1
[fft_size_2
+i
].im
= -tab1
[fft_size_2
-i
].im
;
388 memcpy(tab2
, tab1
, fft_size
* sizeof(FFTSample
));
389 tab2
[1] = tab1
[fft_size_2
].re
;
391 r
->rdft_calc(r
, tab2
);
392 fft_ref(tab_ref
, tab1
, fft_nbits
);
393 for (i
= 0; i
< fft_size
; i
++) {
397 err
= check_diff((float *)tab_ref
, (float *)tab
, fft_size
* 2, 0.5);
399 for (i
= 0; i
< fft_size
; i
++) {
400 tab2
[i
] = tab1
[i
].re
;
403 r
->rdft_calc(r
, tab2
);
404 fft_ref(tab_ref
, tab1
, fft_nbits
);
405 tab_ref
[0].im
= tab_ref
[fft_size_2
].re
;
406 err
= check_diff((float *)tab_ref
, (float *)tab2
, fft_size
, 1.0);
410 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
413 idct_ref(tab_ref
, tab1
, fft_nbits
);
415 dct_ref(tab_ref
, tab1
, fft_nbits
);
417 err
= check_diff((float *)tab_ref
, (float *)tab
, fft_size
, 1.0);
422 /* do a speed test */
425 int64_t time_start
, duration
;
428 av_log(NULL
, AV_LOG_INFO
,"Speed test...\n");
429 /* we measure during about 1 seconds */
432 time_start
= av_gettime();
433 for (it
= 0; it
< nb_its
; it
++) {
437 m
->imdct_calc(m
, (FFTSample
*)tab
, (FFTSample
*)tab1
);
439 m
->mdct_calc(m
, (FFTSample
*)tab
, (FFTSample
*)tab1
);
443 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
448 memcpy(tab2
, tab1
, fft_size
* sizeof(FFTSample
));
449 r
->rdft_calc(r
, tab2
);
452 memcpy(tab2
, tab1
, fft_size
* sizeof(FFTSample
));
453 d
->dct_calc(d
, tab2
);
458 duration
= av_gettime() - time_start
;
459 if (duration
>= 1000000)
463 av_log(NULL
, AV_LOG_INFO
,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
464 (double)duration
/ nb_its
,
465 (double)duration
/ 1000000.0,