ffmpeg / libavcodec / ffttest.c @ 2912e87a
History  View  Annotate  Download (11.9 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 
#include <math.h> 
31 
#include <unistd.h> 
32 
#include <sys/time.h> 
33 
#include <stdlib.h> 
34 
#include <string.h> 
35  
36 
#undef exit

37  
38 
/* reference fft */

39  
40 
#define MUL16(a,b) ((a) * (b))

41  
42 
#define CMAC(pre, pim, are, aim, bre, bim) \

43 
{\ 
44 
pre += (MUL16(are, bre)  MUL16(aim, bim));\ 
45 
pim += (MUL16(are, bim) + MUL16(bre, aim));\ 
46 
} 
47  
48 
FFTComplex *exptab; 
49  
50 
static void fft_ref_init(int nbits, int inverse) 
51 
{ 
52 
int n, i;

53 
double c1, s1, alpha;

54  
55 
n = 1 << nbits;

56 
exptab = av_malloc((n / 2) * sizeof(FFTComplex)); 
57  
58 
for (i = 0; i < (n/2); i++) { 
59 
alpha = 2 * M_PI * (float)i / (float)n; 
60 
c1 = cos(alpha); 
61 
s1 = sin(alpha); 
62 
if (!inverse)

63 
s1 = s1; 
64 
exptab[i].re = c1; 
65 
exptab[i].im = s1; 
66 
} 
67 
} 
68  
69 
static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) 
70 
{ 
71 
int n, i, j, k, n2;

72 
double tmp_re, tmp_im, s, c;

73 
FFTComplex *q; 
74  
75 
n = 1 << nbits;

76 
n2 = n >> 1;

77 
for (i = 0; i < n; i++) { 
78 
tmp_re = 0;

79 
tmp_im = 0;

80 
q = tab; 
81 
for (j = 0; j < n; j++) { 
82 
k = (i * j) & (n  1);

83 
if (k >= n2) {

84 
c = exptab[k  n2].re; 
85 
s = exptab[k  n2].im; 
86 
} else {

87 
c = exptab[k].re; 
88 
s = exptab[k].im; 
89 
} 
90 
CMAC(tmp_re, tmp_im, c, s, q>re, q>im); 
91 
q++; 
92 
} 
93 
tabr[i].re = tmp_re; 
94 
tabr[i].im = tmp_im; 
95 
} 
96 
} 
97  
98 
static void imdct_ref(float *out, float *in, int nbits) 
99 
{ 
100 
int n = 1<<nbits; 
101 
int k, i, a;

102 
double sum, f;

103  
104 
for (i = 0; i < n; i++) { 
105 
sum = 0;

106 
for (k = 0; k < n/2; k++) { 
107 
a = (2 * i + 1 + (n / 2)) * (2 * k + 1); 
108 
f = cos(M_PI * a / (double)(2 * n)); 
109 
sum += f * in[k]; 
110 
} 
111 
out[i] = sum; 
112 
} 
113 
} 
114  
115 
/* NOTE: no normalisation by 1 / N is done */

116 
static void mdct_ref(float *output, float *input, int nbits) 
117 
{ 
118 
int n = 1<<nbits; 
119 
int k, i;

120 
double a, s;

121  
122 
/* do it by hand */

123 
for (k = 0; k < n/2; k++) { 
124 
s = 0;

125 
for (i = 0; i < n; i++) { 
126 
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); 
127 
s += input[i] * cos(a); 
128 
} 
129 
output[k] = s; 
130 
} 
131 
} 
132  
133 
static void idct_ref(float *output, float *input, int nbits) 
134 
{ 
135 
int n = 1<<nbits; 
136 
int k, i;

137 
double a, s;

138  
139 
/* do it by hand */

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

147 
} 
148 
} 
149 
static void dct_ref(float *output, float *input, int nbits) 
150 
{ 
151 
int n = 1<<nbits; 
152 
int k, i;

153 
double a, s;

154  
155 
/* do it by hand */

156 
for (k = 0; k < n; k++) { 
157 
s = 0;

158 
for (i = 0; i < n; i++) { 
159 
a = M_PI*k*(i+0.5) / n; 
160 
s += input[i] * cos(a); 
161 
} 
162 
output[k] = s; 
163 
} 
164 
} 
165  
166  
167 
static float frandom(AVLFG *prng) 
168 
{ 
169 
return (int16_t)av_lfg_get(prng) / 32768.0; 
170 
} 
171  
172 
static int64_t gettime(void) 
173 
{ 
174 
struct timeval tv;

175 
gettimeofday(&tv,NULL);

176 
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; 
177 
} 
178  
179 
static int check_diff(float *tab1, float *tab2, int n, double scale) 
180 
{ 
181 
int i;

182 
double max= 0; 
183 
double error= 0; 
184 
int err = 0; 
185  
186 
for (i = 0; i < n; i++) { 
187 
double e= fabsf(tab1[i]  (tab2[i] / scale));

188 
if (e >= 1e3) { 
189 
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: %10.6f %10.6f\n", 
190 
i, tab1[i], tab2[i]); 
191 
err = 1;

192 
} 
193 
error+= e*e; 
194 
if(e>max) max= e;

195 
} 
196 
av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n); 
197 
return err;

198 
} 
199  
200  
201 
static void help(void) 
202 
{ 
203 
av_log(NULL, AV_LOG_INFO,"usage: ffttest [h] [s] [i] [n b]\n" 
204 
"h print this help\n"

205 
"s speed test\n"

206 
"m (I)MDCT test\n"

207 
"d (I)DCT test\n"

208 
"r (I)RDFT test\n"

209 
"i inverse transform test\n"

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

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

212 
); 
213 
exit(1);

214 
} 
215  
216 
enum tf_transform {

217 
TRANSFORM_FFT, 
218 
TRANSFORM_MDCT, 
219 
TRANSFORM_RDFT, 
220 
TRANSFORM_DCT, 
221 
}; 
222  
223 
int main(int argc, char **argv) 
224 
{ 
225 
FFTComplex *tab, *tab1, *tab_ref; 
226 
FFTSample *tab2; 
227 
int it, i, c;

228 
int do_speed = 0; 
229 
int err = 1; 
230 
enum tf_transform transform = TRANSFORM_FFT;

231 
int do_inverse = 0; 
232 
FFTContext s1, *s = &s1; 
233 
FFTContext m1, *m = &m1; 
234 
RDFTContext r1, *r = &r1; 
235 
DCTContext d1, *d = &d1; 
236 
int fft_nbits, fft_size, fft_size_2;

237 
double scale = 1.0; 
238 
AVLFG prng; 
239 
av_lfg_init(&prng, 1);

240  
241 
fft_nbits = 9;

242 
for(;;) {

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

244 
if (c == 1) 
245 
break;

246 
switch(c) {

247 
case 'h': 
248 
help(); 
249 
break;

250 
case 's': 
251 
do_speed = 1;

252 
break;

253 
case 'i': 
254 
do_inverse = 1;

255 
break;

256 
case 'm': 
257 
transform = TRANSFORM_MDCT; 
258 
break;

259 
case 'r': 
260 
transform = TRANSFORM_RDFT; 
261 
break;

262 
case 'd': 
263 
transform = TRANSFORM_DCT; 
264 
break;

265 
case 'n': 
266 
fft_nbits = atoi(optarg); 
267 
break;

268 
case 'f': 
269 
scale = atof(optarg); 
270 
break;

271 
} 
272 
} 
273  
274 
fft_size = 1 << fft_nbits;

275 
fft_size_2 = fft_size >> 1;

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

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

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

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

280  
281 
switch (transform) {

282 
case TRANSFORM_MDCT:

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

285 
av_log(NULL, AV_LOG_INFO,"IMDCT"); 
286 
else

287 
av_log(NULL, AV_LOG_INFO,"MDCT"); 
288 
ff_mdct_init(m, fft_nbits, do_inverse, scale); 
289 
break;

290 
case TRANSFORM_FFT:

291 
if (do_inverse)

292 
av_log(NULL, AV_LOG_INFO,"IFFT"); 
293 
else

294 
av_log(NULL, AV_LOG_INFO,"FFT"); 
295 
ff_fft_init(s, fft_nbits, do_inverse); 
296 
fft_ref_init(fft_nbits, do_inverse); 
297 
break;

298 
case TRANSFORM_RDFT:

299 
if (do_inverse)

300 
av_log(NULL, AV_LOG_INFO,"IDFT_C2R"); 
301 
else

302 
av_log(NULL, AV_LOG_INFO,"DFT_R2C"); 
303 
ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C); 
304 
fft_ref_init(fft_nbits, do_inverse); 
305 
break;

306 
case TRANSFORM_DCT:

307 
if (do_inverse)

308 
av_log(NULL, AV_LOG_INFO,"DCT_III"); 
309 
else

310 
av_log(NULL, AV_LOG_INFO,"DCT_II"); 
311 
ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II); 
312 
break;

313 
} 
314 
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 
315  
316 
/* generate random data */

317  
318 
for (i = 0; i < fft_size; i++) { 
319 
tab1[i].re = frandom(&prng); 
320 
tab1[i].im = frandom(&prng); 
321 
} 
322  
323 
/* checking result */

324 
av_log(NULL, AV_LOG_INFO,"Checking...\n"); 
325  
326 
switch (transform) {

327 
case TRANSFORM_MDCT:

328 
if (do_inverse) {

329 
imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); 
330 
ff_imdct_calc(m, tab2, (float *)tab1);

331 
err = check_diff((float *)tab_ref, tab2, fft_size, scale);

332 
} else {

333 
mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); 
334  
335 
ff_mdct_calc(m, tab2, (float *)tab1);

336  
337 
err = check_diff((float *)tab_ref, tab2, fft_size / 2, scale); 
338 
} 
339 
break;

340 
case TRANSFORM_FFT:

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

342 
ff_fft_permute(s, tab); 
343 
ff_fft_calc(s, tab); 
344  
345 
fft_ref(tab_ref, tab1, fft_nbits); 
346 
err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 1.0); 
347 
break;

348 
case TRANSFORM_RDFT:

349 
if (do_inverse) {

350 
tab1[ 0].im = 0; 
351 
tab1[fft_size_2].im = 0;

352 
for (i = 1; i < fft_size_2; i++) { 
353 
tab1[fft_size_2+i].re = tab1[fft_size_2i].re; 
354 
tab1[fft_size_2+i].im = tab1[fft_size_2i].im; 
355 
} 
356  
357 
memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

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

359  
360 
ff_rdft_calc(r, tab2); 
361 
fft_ref(tab_ref, tab1, fft_nbits); 
362 
for (i = 0; i < fft_size; i++) { 
363 
tab[i].re = tab2[i]; 
364 
tab[i].im = 0;

365 
} 
366 
err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5); 
367 
} else {

368 
for (i = 0; i < fft_size; i++) { 
369 
tab2[i] = tab1[i].re; 
370 
tab1[i].im = 0;

371 
} 
372 
ff_rdft_calc(r, tab2); 
373 
fft_ref(tab_ref, tab1, fft_nbits); 
374 
tab_ref[0].im = tab_ref[fft_size_2].re;

375 
err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0); 
376 
} 
377 
break;

378 
case TRANSFORM_DCT:

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

380 
ff_dct_calc(d, tab); 
381 
if (do_inverse) {

382 
idct_ref(tab_ref, tab1, fft_nbits); 
383 
} else {

384 
dct_ref(tab_ref, tab1, fft_nbits); 
385 
} 
386 
err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0); 
387 
break;

388 
} 
389  
390 
/* do a speed test */

391  
392 
if (do_speed) {

393 
int64_t time_start, duration; 
394 
int nb_its;

395  
396 
av_log(NULL, AV_LOG_INFO,"Speed test...\n"); 
397 
/* we measure during about 1 seconds */

398 
nb_its = 1;

399 
for(;;) {

400 
time_start = gettime(); 
401 
for (it = 0; it < nb_its; it++) { 
402 
switch (transform) {

403 
case TRANSFORM_MDCT:

404 
if (do_inverse) {

405 
ff_imdct_calc(m, (float *)tab, (float *)tab1); 
406 
} else {

407 
ff_mdct_calc(m, (float *)tab, (float *)tab1); 
408 
} 
409 
break;

410 
case TRANSFORM_FFT:

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

412 
ff_fft_calc(s, tab); 
413 
break;

414 
case TRANSFORM_RDFT:

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

416 
ff_rdft_calc(r, tab2); 
417 
break;

418 
case TRANSFORM_DCT:

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

420 
ff_dct_calc(d, tab2); 
421 
break;

422 
} 
423 
} 
424 
duration = gettime()  time_start; 
425 
if (duration >= 1000000) 
426 
break;

427 
nb_its *= 2;

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

431 
(double)duration / 1000000.0, 
432 
nb_its); 
433 
} 
434  
435 
switch (transform) {

436 
case TRANSFORM_MDCT:

437 
ff_mdct_end(m); 
438 
break;

439 
case TRANSFORM_FFT:

440 
ff_fft_end(s); 
441 
break;

442 
case TRANSFORM_RDFT:

443 
ff_rdft_end(r); 
444 
break;

445 
case TRANSFORM_DCT:

446 
ff_dct_end(d); 
447 
break;

448 
} 
449  
450 
av_free(tab); 
451 
av_free(tab1); 
452 
av_free(tab2); 
453 
av_free(tab_ref); 
454 
av_free(exptab); 
455  
456 
return err;

457 
} 