ffmpeg / libavcodec / ffttest.c @ 6efe6028
History  View  Annotate  Download (12.6 KB)
1 
/*


2 
* (c) 2002 Fabrice Bellard

3 
*

4 
* This file is part of Libav.

5 
*

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.

10 
*

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.

15 
*

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 021101301 USA

19 
*/

20  
21 
/**

22 
* @file

23 
* FFT and MDCT tests.

24 
*/

25  
26 
#include "libavutil/mathematics.h" 
27 
#include "libavutil/lfg.h" 
28 
#include "libavutil/log.h" 
29 
#include "fft.h" 
30 
#if CONFIG_FFT_FLOAT

31 
#include "dct.h" 
32 
#include "rdft.h" 
33 
#endif

34 
#include <math.h> 
35 
#include <unistd.h> 
36 
#include <sys/time.h> 
37 
#include <stdlib.h> 
38 
#include <string.h> 
39  
40 
#undef exit

41  
42 
/* reference fft */

43  
44 
#define MUL16(a,b) ((a) * (b))

45  
46 
#define CMAC(pre, pim, are, aim, bre, bim) \

47 
{\ 
48 
pre += (MUL16(are, bre)  MUL16(aim, bim));\ 
49 
pim += (MUL16(are, bim) + MUL16(bre, aim));\ 
50 
} 
51  
52 
#if CONFIG_FFT_FLOAT

53 
# define RANGE 1.0 
54 
# define REF_SCALE(x, bits) (x)

55 
# define FMT "%10.6f" 
56 
#else

57 
# define RANGE 16384 
58 
# define REF_SCALE(x, bits) ((x) / (1<<(bits))) 
59 
# define FMT "%6d" 
60 
#endif

61  
62 
struct {

63 
float re, im;

64 
} *exptab; 
65  
66 
static void fft_ref_init(int nbits, int inverse) 
67 
{ 
68 
int n, i;

69 
double c1, s1, alpha;

70  
71 
n = 1 << nbits;

72 
exptab = av_malloc((n / 2) * sizeof(*exptab)); 
73  
74 
for (i = 0; i < (n/2); i++) { 
75 
alpha = 2 * M_PI * (float)i / (float)n; 
76 
c1 = cos(alpha); 
77 
s1 = sin(alpha); 
78 
if (!inverse)

79 
s1 = s1; 
80 
exptab[i].re = c1; 
81 
exptab[i].im = s1; 
82 
} 
83 
} 
84  
85 
static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) 
86 
{ 
87 
int n, i, j, k, n2;

88 
double tmp_re, tmp_im, s, c;

89 
FFTComplex *q; 
90  
91 
n = 1 << nbits;

92 
n2 = n >> 1;

93 
for (i = 0; i < n; i++) { 
94 
tmp_re = 0;

95 
tmp_im = 0;

96 
q = tab; 
97 
for (j = 0; j < n; j++) { 
98 
k = (i * j) & (n  1);

99 
if (k >= n2) {

100 
c = exptab[k  n2].re; 
101 
s = exptab[k  n2].im; 
102 
} else {

103 
c = exptab[k].re; 
104 
s = exptab[k].im; 
105 
} 
106 
CMAC(tmp_re, tmp_im, c, s, q>re, q>im); 
107 
q++; 
108 
} 
109 
tabr[i].re = REF_SCALE(tmp_re, nbits); 
110 
tabr[i].im = REF_SCALE(tmp_im, nbits); 
111 
} 
112 
} 
113  
114 
static void imdct_ref(FFTSample *out, FFTSample *in, int nbits) 
115 
{ 
116 
int n = 1<<nbits; 
117 
int k, i, a;

118 
double sum, f;

119  
120 
for (i = 0; i < n; i++) { 
121 
sum = 0;

122 
for (k = 0; k < n/2; k++) { 
123 
a = (2 * i + 1 + (n / 2)) * (2 * k + 1); 
124 
f = cos(M_PI * a / (double)(2 * n)); 
125 
sum += f * in[k]; 
126 
} 
127 
out[i] = REF_SCALE(sum, nbits  2);

128 
} 
129 
} 
130  
131 
/* NOTE: no normalisation by 1 / N is done */

132 
static void mdct_ref(FFTSample *output, FFTSample *input, int nbits) 
133 
{ 
134 
int n = 1<<nbits; 
135 
int k, i;

136 
double a, s;

137  
138 
/* do it by hand */

139 
for (k = 0; k < n/2; k++) { 
140 
s = 0;

141 
for (i = 0; i < n; i++) { 
142 
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); 
143 
s += input[i] * cos(a); 
144 
} 
145 
output[k] = REF_SCALE(s, nbits  1);

146 
} 
147 
} 
148  
149 
#if CONFIG_FFT_FLOAT

150 
static void idct_ref(float *output, float *input, int nbits) 
151 
{ 
152 
int n = 1<<nbits; 
153 
int k, i;

154 
double a, s;

155  
156 
/* do it by hand */

157 
for (i = 0; i < n; i++) { 
158 
s = 0.5 * input[0]; 
159 
for (k = 1; k < n; k++) { 
160 
a = M_PI*k*(i+0.5) / n; 
161 
s += input[k] * cos(a); 
162 
} 
163 
output[i] = 2 * s / n;

164 
} 
165 
} 
166 
static void dct_ref(float *output, float *input, int nbits) 
167 
{ 
168 
int n = 1<<nbits; 
169 
int k, i;

170 
double a, s;

171  
172 
/* do it by hand */

173 
for (k = 0; k < n; k++) { 
174 
s = 0;

175 
for (i = 0; i < n; i++) { 
176 
a = M_PI*k*(i+0.5) / n; 
177 
s += input[i] * cos(a); 
178 
} 
179 
output[k] = s; 
180 
} 
181 
} 
182 
#endif

183  
184  
185 
static FFTSample frandom(AVLFG *prng)

186 
{ 
187 
return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE; 
188 
} 
189  
190 
static int64_t gettime(void) 
191 
{ 
192 
struct timeval tv;

193 
gettimeofday(&tv,NULL);

194 
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; 
195 
} 
196  
197 
static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale) 
198 
{ 
199 
int i;

200 
double max= 0; 
201 
double error= 0; 
202 
int err = 0; 
203  
204 
for (i = 0; i < n; i++) { 
205 
double e = fabsf(tab1[i]  (tab2[i] / scale)) / RANGE;

206 
if (e >= 1e3) { 
207 
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n", 
208 
i, tab1[i], tab2[i]); 
209 
err = 1;

210 
} 
211 
error+= e*e; 
212 
if(e>max) max= e;

213 
} 
214 
av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n); 
215 
return err;

216 
} 
217  
218  
219 
static void help(void) 
220 
{ 
221 
av_log(NULL, AV_LOG_INFO,"usage: ffttest [h] [s] [i] [n b]\n" 
222 
"h print this help\n"

223 
"s speed test\n"

224 
"m (I)MDCT test\n"

225 
"d (I)DCT test\n"

226 
"r (I)RDFT test\n"

227 
"i inverse transform test\n"

228 
"n b set the transform size to 2^b\n"

229 
"f x set scale factor for output data of (I)MDCT to x\n"

230 
); 
231 
exit(1);

232 
} 
233  
234 
enum tf_transform {

235 
TRANSFORM_FFT, 
236 
TRANSFORM_MDCT, 
237 
TRANSFORM_RDFT, 
238 
TRANSFORM_DCT, 
239 
}; 
240  
241 
int main(int argc, char **argv) 
242 
{ 
243 
FFTComplex *tab, *tab1, *tab_ref; 
244 
FFTSample *tab2; 
245 
int it, i, c;

246 
int do_speed = 0; 
247 
int err = 1; 
248 
enum tf_transform transform = TRANSFORM_FFT;

249 
int do_inverse = 0; 
250 
FFTContext s1, *s = &s1; 
251 
FFTContext m1, *m = &m1; 
252 
#if CONFIG_FFT_FLOAT

253 
RDFTContext r1, *r = &r1; 
254 
DCTContext d1, *d = &d1; 
255 
#endif

256 
int fft_nbits, fft_size, fft_size_2;

257 
double scale = 1.0; 
258 
AVLFG prng; 
259 
av_lfg_init(&prng, 1);

260  
261 
fft_nbits = 9;

262 
for(;;) {

263 
c = getopt(argc, argv, "hsimrdn:f:");

264 
if (c == 1) 
265 
break;

266 
switch(c) {

267 
case 'h': 
268 
help(); 
269 
break;

270 
case 's': 
271 
do_speed = 1;

272 
break;

273 
case 'i': 
274 
do_inverse = 1;

275 
break;

276 
case 'm': 
277 
transform = TRANSFORM_MDCT; 
278 
break;

279 
case 'r': 
280 
transform = TRANSFORM_RDFT; 
281 
break;

282 
case 'd': 
283 
transform = TRANSFORM_DCT; 
284 
break;

285 
case 'n': 
286 
fft_nbits = atoi(optarg); 
287 
break;

288 
case 'f': 
289 
scale = atof(optarg); 
290 
break;

291 
} 
292 
} 
293  
294 
fft_size = 1 << fft_nbits;

295 
fft_size_2 = fft_size >> 1;

296 
tab = av_malloc(fft_size * sizeof(FFTComplex));

297 
tab1 = av_malloc(fft_size * sizeof(FFTComplex));

298 
tab_ref = av_malloc(fft_size * sizeof(FFTComplex));

299 
tab2 = av_malloc(fft_size * sizeof(FFTSample));

300  
301 
switch (transform) {

302 
case TRANSFORM_MDCT:

303 
av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale); 
304 
if (do_inverse)

305 
av_log(NULL, AV_LOG_INFO,"IMDCT"); 
306 
else

307 
av_log(NULL, AV_LOG_INFO,"MDCT"); 
308 
ff_mdct_init(m, fft_nbits, do_inverse, scale); 
309 
break;

310 
case TRANSFORM_FFT:

311 
if (do_inverse)

312 
av_log(NULL, AV_LOG_INFO,"IFFT"); 
313 
else

314 
av_log(NULL, AV_LOG_INFO,"FFT"); 
315 
ff_fft_init(s, fft_nbits, do_inverse); 
316 
fft_ref_init(fft_nbits, do_inverse); 
317 
break;

318 
#if CONFIG_FFT_FLOAT

319 
case TRANSFORM_RDFT:

320 
if (do_inverse)

321 
av_log(NULL, AV_LOG_INFO,"IDFT_C2R"); 
322 
else

323 
av_log(NULL, AV_LOG_INFO,"DFT_R2C"); 
324 
ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C); 
325 
fft_ref_init(fft_nbits, do_inverse); 
326 
break;

327 
case TRANSFORM_DCT:

328 
if (do_inverse)

329 
av_log(NULL, AV_LOG_INFO,"DCT_III"); 
330 
else

331 
av_log(NULL, AV_LOG_INFO,"DCT_II"); 
332 
ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II); 
333 
break;

334 
#endif

335 
default:

336 
av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n"); 
337 
return 1; 
338 
} 
339 
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 
340  
341 
/* generate random data */

342  
343 
for (i = 0; i < fft_size; i++) { 
344 
tab1[i].re = frandom(&prng); 
345 
tab1[i].im = frandom(&prng); 
346 
} 
347  
348 
/* checking result */

349 
av_log(NULL, AV_LOG_INFO,"Checking...\n"); 
350  
351 
switch (transform) {

352 
case TRANSFORM_MDCT:

353 
if (do_inverse) {

354 
imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits); 
355 
m>imdct_calc(m, tab2, (FFTSample *)tab1); 
356 
err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale); 
357 
} else {

358 
mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits); 
359  
360 
m>mdct_calc(m, tab2, (FFTSample *)tab1); 
361  
362 
err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);

363 
} 
364 
break;

365 
case TRANSFORM_FFT:

366 
memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

367 
s>fft_permute(s, tab); 
368 
s>fft_calc(s, tab); 
369  
370 
fft_ref(tab_ref, tab1, fft_nbits); 
371 
err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0); 
372 
break;

373 
#if CONFIG_FFT_FLOAT

374 
case TRANSFORM_RDFT:

375 
if (do_inverse) {

376 
tab1[ 0].im = 0; 
377 
tab1[fft_size_2].im = 0;

378 
for (i = 1; i < fft_size_2; i++) { 
379 
tab1[fft_size_2+i].re = tab1[fft_size_2i].re; 
380 
tab1[fft_size_2+i].im = tab1[fft_size_2i].im; 
381 
} 
382  
383 
memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

384 
tab2[1] = tab1[fft_size_2].re;

385  
386 
r>rdft_calc(r, tab2); 
387 
fft_ref(tab_ref, tab1, fft_nbits); 
388 
for (i = 0; i < fft_size; i++) { 
389 
tab[i].re = tab2[i]; 
390 
tab[i].im = 0;

391 
} 
392 
err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5); 
393 
} else {

394 
for (i = 0; i < fft_size; i++) { 
395 
tab2[i] = tab1[i].re; 
396 
tab1[i].im = 0;

397 
} 
398 
r>rdft_calc(r, tab2); 
399 
fft_ref(tab_ref, tab1, fft_nbits); 
400 
tab_ref[0].im = tab_ref[fft_size_2].re;

401 
err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0); 
402 
} 
403 
break;

404 
case TRANSFORM_DCT:

405 
memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

406 
d>dct_calc(d, tab); 
407 
if (do_inverse) {

408 
idct_ref(tab_ref, tab1, fft_nbits); 
409 
} else {

410 
dct_ref(tab_ref, tab1, fft_nbits); 
411 
} 
412 
err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0); 
413 
break;

414 
#endif

415 
} 
416  
417 
/* do a speed test */

418  
419 
if (do_speed) {

420 
int64_t time_start, duration; 
421 
int nb_its;

422  
423 
av_log(NULL, AV_LOG_INFO,"Speed test...\n"); 
424 
/* we measure during about 1 seconds */

425 
nb_its = 1;

426 
for(;;) {

427 
time_start = gettime(); 
428 
for (it = 0; it < nb_its; it++) { 
429 
switch (transform) {

430 
case TRANSFORM_MDCT:

431 
if (do_inverse) {

432 
m>imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1); 
433 
} else {

434 
m>mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1); 
435 
} 
436 
break;

437 
case TRANSFORM_FFT:

438 
memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

439 
s>fft_calc(s, tab); 
440 
break;

441 
#if CONFIG_FFT_FLOAT

442 
case TRANSFORM_RDFT:

443 
memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

444 
r>rdft_calc(r, tab2); 
445 
break;

446 
case TRANSFORM_DCT:

447 
memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

448 
d>dct_calc(d, tab2); 
449 
break;

450 
#endif

451 
} 
452 
} 
453 
duration = gettime()  time_start; 
454 
if (duration >= 1000000) 
455 
break;

456 
nb_its *= 2;

457 
} 
458 
av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n", 
459 
(double)duration / nb_its,

460 
(double)duration / 1000000.0, 
461 
nb_its); 
462 
} 
463  
464 
switch (transform) {

465 
case TRANSFORM_MDCT:

466 
ff_mdct_end(m); 
467 
break;

468 
case TRANSFORM_FFT:

469 
ff_fft_end(s); 
470 
break;

471 
#if CONFIG_FFT_FLOAT

472 
case TRANSFORM_RDFT:

473 
ff_rdft_end(r); 
474 
break;

475 
case TRANSFORM_DCT:

476 
ff_dct_end(d); 
477 
break;

478 
#endif

479 
} 
480  
481 
av_free(tab); 
482 
av_free(tab1); 
483 
av_free(tab2); 
484 
av_free(tab_ref); 
485 
av_free(exptab); 
486  
487 
return err;

488 
} 