2 * (c) 2002 Fabrice Bellard
4 * This file is part of FFmpeg.
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.
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.
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 02110-1301 USA
22 * @file libavcodec/fft-test.c
26 #include "libavutil/lfg.h"
38 #define MUL16(a,b) ((a) * (b))
40 #define CMAC(pre, pim, are, aim, bre, bim) \
42 pre += (MUL16(are, bre) - MUL16(aim, bim));\
43 pim += (MUL16(are, bim) + MUL16(bre, aim));\
48 static void fft_ref_init(int nbits
, int inverse
)
54 exptab
= av_malloc((n
/ 2) * sizeof(FFTComplex
));
56 for(i
=0;i
<(n
/2);i
++) {
57 alpha
= 2 * M_PI
* (float)i
/ (float)n
;
67 static void fft_ref(FFTComplex
*tabr
, FFTComplex
*tab
, int nbits
)
70 double tmp_re
, tmp_im
, s
, c
;
80 k
= (i
* j
) & (n
- 1);
82 c
= -exptab
[k
- n2
].re
;
83 s
= -exptab
[k
- n2
].im
;
88 CMAC(tmp_re
, tmp_im
, c
, s
, q
->re
, q
->im
);
96 static void imdct_ref(float *out
, float *in
, int nbits
)
105 a
= (2 * i
+ 1 + (n
/ 2)) * (2 * k
+ 1);
106 f
= cos(M_PI
* a
/ (double)(2 * n
));
113 /* NOTE: no normalisation by 1 / N is done */
114 static void mdct_ref(float *output
, float *input
, int nbits
)
124 a
= (2*M_PI
*(2*i
+1+n
/2)*(2*k
+1) / (4 * n
));
125 s
+= input
[i
] * cos(a
);
132 static float frandom(void)
135 av_lfg_init(&prng
, 1);
136 return (float)((av_lfg_get(&prng
) & 0xffff) - 32768) / 32768.0;
139 static int64_t gettime(void)
142 gettimeofday(&tv
,NULL
);
143 return (int64_t)tv
.tv_sec
* 1000000 + tv
.tv_usec
;
146 static void check_diff(float *tab1
, float *tab2
, int n
)
153 double e
= fabsf(tab1
[i
] - tab2
[i
]);
155 av_log(NULL
, AV_LOG_ERROR
, "ERROR %d: %f %f\n",
156 i
, tab1
[i
], tab2
[i
]);
161 av_log(NULL
, AV_LOG_INFO
, "max:%f e:%g\n", max
, sqrt(error
)/n
);
165 static void help(void)
167 av_log(NULL
, AV_LOG_INFO
,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
168 "-h print this help\n"
171 "-i inverse transform test\n"
172 "-n b set the transform size to 2^b\n"
179 int main(int argc
, char **argv
)
181 FFTComplex
*tab
, *tab1
, *tab_ref
;
187 FFTContext s1
, *s
= &s1
;
188 MDCTContext m1
, *m
= &m1
;
189 int fft_nbits
, fft_size
;
193 c
= getopt(argc
, argv
, "hsimn:");
210 fft_nbits
= atoi(optarg
);
215 fft_size
= 1 << fft_nbits
;
216 tab
= av_malloc(fft_size
* sizeof(FFTComplex
));
217 tab1
= av_malloc(fft_size
* sizeof(FFTComplex
));
218 tab_ref
= av_malloc(fft_size
* sizeof(FFTComplex
));
219 tab2
= av_malloc(fft_size
* sizeof(FFTSample
));
223 av_log(NULL
, AV_LOG_INFO
,"IMDCT");
225 av_log(NULL
, AV_LOG_INFO
,"MDCT");
226 ff_mdct_init(m
, fft_nbits
, do_inverse
);
229 av_log(NULL
, AV_LOG_INFO
,"IFFT");
231 av_log(NULL
, AV_LOG_INFO
,"FFT");
232 ff_fft_init(s
, fft_nbits
, do_inverse
);
233 fft_ref_init(fft_nbits
, do_inverse
);
235 av_log(NULL
, AV_LOG_INFO
," %d test\n", fft_size
);
237 /* generate random data */
239 for(i
=0;i
<fft_size
;i
++) {
240 tab1
[i
].re
= frandom();
241 tab1
[i
].im
= frandom();
244 /* checking result */
245 av_log(NULL
, AV_LOG_INFO
,"Checking...\n");
249 imdct_ref((float *)tab_ref
, (float *)tab1
, fft_nbits
);
250 ff_imdct_calc(m
, tab2
, (float *)tab1
);
251 check_diff((float *)tab_ref
, tab2
, fft_size
);
253 mdct_ref((float *)tab_ref
, (float *)tab1
, fft_nbits
);
255 ff_mdct_calc(m
, tab2
, (float *)tab1
);
257 check_diff((float *)tab_ref
, tab2
, fft_size
/ 2);
260 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
261 ff_fft_permute(s
, tab
);
264 fft_ref(tab_ref
, tab1
, fft_nbits
);
265 check_diff((float *)tab_ref
, (float *)tab
, fft_size
* 2);
268 /* do a speed test */
271 int64_t time_start
, duration
;
274 av_log(NULL
, AV_LOG_INFO
,"Speed test...\n");
275 /* we measure during about 1 seconds */
278 time_start
= gettime();
279 for(it
=0;it
<nb_its
;it
++) {
282 ff_imdct_calc(m
, (float *)tab
, (float *)tab1
);
284 ff_mdct_calc(m
, (float *)tab
, (float *)tab1
);
287 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
291 duration
= gettime() - time_start
;
292 if (duration
>= 1000000)
296 av_log(NULL
, AV_LOG_INFO
,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
297 (double)duration
/ nb_its
,
298 (double)duration
/ 1000000.0,