ffmpeg / libavcodec / ffttest.c @ cc61f96f
History  View  Annotate  Download (7.28 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 
float 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 
float 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 n) 
99 
{ 
100 
int k, i, a;

101 
float sum, f;

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

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

115 
void mdct_ref(float *output, float *input, int n) 
116 
{ 
117 
int k, i;

118 
float a, s;

119  
120 
/* do it by hand */

121 
for(k=0;k<n/2;k++) { 
122 
s = 0;

123 
for(i=0;i<n;i++) { 
124 
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); 
125 
s += input[i] * cos(a); 
126 
} 
127 
output[k] = s; 
128 
} 
129 
} 
130  
131  
132 
float frandom(void) 
133 
{ 
134 
return (float)((random() & 0xffff)  32768) / 32768.0; 
135 
} 
136  
137 
int64_t gettime(void)

138 
{ 
139 
struct timeval tv;

140 
gettimeofday(&tv,NULL);

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

147  
148 
for(i=0;i<n;i++) { 
149 
if (fabsf(tab1[i]  tab2[i]) >= 1e3) { 
150 
av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n", 
151 
i, tab1[i], tab2[i]); 
152 
} 
153 
} 
154 
} 
155  
156  
157 
void help(void) 
158 
{ 
159 
av_log(NULL, AV_LOG_INFO,"usage: ffttest [h] [s] [i] [n b]\n" 
160 
"h print this help\n"

161 
"s speed test\n"

162 
"m (I)MDCT test\n"

163 
"i inverse transform test\n"

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

165 
); 
166 
exit(1);

167 
} 
168  
169  
170  
171 
int main(int argc, char **argv) 
172 
{ 
173 
FFTComplex *tab, *tab1, *tab_ref; 
174 
FFTSample *tabtmp, *tab2; 
175 
int it, i, c;

176 
int do_speed = 0; 
177 
int do_mdct = 0; 
178 
int do_inverse = 0; 
179 
FFTContext s1, *s = &s1; 
180 
MDCTContext m1, *m = &m1; 
181 
int fft_nbits, fft_size;

182  
183 
mm_flags = 0;

184 
fft_nbits = 9;

185 
for(;;) {

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

187 
if (c == 1) 
188 
break;

189 
switch(c) {

190 
case 'h': 
191 
help(); 
192 
break;

193 
case 's': 
194 
do_speed = 1;

195 
break;

196 
case 'i': 
197 
do_inverse = 1;

198 
break;

199 
case 'm': 
200 
do_mdct = 1;

201 
break;

202 
case 'n': 
203 
fft_nbits = atoi(optarg); 
204 
break;

205 
} 
206 
} 
207  
208 
fft_size = 1 << fft_nbits;

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

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

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

212 
tabtmp = av_malloc(fft_size / 2 * sizeof(FFTSample)); 
213 
tab2 = av_malloc(fft_size * sizeof(FFTSample));

214  
215 
if (do_mdct) {

216 
if (do_inverse)

217 
av_log(NULL, AV_LOG_INFO,"IMDCT"); 
218 
else

219 
av_log(NULL, AV_LOG_INFO,"MDCT"); 
220 
ff_mdct_init(m, fft_nbits, do_inverse); 
221 
} else {

222 
if (do_inverse)

223 
av_log(NULL, AV_LOG_INFO,"IFFT"); 
224 
else

225 
av_log(NULL, AV_LOG_INFO,"FFT"); 
226 
ff_fft_init(s, fft_nbits, do_inverse); 
227 
fft_ref_init(fft_nbits, do_inverse); 
228 
} 
229 
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 
230  
231 
/* generate random data */

232  
233 
for(i=0;i<fft_size;i++) { 
234 
tab1[i].re = frandom(); 
235 
tab1[i].im = frandom(); 
236 
} 
237  
238 
/* checking result */

239 
av_log(NULL, AV_LOG_INFO,"Checking...\n"); 
240  
241 
if (do_mdct) {

242 
if (do_inverse) {

243 
imdct_ref((float *)tab_ref, (float *)tab1, fft_size); 
244 
ff_imdct_calc(m, tab2, (float *)tab1, tabtmp);

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

246 
} else {

247 
mdct_ref((float *)tab_ref, (float *)tab1, fft_size); 
248  
249 
ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);

250  
251 
check_diff((float *)tab_ref, tab2, fft_size / 2); 
252 
} 
253 
} else {

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

255 
ff_fft_permute(s, tab); 
256 
ff_fft_calc(s, tab); 
257  
258 
fft_ref(tab_ref, tab1, fft_nbits); 
259 
check_diff((float *)tab_ref, (float *)tab, fft_size * 2); 
260 
} 
261  
262 
/* do a speed test */

263  
264 
if (do_speed) {

265 
int64_t time_start, duration; 
266 
int nb_its;

267  
268 
av_log(NULL, AV_LOG_INFO,"Speed test...\n"); 
269 
/* we measure during about 1 seconds */

270 
nb_its = 1;

271 
for(;;) {

272 
time_start = gettime(); 
273 
for(it=0;it<nb_its;it++) { 
274 
if (do_mdct) {

275 
if (do_inverse) {

276 
ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp); 
277 
} else {

278 
ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp); 
279 
} 
280 
} else {

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

282 
ff_fft_calc(s, tab); 
283 
} 
284 
} 
285 
duration = gettime()  time_start; 
286 
if (duration >= 1000000) 
287 
break;

288 
nb_its *= 2;

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

292 
(double)duration / 1000000.0, 
293 
nb_its); 
294 
} 
295  
296 
if (do_mdct) {

297 
ff_mdct_end(m); 
298 
} else {

299 
ff_fft_end(s); 
300 
} 
301 
return 0; 
302 
} 