Mercurial > libavcodec.hg
annotate fft-test.c @ 1883:12f260596631 libavcodec
don't fetch the size of SEQH unless SEQH was found
author | melanson |
---|---|
date | Sun, 14 Mar 2004 15:04:21 +0000 |
parents | dd63cb7e5080 |
children | 786ccf72ccd5 |
rev | line source |
---|---|
1106 | 1 /** |
2 * @file fft-test.c | |
3 * FFT and MDCT tests. | |
4 */ | |
5 | |
781 | 6 #include "dsputil.h" |
7 #include <math.h> | |
980 | 8 #include <unistd.h> |
781 | 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 | |
1064 | 112 int64_t gettime(void) |
781 | 113 { |
114 struct timeval tv; | |
115 gettimeofday(&tv,NULL); | |
1064 | 116 return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; |
781 | 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]) >= 1e-3) { | |
125 printf("ERROR %d: %f %f\n", | |
126 i, tab1[i], tab2[i]); | |
127 } | |
128 } | |
129 } | |
130 | |
131 | |
132 void help(void) | |
133 { | |
134 printf("usage: fft-test [-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 printf("IMDCT"); | |
193 else | |
194 printf("MDCT"); | |
969 | 195 ff_mdct_init(m, fft_nbits, do_inverse); |
781 | 196 } else { |
197 if (do_inverse) | |
198 printf("IFFT"); | |
199 else | |
200 printf("FFT"); | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
201 ff_fft_init(s, fft_nbits, do_inverse); |
781 | 202 fft_ref_init(fft_nbits, do_inverse); |
203 } | |
204 printf(" %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 printf("Checking...\n"); | |
215 | |
216 if (do_mdct) { | |
217 if (do_inverse) { | |
218 imdct_ref((float *)tab_ref, (float *)tab1, fft_size); | |
969 | 219 ff_imdct_calc(m, tab2, (float *)tab1, tabtmp); |
781 | 220 check_diff((float *)tab_ref, tab2, fft_size); |
221 } else { | |
222 mdct_ref((float *)tab_ref, (float *)tab1, fft_size); | |
223 | |
969 | 224 ff_mdct_calc(m, tab2, (float *)tab1, tabtmp); |
781 | 225 |
226 check_diff((float *)tab_ref, tab2, fft_size / 2); | |
227 } | |
228 } else { | |
229 memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
230 ff_fft_permute(s, tab); |
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
231 ff_fft_calc(s, tab); |
781 | 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) { | |
1064 | 240 int64_t time_start, duration; |
781 | 241 int nb_its; |
242 | |
243 printf("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) { | |
969 | 251 ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp); |
781 | 252 } else { |
969 | 253 ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp); |
781 | 254 } |
255 } else { | |
256 memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); | |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
257 ff_fft_calc(s, tab); |
781 | 258 } |
259 } | |
260 duration = gettime() - time_start; | |
261 if (duration >= 1000000) | |
262 break; | |
263 nb_its *= 2; | |
264 } | |
265 printf("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) { | |
969 | 272 ff_mdct_end(m); |
781 | 273 } else { |
1879
dd63cb7e5080
fft_*() renamed into ff_fft_*() patch by (Gildas Bazin <gbazin at altern dot org>)
michael
parents:
1106
diff
changeset
|
274 ff_fft_end(s); |
781 | 275 } |
276 return 0; | |
277 } |