Add official LGPL license headers to the files that were missing them.
[ffmpeg.git] / libavcodec / fft-test.c
1 /*
2  * (c) 2002 Fabrice Bellard
3  *
4  * This library is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU Lesser General Public
6  * License as published by the Free Software Foundation; either
7  * version 2 of the License, or (at your option) any later version.
8  *
9  * This library is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12  * Lesser General Public License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with this library; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18
19 /**
20  * @file fft-test.c
21  * FFT and MDCT tests.
22  */
23
24 #include "dsputil.h"
25 #include <math.h>
26 #include <unistd.h>
27 #include <sys/time.h>
28
29 int mm_flags;
30
31 /* reference fft */
32
33 #define MUL16(a,b) ((a) * (b))
34
35 #define CMAC(pre, pim, are, aim, bre, bim) \
36 {\
37    pre += (MUL16(are, bre) - MUL16(aim, bim));\
38    pim += (MUL16(are, bim) + MUL16(bre, aim));\
39 }
40
41 FFTComplex *exptab;
42
43 void fft_ref_init(int nbits, int inverse)
44 {
45     int n, i;
46     float c1, s1, alpha;
47
48     n = 1 << nbits;
49     exptab = av_malloc((n / 2) * sizeof(FFTComplex));
50
51     for(i=0;i<(n/2);i++) {
52         alpha = 2 * M_PI * (float)i / (float)n;
53         c1 = cos(alpha);
54         s1 = sin(alpha);
55         if (!inverse)
56             s1 = -s1;
57         exptab[i].re = c1;
58         exptab[i].im = s1;
59     }
60 }
61
62 void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
63 {
64     int n, i, j, k, n2;
65     float tmp_re, tmp_im, s, c;
66     FFTComplex *q;
67
68     n = 1 << nbits;
69     n2 = n >> 1;
70     for(i=0;i<n;i++) {
71         tmp_re = 0;
72         tmp_im = 0;
73         q = tab;
74         for(j=0;j<n;j++) {
75             k = (i * j) & (n - 1);
76             if (k >= n2) {
77                 c = -exptab[k - n2].re;
78                 s = -exptab[k - n2].im;
79             } else {
80                 c = exptab[k].re;
81                 s = exptab[k].im;
82             }
83             CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
84             q++;
85         }
86         tabr[i].re = tmp_re;
87         tabr[i].im = tmp_im;
88     }
89 }
90
91 void imdct_ref(float *out, float *in, int n)
92 {
93     int k, i, a;
94     float sum, f;
95
96     for(i=0;i<n;i++) {
97         sum = 0;
98         for(k=0;k<n/2;k++) {
99             a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
100             f = cos(M_PI * a / (double)(2 * n));
101             sum += f * in[k];
102         }
103         out[i] = -sum;
104     }
105 }
106
107 /* NOTE: no normalisation by 1 / N is done */
108 void mdct_ref(float *output, float *input, int n)
109 {
110     int k, i;
111     float a, s;
112
113     /* do it by hand */
114     for(k=0;k<n/2;k++) {
115         s = 0;
116         for(i=0;i<n;i++) {
117             a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
118             s += input[i] * cos(a);
119         }
120         output[k] = s;
121     }
122 }
123
124
125 float frandom(void)
126 {
127     return (float)((random() & 0xffff) - 32768) / 32768.0;
128 }
129
130 int64_t gettime(void)
131 {
132     struct timeval tv;
133     gettimeofday(&tv,NULL);
134     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
135 }
136
137 void check_diff(float *tab1, float *tab2, int n)
138 {
139     int i;
140
141     for(i=0;i<n;i++) {
142         if (fabsf(tab1[i] - tab2[i]) >= 1e-3) {
143             av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n",
144                    i, tab1[i], tab2[i]);
145         }
146     }
147 }
148
149
150 void help(void)
151 {
152     av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
153            "-h     print this help\n"
154            "-s     speed test\n"
155            "-m     (I)MDCT test\n"
156            "-i     inverse transform test\n"
157            "-n b   set the transform size to 2^b\n"
158            );
159     exit(1);
160 }
161
162
163
164 int main(int argc, char **argv)
165 {
166     FFTComplex *tab, *tab1, *tab_ref;
167     FFTSample *tabtmp, *tab2;
168     int it, i, c;
169     int do_speed = 0;
170     int do_mdct = 0;
171     int do_inverse = 0;
172     FFTContext s1, *s = &s1;
173     MDCTContext m1, *m = &m1;
174     int fft_nbits, fft_size;
175
176     mm_flags = 0;
177     fft_nbits = 9;
178     for(;;) {
179         c = getopt(argc, argv, "hsimn:");
180         if (c == -1)
181             break;
182         switch(c) {
183         case 'h':
184             help();
185             break;
186         case 's':
187             do_speed = 1;
188             break;
189         case 'i':
190             do_inverse = 1;
191             break;
192         case 'm':
193             do_mdct = 1;
194             break;
195         case 'n':
196             fft_nbits = atoi(optarg);
197             break;
198         }
199     }
200
201     fft_size = 1 << fft_nbits;
202     tab = av_malloc(fft_size * sizeof(FFTComplex));
203     tab1 = av_malloc(fft_size * sizeof(FFTComplex));
204     tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
205     tabtmp = av_malloc(fft_size / 2 * sizeof(FFTSample));
206     tab2 = av_malloc(fft_size * sizeof(FFTSample));
207
208     if (do_mdct) {
209         if (do_inverse)
210             av_log(NULL, AV_LOG_INFO,"IMDCT");
211         else
212             av_log(NULL, AV_LOG_INFO,"MDCT");
213         ff_mdct_init(m, fft_nbits, do_inverse);
214     } else {
215         if (do_inverse)
216             av_log(NULL, AV_LOG_INFO,"IFFT");
217         else
218             av_log(NULL, AV_LOG_INFO,"FFT");
219         ff_fft_init(s, fft_nbits, do_inverse);
220         fft_ref_init(fft_nbits, do_inverse);
221     }
222     av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
223
224     /* generate random data */
225
226     for(i=0;i<fft_size;i++) {
227         tab1[i].re = frandom();
228         tab1[i].im = frandom();
229     }
230
231     /* checking result */
232     av_log(NULL, AV_LOG_INFO,"Checking...\n");
233
234     if (do_mdct) {
235         if (do_inverse) {
236             imdct_ref((float *)tab_ref, (float *)tab1, fft_size);
237             ff_imdct_calc(m, tab2, (float *)tab1, tabtmp);
238             check_diff((float *)tab_ref, tab2, fft_size);
239         } else {
240             mdct_ref((float *)tab_ref, (float *)tab1, fft_size);
241
242             ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);
243
244             check_diff((float *)tab_ref, tab2, fft_size / 2);
245         }
246     } else {
247         memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
248         ff_fft_permute(s, tab);
249         ff_fft_calc(s, tab);
250
251         fft_ref(tab_ref, tab1, fft_nbits);
252         check_diff((float *)tab_ref, (float *)tab, fft_size * 2);
253     }
254
255     /* do a speed test */
256
257     if (do_speed) {
258         int64_t time_start, duration;
259         int nb_its;
260
261         av_log(NULL, AV_LOG_INFO,"Speed test...\n");
262         /* we measure during about 1 seconds */
263         nb_its = 1;
264         for(;;) {
265             time_start = gettime();
266             for(it=0;it<nb_its;it++) {
267                 if (do_mdct) {
268                     if (do_inverse) {
269                         ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp);
270                     } else {
271                         ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp);
272                     }
273                 } else {
274                     memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
275                     ff_fft_calc(s, tab);
276                 }
277             }
278             duration = gettime() - time_start;
279             if (duration >= 1000000)
280                 break;
281             nb_its *= 2;
282         }
283         av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
284                (double)duration / nb_its,
285                (double)duration / 1000000.0,
286                nb_its);
287     }
288
289     if (do_mdct) {
290         ff_mdct_end(m);
291     } else {
292         ff_fft_end(s);
293     }
294     return 0;
295 }