Revision 6efe6028 libavcodec/ffttest.c
libavcodec/ffttest.c  

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

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

32  34 
#include <math.h> 
33  35 
#include <unistd.h> 
34  36 
#include <sys/time.h> 
...  ...  
47  49 
pim += (MUL16(are, bim) + MUL16(bre, aim));\ 
48  50 
} 
49  51  
50 
FFTComplex *exptab; 

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; 

51  65  
52  66 
static void fft_ref_init(int nbits, int inverse) 
53  67 
{ 
...  ...  
55  69 
double c1, s1, alpha; 
56  70  
57  71 
n = 1 << nbits; 
58 
exptab = av_malloc((n / 2) * sizeof(FFTComplex));


72 
exptab = av_malloc((n / 2) * sizeof(*exptab));


59  73  
60  74 
for (i = 0; i < (n/2); i++) { 
61  75 
alpha = 2 * M_PI * (float)i / (float)n; 
...  ...  
92  106 
CMAC(tmp_re, tmp_im, c, s, q>re, q>im); 
93  107 
q++; 
94  108 
} 
95 
tabr[i].re = tmp_re;


96 
tabr[i].im = tmp_im;


109 
tabr[i].re = REF_SCALE(tmp_re, nbits);


110 
tabr[i].im = REF_SCALE(tmp_im, nbits);


97  111 
} 
98  112 
} 
99  113  
100 
static void imdct_ref(float *out, float *in, int nbits)


114 
static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)


101  115 
{ 
102  116 
int n = 1<<nbits; 
103  117 
int k, i, a; 
...  ...  
110  124 
f = cos(M_PI * a / (double)(2 * n)); 
111  125 
sum += f * in[k]; 
112  126 
} 
113 
out[i] = sum;


127 
out[i] = REF_SCALE(sum, nbits  2);


114  128 
} 
115  129 
} 
116  130  
117  131 
/* NOTE: no normalisation by 1 / N is done */ 
118 
static void mdct_ref(float *output, float *input, int nbits)


132 
static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)


119  133 
{ 
120  134 
int n = 1<<nbits; 
121  135 
int k, i; 
...  ...  
128  142 
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); 
129  143 
s += input[i] * cos(a); 
130  144 
} 
131 
output[k] = s;


145 
output[k] = REF_SCALE(s, nbits  1);


132  146 
} 
133  147 
} 
134  148  
149 
#if CONFIG_FFT_FLOAT 

135  150 
static void idct_ref(float *output, float *input, int nbits) 
136  151 
{ 
137  152 
int n = 1<<nbits; 
...  ...  
164  179 
output[k] = s; 
165  180 
} 
166  181 
} 
182 
#endif 

167  183  
168  184  
169 
static float frandom(AVLFG *prng)


185 
static FFTSample frandom(AVLFG *prng)


170  186 
{ 
171 
return (int16_t)av_lfg_get(prng) / 32768.0; 

187 
return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;


172  188 
} 
173  189  
174  190 
static int64_t gettime(void) 
...  ...  
178  194 
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; 
179  195 
} 
180  196  
181 
static int check_diff(float *tab1, float *tab2, int n, double scale)


197 
static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)


182  198 
{ 
183  199 
int i; 
184  200 
double max= 0; 
...  ...  
186  202 
int err = 0; 
187  203  
188  204 
for (i = 0; i < n; i++) { 
189 
double e= fabsf(tab1[i]  (tab2[i] / scale));


205 
double e = fabsf(tab1[i]  (tab2[i] / scale)) / RANGE;


190  206 
if (e >= 1e3) { 
191 
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: %10.6f %10.6f\n",


207 
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",


192  208 
i, tab1[i], tab2[i]); 
193  209 
err = 1; 
194  210 
} 
...  ...  
233  249 
int do_inverse = 0; 
234  250 
FFTContext s1, *s = &s1; 
235  251 
FFTContext m1, *m = &m1; 
252 
#if CONFIG_FFT_FLOAT 

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

238  256 
int fft_nbits, fft_size, fft_size_2; 
239  257 
double scale = 1.0; 
240  258 
AVLFG prng; 
...  ...  
297  315 
ff_fft_init(s, fft_nbits, do_inverse); 
298  316 
fft_ref_init(fft_nbits, do_inverse); 
299  317 
break; 
318 
#if CONFIG_FFT_FLOAT 

300  319 
case TRANSFORM_RDFT: 
301  320 
if (do_inverse) 
302  321 
av_log(NULL, AV_LOG_INFO,"IDFT_C2R"); 
...  ...  
312  331 
av_log(NULL, AV_LOG_INFO,"DCT_II"); 
313  332 
ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II); 
314  333 
break; 
334 
#endif 

335 
default: 

336 
av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n"); 

337 
return 1; 

315  338 
} 
316  339 
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 
317  340  
...  ...  
328  351 
switch (transform) { 
329  352 
case TRANSFORM_MDCT: 
330  353 
if (do_inverse) { 
331 
imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits);


332 
m>imdct_calc(m, tab2, (float *)tab1);


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


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);


334  357 
} else { 
335 
mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits);


358 
mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);


336  359  
337 
m>mdct_calc(m, tab2, (float *)tab1);


360 
m>mdct_calc(m, tab2, (FFTSample *)tab1);


338  361  
339 
err = check_diff((float *)tab_ref, tab2, fft_size / 2, scale);


362 
err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);


340  363 
} 
341  364 
break; 
342  365 
case TRANSFORM_FFT: 
...  ...  
345  368 
s>fft_calc(s, tab); 
346  369  
347  370 
fft_ref(tab_ref, tab1, fft_nbits); 
348 
err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 1.0);


371 
err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);


349  372 
break; 
373 
#if CONFIG_FFT_FLOAT 

350  374 
case TRANSFORM_RDFT: 
351  375 
if (do_inverse) { 
352  376 
tab1[ 0].im = 0; 
...  ...  
387  411 
} 
388  412 
err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0); 
389  413 
break; 
414 
#endif 

390  415 
} 
391  416  
392  417 
/* do a speed test */ 
...  ...  
404  429 
switch (transform) { 
405  430 
case TRANSFORM_MDCT: 
406  431 
if (do_inverse) { 
407 
m>imdct_calc(m, (float *)tab, (float *)tab1);


432 
m>imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);


408  433 
} else { 
409 
m>mdct_calc(m, (float *)tab, (float *)tab1);


434 
m>mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);


410  435 
} 
411  436 
break; 
412  437 
case TRANSFORM_FFT: 
413  438 
memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); 
414  439 
s>fft_calc(s, tab); 
415  440 
break; 
441 
#if CONFIG_FFT_FLOAT 

416  442 
case TRANSFORM_RDFT: 
417  443 
memcpy(tab2, tab1, fft_size * sizeof(FFTSample)); 
418  444 
r>rdft_calc(r, tab2); 
...  ...  
421  447 
memcpy(tab2, tab1, fft_size * sizeof(FFTSample)); 
422  448 
d>dct_calc(d, tab2); 
423  449 
break; 
450 
#endif 

424  451 
} 
425  452 
} 
426  453 
duration = gettime()  time_start; 
...  ...  
441  468 
case TRANSFORM_FFT: 
442  469 
ff_fft_end(s); 
443  470 
break; 
471 
#if CONFIG_FFT_FLOAT 

444  472 
case TRANSFORM_RDFT: 
445  473 
ff_rdft_end(r); 
446  474 
break; 
447  475 
case TRANSFORM_DCT: 
448  476 
ff_dct_end(d); 
449  477 
break; 
478 
#endif 

450  479 
} 
451  480  
452  481 
av_free(tab); 
Also available in: Unified diff