ffmpeg / libavcodec / ffttest.c @ dc1ecb6a
History  View  Annotate  Download (7.5 KB)
1 
/*


2 
* (c) 2002 Fabrice Bellard

3 
*

4 
* This file is part of FFmpeg.

5 
*

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.

10 
*

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.

15 
*

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

19 
*/

20  
21 
/**

22 
* @file ffttest.c

23 
* FFT and MDCT tests.

24 
*/

25  
26 
#include "dsputil.h" 
27 
#include <math.h> 
28 
#include <unistd.h> 
29 
#include <sys/time.h> 
30 
#include <stdlib.h> 
31 
#include <string.h> 
32  
33 
#undef exit

34 
#undef random

35  
36 
int mm_flags;

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 
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 
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 
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 
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  
134 
float frandom(void) 
135 
{ 
136 
return (float)((random() & 0xffff)  32768) / 32768.0; 
137 
} 
138  
139 
int64_t gettime(void)

140 
{ 
141 
struct timeval tv;

142 
gettimeofday(&tv,NULL);

143 
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; 
144 
} 
145  
146 
void check_diff(float *tab1, float *tab2, int n) 
147 
{ 
148 
int i;

149 
double max= 0; 
150 
double error= 0; 
151  
152 
for(i=0;i<n;i++) { 
153 
double e= fabsf(tab1[i]  tab2[i]);

154 
if (e >= 1e3) { 
155 
av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n", 
156 
i, tab1[i], tab2[i]); 
157 
} 
158 
error+= e*e; 
159 
if(e>max) max= e;

160 
} 
161 
av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n); 
162 
} 
163  
164  
165 
void help(void) 
166 
{ 
167 
av_log(NULL, AV_LOG_INFO,"usage: ffttest [h] [s] [i] [n b]\n" 
168 
"h print this help\n"

169 
"s speed test\n"

170 
"m (I)MDCT test\n"

171 
"i inverse transform test\n"

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

173 
); 
174 
exit(1);

175 
} 
176  
177  
178  
179 
int main(int argc, char **argv) 
180 
{ 
181 
FFTComplex *tab, *tab1, *tab_ref; 
182 
FFTSample *tabtmp, *tab2; 
183 
int it, i, c;

184 
int do_speed = 0; 
185 
int do_mdct = 0; 
186 
int do_inverse = 0; 
187 
FFTContext s1, *s = &s1; 
188 
MDCTContext m1, *m = &m1; 
189 
int fft_nbits, fft_size;

190  
191 
mm_flags = 0;

192 
fft_nbits = 9;

193 
for(;;) {

194 
c = getopt(argc, argv, "hsimn:");

195 
if (c == 1) 
196 
break;

197 
switch(c) {

198 
case 'h': 
199 
help(); 
200 
break;

201 
case 's': 
202 
do_speed = 1;

203 
break;

204 
case 'i': 
205 
do_inverse = 1;

206 
break;

207 
case 'm': 
208 
do_mdct = 1;

209 
break;

210 
case 'n': 
211 
fft_nbits = atoi(optarg); 
212 
break;

213 
} 
214 
} 
215  
216 
fft_size = 1 << fft_nbits;

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

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

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

220 
tabtmp = av_malloc(fft_size / 2 * sizeof(FFTSample)); 
221 
tab2 = av_malloc(fft_size * sizeof(FFTSample));

222  
223 
if (do_mdct) {

224 
if (do_inverse)

225 
av_log(NULL, AV_LOG_INFO,"IMDCT"); 
226 
else

227 
av_log(NULL, AV_LOG_INFO,"MDCT"); 
228 
ff_mdct_init(m, fft_nbits, do_inverse); 
229 
} else {

230 
if (do_inverse)

231 
av_log(NULL, AV_LOG_INFO,"IFFT"); 
232 
else

233 
av_log(NULL, AV_LOG_INFO,"FFT"); 
234 
ff_fft_init(s, fft_nbits, do_inverse); 
235 
fft_ref_init(fft_nbits, do_inverse); 
236 
} 
237 
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 
238  
239 
/* generate random data */

240  
241 
for(i=0;i<fft_size;i++) { 
242 
tab1[i].re = frandom(); 
243 
tab1[i].im = frandom(); 
244 
} 
245  
246 
/* checking result */

247 
av_log(NULL, AV_LOG_INFO,"Checking...\n"); 
248  
249 
if (do_mdct) {

250 
if (do_inverse) {

251 
imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); 
252 
ff_imdct_calc(m, tab2, (float *)tab1, tabtmp);

253 
check_diff((float *)tab_ref, tab2, fft_size);

254 
} else {

255 
mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); 
256  
257 
ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);

258  
259 
check_diff((float *)tab_ref, tab2, fft_size / 2); 
260 
} 
261 
} else {

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

263 
ff_fft_permute(s, tab); 
264 
ff_fft_calc(s, tab); 
265  
266 
fft_ref(tab_ref, tab1, fft_nbits); 
267 
check_diff((float *)tab_ref, (float *)tab, fft_size * 2); 
268 
} 
269  
270 
/* do a speed test */

271  
272 
if (do_speed) {

273 
int64_t time_start, duration; 
274 
int nb_its;

275  
276 
av_log(NULL, AV_LOG_INFO,"Speed test...\n"); 
277 
/* we measure during about 1 seconds */

278 
nb_its = 1;

279 
for(;;) {

280 
time_start = gettime(); 
281 
for(it=0;it<nb_its;it++) { 
282 
if (do_mdct) {

283 
if (do_inverse) {

284 
ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp); 
285 
} else {

286 
ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp); 
287 
} 
288 
} else {

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

290 
ff_fft_calc(s, tab); 
291 
} 
292 
} 
293 
duration = gettime()  time_start; 
294 
if (duration >= 1000000) 
295 
break;

296 
nb_its *= 2;

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

300 
(double)duration / 1000000.0, 
301 
nb_its); 
302 
} 
303  
304 
if (do_mdct) {

305 
ff_mdct_end(m); 
306 
} else {

307 
ff_fft_end(s); 
308 
} 
309 
return 0; 
310 
} 