ffmpeg / libavcodec / ffttest.c @ 94d85eaf
History  View  Annotate  Download (6.43 KB)
1 
/**


2 
* @file ffttest.c

3 
* FFT and MDCT tests.

4 
*/

5  
6 
#include "dsputil.h" 
7 
#include <math.h> 
8 
#include <unistd.h> 
9 
#include <sys/time.h> 
10  
11 
int mm_flags;

12  
13 
/* reference fft */

14  
15 
#define MUL16(a,b) ((a) * (b))

16  
17 
#define CMAC(pre, pim, are, aim, bre, bim) \

18 
{\ 
19 
pre += (MUL16(are, bre)  MUL16(aim, bim));\ 
20 
pim += (MUL16(are, bim) + MUL16(bre, aim));\ 
21 
} 
22  
23 
FFTComplex *exptab; 
24  
25 
void fft_ref_init(int nbits, int inverse) 
26 
{ 
27 
int n, i;

28 
float c1, s1, alpha;

29  
30 
n = 1 << nbits;

31 
exptab = av_malloc((n / 2) * sizeof(FFTComplex)); 
32  
33 
for(i=0;i<(n/2);i++) { 
34 
alpha = 2 * M_PI * (float)i / (float)n; 
35 
c1 = cos(alpha); 
36 
s1 = sin(alpha); 
37 
if (!inverse)

38 
s1 = s1; 
39 
exptab[i].re = c1; 
40 
exptab[i].im = s1; 
41 
} 
42 
} 
43  
44 
void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) 
45 
{ 
46 
int n, i, j, k, n2;

47 
float tmp_re, tmp_im, s, c;

48 
FFTComplex *q; 
49  
50 
n = 1 << nbits;

51 
n2 = n >> 1;

52 
for(i=0;i<n;i++) { 
53 
tmp_re = 0;

54 
tmp_im = 0;

55 
q = tab; 
56 
for(j=0;j<n;j++) { 
57 
k = (i * j) & (n  1);

58 
if (k >= n2) {

59 
c = exptab[k  n2].re; 
60 
s = exptab[k  n2].im; 
61 
} else {

62 
c = exptab[k].re; 
63 
s = exptab[k].im; 
64 
} 
65 
CMAC(tmp_re, tmp_im, c, s, q>re, q>im); 
66 
q++; 
67 
} 
68 
tabr[i].re = tmp_re; 
69 
tabr[i].im = tmp_im; 
70 
} 
71 
} 
72  
73 
void imdct_ref(float *out, float *in, int n) 
74 
{ 
75 
int k, i, a;

76 
float sum, f;

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

80 
for(k=0;k<n/2;k++) { 
81 
a = (2 * i + 1 + (n / 2)) * (2 * k + 1); 
82 
f = cos(M_PI * a / (double)(2 * n)); 
83 
sum += f * in[k]; 
84 
} 
85 
out[i] = sum; 
86 
} 
87 
} 
88  
89 
/* NOTE: no normalisation by 1 / N is done */

90 
void mdct_ref(float *output, float *input, int n) 
91 
{ 
92 
int k, i;

93 
float a, s;

94  
95 
/* do it by hand */

96 
for(k=0;k<n/2;k++) { 
97 
s = 0;

98 
for(i=0;i<n;i++) { 
99 
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); 
100 
s += input[i] * cos(a); 
101 
} 
102 
output[k] = s; 
103 
} 
104 
} 
105  
106  
107 
float frandom(void) 
108 
{ 
109 
return (float)((random() & 0xffff)  32768) / 32768.0; 
110 
} 
111  
112 
int64_t gettime(void)

113 
{ 
114 
struct timeval tv;

115 
gettimeofday(&tv,NULL);

116 
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; 
117 
} 
118  
119 
void check_diff(float *tab1, float *tab2, int n) 
120 
{ 
121 
int i;

122  
123 
for(i=0;i<n;i++) { 
124 
if (fabsf(tab1[i]  tab2[i]) >= 1e3) { 
125 
av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n", 
126 
i, tab1[i], tab2[i]); 
127 
} 
128 
} 
129 
} 
130  
131  
132 
void help(void) 
133 
{ 
134 
av_log(NULL, AV_LOG_INFO,"usage: ffttest [h] [s] [i] [n b]\n" 
135 
"h print this help\n"

136 
"s speed test\n"

137 
"m (I)MDCT test\n"

138 
"i inverse transform test\n"

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

140 
); 
141 
exit(1);

142 
} 
143  
144  
145  
146 
int main(int argc, char **argv) 
147 
{ 
148 
FFTComplex *tab, *tab1, *tab_ref; 
149 
FFTSample *tabtmp, *tab2; 
150 
int it, i, c;

151 
int do_speed = 0; 
152 
int do_mdct = 0; 
153 
int do_inverse = 0; 
154 
FFTContext s1, *s = &s1; 
155 
MDCTContext m1, *m = &m1; 
156 
int fft_nbits, fft_size;

157  
158 
mm_flags = 0;

159 
fft_nbits = 9;

160 
for(;;) {

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

162 
if (c == 1) 
163 
break;

164 
switch(c) {

165 
case 'h': 
166 
help(); 
167 
break;

168 
case 's': 
169 
do_speed = 1;

170 
break;

171 
case 'i': 
172 
do_inverse = 1;

173 
break;

174 
case 'm': 
175 
do_mdct = 1;

176 
break;

177 
case 'n': 
178 
fft_nbits = atoi(optarg); 
179 
break;

180 
} 
181 
} 
182  
183 
fft_size = 1 << fft_nbits;

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

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

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

187 
tabtmp = av_malloc(fft_size / 2 * sizeof(FFTSample)); 
188 
tab2 = av_malloc(fft_size * sizeof(FFTSample));

189  
190 
if (do_mdct) {

191 
if (do_inverse)

192 
av_log(NULL, AV_LOG_INFO,"IMDCT"); 
193 
else

194 
av_log(NULL, AV_LOG_INFO,"MDCT"); 
195 
ff_mdct_init(m, fft_nbits, do_inverse); 
196 
} else {

197 
if (do_inverse)

198 
av_log(NULL, AV_LOG_INFO,"IFFT"); 
199 
else

200 
av_log(NULL, AV_LOG_INFO,"FFT"); 
201 
ff_fft_init(s, fft_nbits, do_inverse); 
202 
fft_ref_init(fft_nbits, do_inverse); 
203 
} 
204 
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 
205  
206 
/* generate random data */

207  
208 
for(i=0;i<fft_size;i++) { 
209 
tab1[i].re = frandom(); 
210 
tab1[i].im = frandom(); 
211 
} 
212  
213 
/* checking result */

214 
av_log(NULL, AV_LOG_INFO,"Checking...\n"); 
215  
216 
if (do_mdct) {

217 
if (do_inverse) {

218 
imdct_ref((float *)tab_ref, (float *)tab1, fft_size); 
219 
ff_imdct_calc(m, tab2, (float *)tab1, tabtmp);

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

221 
} else {

222 
mdct_ref((float *)tab_ref, (float *)tab1, fft_size); 
223  
224 
ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);

225  
226 
check_diff((float *)tab_ref, tab2, fft_size / 2); 
227 
} 
228 
} else {

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

230 
ff_fft_permute(s, tab); 
231 
ff_fft_calc(s, tab); 
232  
233 
fft_ref(tab_ref, tab1, fft_nbits); 
234 
check_diff((float *)tab_ref, (float *)tab, fft_size * 2); 
235 
} 
236  
237 
/* do a speed test */

238  
239 
if (do_speed) {

240 
int64_t time_start, duration; 
241 
int nb_its;

242  
243 
av_log(NULL, AV_LOG_INFO,"Speed test...\n"); 
244 
/* we measure during about 1 seconds */

245 
nb_its = 1;

246 
for(;;) {

247 
time_start = gettime(); 
248 
for(it=0;it<nb_its;it++) { 
249 
if (do_mdct) {

250 
if (do_inverse) {

251 
ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp); 
252 
} else {

253 
ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp); 
254 
} 
255 
} else {

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

257 
ff_fft_calc(s, tab); 
258 
} 
259 
} 
260 
duration = gettime()  time_start; 
261 
if (duration >= 1000000) 
262 
break;

263 
nb_its *= 2;

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

267 
(double)duration / 1000000.0, 
268 
nb_its); 
269 
} 
270  
271 
if (do_mdct) {

272 
ff_mdct_end(m); 
273 
} else {

274 
ff_fft_end(s); 
275 
} 
276 
return 0; 
277 
} 