Mercurial > libavcodec.hg
comparison fft-test.c @ 781:6f5e87957bcb libavcodec
new generic FFT/MDCT code for audio codecs
author | bellard |
---|---|
date | Mon, 28 Oct 2002 00:34:08 +0000 |
parents | |
children | d12d0c93d1e2 |
comparison
equal
deleted
inserted
replaced
780:a48bb8bc63dd | 781:6f5e87957bcb |
---|---|
1 /* FFT and MDCT tests */ | |
2 #include "dsputil.h" | |
3 #include <math.h> | |
4 #include <getopt.h> | |
5 #include <sys/time.h> | |
6 | |
7 int mm_flags; | |
8 | |
9 void *av_malloc(int size) | |
10 { | |
11 void *ptr; | |
12 ptr = malloc(size); | |
13 return ptr; | |
14 } | |
15 | |
16 void av_free(void *ptr) | |
17 { | |
18 /* XXX: this test should not be needed on most libcs */ | |
19 if (ptr) | |
20 free(ptr); | |
21 } | |
22 | |
23 /* cannot call it directly because of 'void **' casting is not automatic */ | |
24 void __av_freep(void **ptr) | |
25 { | |
26 av_free(*ptr); | |
27 *ptr = NULL; | |
28 } | |
29 | |
30 /* reference fft */ | |
31 | |
32 #define MUL16(a,b) ((a) * (b)) | |
33 | |
34 #define CMAC(pre, pim, are, aim, bre, bim) \ | |
35 {\ | |
36 pre += (MUL16(are, bre) - MUL16(aim, bim));\ | |
37 pim += (MUL16(are, bim) + MUL16(bre, aim));\ | |
38 } | |
39 | |
40 FFTComplex *exptab; | |
41 | |
42 void fft_ref_init(int nbits, int inverse) | |
43 { | |
44 int n, i; | |
45 float c1, s1, alpha; | |
46 | |
47 n = 1 << nbits; | |
48 exptab = av_malloc((n / 2) * sizeof(FFTComplex)); | |
49 | |
50 for(i=0;i<(n/2);i++) { | |
51 alpha = 2 * M_PI * (float)i / (float)n; | |
52 c1 = cos(alpha); | |
53 s1 = sin(alpha); | |
54 if (!inverse) | |
55 s1 = -s1; | |
56 exptab[i].re = c1; | |
57 exptab[i].im = s1; | |
58 } | |
59 } | |
60 | |
61 void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) | |
62 { | |
63 int n, i, j, k, n2; | |
64 float tmp_re, tmp_im, s, c; | |
65 FFTComplex *q; | |
66 | |
67 n = 1 << nbits; | |
68 n2 = n >> 1; | |
69 for(i=0;i<n;i++) { | |
70 tmp_re = 0; | |
71 tmp_im = 0; | |
72 q = tab; | |
73 for(j=0;j<n;j++) { | |
74 k = (i * j) & (n - 1); | |
75 if (k >= n2) { | |
76 c = -exptab[k - n2].re; | |
77 s = -exptab[k - n2].im; | |
78 } else { | |
79 c = exptab[k].re; | |
80 s = exptab[k].im; | |
81 } | |
82 CMAC(tmp_re, tmp_im, c, s, q->re, q->im); | |
83 q++; | |
84 } | |
85 tabr[i].re = tmp_re; | |
86 tabr[i].im = tmp_im; | |
87 } | |
88 } | |
89 | |
90 void imdct_ref(float *out, float *in, int n) | |
91 { | |
92 int k, i, a; | |
93 float sum, f; | |
94 | |
95 for(i=0;i<n;i++) { | |
96 sum = 0; | |
97 for(k=0;k<n/2;k++) { | |
98 a = (2 * i + 1 + (n / 2)) * (2 * k + 1); | |
99 f = cos(M_PI * a / (double)(2 * n)); | |
100 sum += f * in[k]; | |
101 } | |
102 out[i] = -sum; | |
103 } | |
104 } | |
105 | |
106 /* NOTE: no normalisation by 1 / N is done */ | |
107 void mdct_ref(float *output, float *input, int n) | |
108 { | |
109 int k, i; | |
110 float a, s; | |
111 | |
112 /* do it by hand */ | |
113 for(k=0;k<n/2;k++) { | |
114 s = 0; | |
115 for(i=0;i<n;i++) { | |
116 a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); | |
117 s += input[i] * cos(a); | |
118 } | |
119 output[k] = s; | |
120 } | |
121 } | |
122 | |
123 | |
124 float frandom(void) | |
125 { | |
126 return (float)((random() & 0xffff) - 32768) / 32768.0; | |
127 } | |
128 | |
129 INT64 gettime(void) | |
130 { | |
131 struct timeval tv; | |
132 gettimeofday(&tv,NULL); | |
133 return (INT64)tv.tv_sec * 1000000 + tv.tv_usec; | |
134 } | |
135 | |
136 void check_diff(float *tab1, float *tab2, int n) | |
137 { | |
138 int i; | |
139 | |
140 for(i=0;i<n;i++) { | |
141 if (fabsf(tab1[i] - tab2[i]) >= 1e-3) { | |
142 printf("ERROR %d: %f %f\n", | |
143 i, tab1[i], tab2[i]); | |
144 } | |
145 } | |
146 } | |
147 | |
148 | |
149 void help(void) | |
150 { | |
151 printf("usage: fft-test [-h] [-s] [-i] [-n b]\n" | |
152 "-h print this help\n" | |
153 "-s speed test\n" | |
154 "-m (I)MDCT test\n" | |
155 "-i inverse transform test\n" | |
156 "-n b set the transform size to 2^b\n" | |
157 ); | |
158 exit(1); | |
159 } | |
160 | |
161 | |
162 | |
163 int main(int argc, char **argv) | |
164 { | |
165 FFTComplex *tab, *tab1, *tab_ref; | |
166 FFTSample *tabtmp, *tab2; | |
167 int it, i, c; | |
168 int do_speed = 0; | |
169 int do_mdct = 0; | |
170 int do_inverse = 0; | |
171 FFTContext s1, *s = &s1; | |
172 MDCTContext m1, *m = &m1; | |
173 int fft_nbits, fft_size; | |
174 | |
175 mm_flags = 0; | |
176 fft_nbits = 9; | |
177 for(;;) { | |
178 c = getopt(argc, argv, "hsimn:"); | |
179 if (c == -1) | |
180 break; | |
181 switch(c) { | |
182 case 'h': | |
183 help(); | |
184 break; | |
185 case 's': | |
186 do_speed = 1; | |
187 break; | |
188 case 'i': | |
189 do_inverse = 1; | |
190 break; | |
191 case 'm': | |
192 do_mdct = 1; | |
193 break; | |
194 case 'n': | |
195 fft_nbits = atoi(optarg); | |
196 break; | |
197 } | |
198 } | |
199 | |
200 fft_size = 1 << fft_nbits; | |
201 tab = av_malloc(fft_size * sizeof(FFTComplex)); | |
202 tab1 = av_malloc(fft_size * sizeof(FFTComplex)); | |
203 tab_ref = av_malloc(fft_size * sizeof(FFTComplex)); | |
204 tabtmp = av_malloc(fft_size / 2 * sizeof(FFTSample)); | |
205 tab2 = av_malloc(fft_size * sizeof(FFTSample)); | |
206 | |
207 if (do_mdct) { | |
208 if (do_inverse) | |
209 printf("IMDCT"); | |
210 else | |
211 printf("MDCT"); | |
212 mdct_init(m, fft_nbits, do_inverse); | |
213 } else { | |
214 if (do_inverse) | |
215 printf("IFFT"); | |
216 else | |
217 printf("FFT"); | |
218 fft_init(s, fft_nbits, do_inverse); | |
219 fft_ref_init(fft_nbits, do_inverse); | |
220 } | |
221 printf(" %d test\n", fft_size); | |
222 | |
223 /* generate random data */ | |
224 | |
225 for(i=0;i<fft_size;i++) { | |
226 tab1[i].re = frandom(); | |
227 tab1[i].im = frandom(); | |
228 } | |
229 | |
230 /* checking result */ | |
231 printf("Checking...\n"); | |
232 | |
233 if (do_mdct) { | |
234 if (do_inverse) { | |
235 imdct_ref((float *)tab_ref, (float *)tab1, fft_size); | |
236 imdct_calc(m, tab2, (float *)tab1, tabtmp); | |
237 check_diff((float *)tab_ref, tab2, fft_size); | |
238 } else { | |
239 mdct_ref((float *)tab_ref, (float *)tab1, fft_size); | |
240 | |
241 mdct_calc(m, tab2, (float *)tab1, tabtmp); | |
242 | |
243 check_diff((float *)tab_ref, tab2, fft_size / 2); | |
244 } | |
245 } else { | |
246 memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); | |
247 fft_permute(s, tab); | |
248 fft_calc(s, tab); | |
249 | |
250 fft_ref(tab_ref, tab1, fft_nbits); | |
251 check_diff((float *)tab_ref, (float *)tab, fft_size * 2); | |
252 } | |
253 | |
254 /* do a speed test */ | |
255 | |
256 if (do_speed) { | |
257 INT64 time_start, duration; | |
258 int nb_its; | |
259 | |
260 printf("Speed test...\n"); | |
261 /* we measure during about 1 seconds */ | |
262 nb_its = 1; | |
263 for(;;) { | |
264 time_start = gettime(); | |
265 for(it=0;it<nb_its;it++) { | |
266 if (do_mdct) { | |
267 if (do_inverse) { | |
268 imdct_calc(m, (float *)tab, (float *)tab1, tabtmp); | |
269 } else { | |
270 mdct_calc(m, (float *)tab, (float *)tab1, tabtmp); | |
271 } | |
272 } else { | |
273 memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); | |
274 fft_calc(s, tab); | |
275 } | |
276 } | |
277 duration = gettime() - time_start; | |
278 if (duration >= 1000000) | |
279 break; | |
280 nb_its *= 2; | |
281 } | |
282 printf("time: %0.1f us/transform [total time=%0.2f s its=%d]\n", | |
283 (double)duration / nb_its, | |
284 (double)duration / 1000000.0, | |
285 nb_its); | |
286 } | |
287 | |
288 if (do_mdct) { | |
289 mdct_end(m); | |
290 } else { | |
291 fft_end(s); | |
292 } | |
293 return 0; | |
294 } |