00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00026 #include "libavutil/cpu.h"
00027 #include "libavutil/mathematics.h"
00028 #include "libavutil/lfg.h"
00029 #include "libavutil/log.h"
00030 #include "libavutil/time.h"
00031 #include "fft.h"
00032 #if CONFIG_FFT_FLOAT
00033 #include "dct.h"
00034 #include "rdft.h"
00035 #endif
00036 #include <math.h>
00037 #if HAVE_UNISTD_H
00038 #include <unistd.h>
00039 #endif
00040 #include <stdlib.h>
00041 #include <string.h>
00042
00043
00044
00045 #define MUL16(a,b) ((a) * (b))
00046
00047 #define CMAC(pre, pim, are, aim, bre, bim) \
00048 {\
00049 pre += (MUL16(are, bre) - MUL16(aim, bim));\
00050 pim += (MUL16(are, bim) + MUL16(bre, aim));\
00051 }
00052
00053 #if CONFIG_FFT_FLOAT
00054 # define RANGE 1.0
00055 # define REF_SCALE(x, bits) (x)
00056 # define FMT "%10.6f"
00057 #else
00058 # define RANGE 16384
00059 # define REF_SCALE(x, bits) ((x) / (1<<(bits)))
00060 # define FMT "%6d"
00061 #endif
00062
00063 struct {
00064 float re, im;
00065 } *exptab;
00066
00067 static void fft_ref_init(int nbits, int inverse)
00068 {
00069 int n, i;
00070 double c1, s1, alpha;
00071
00072 n = 1 << nbits;
00073 exptab = av_malloc((n / 2) * sizeof(*exptab));
00074
00075 for (i = 0; i < (n/2); i++) {
00076 alpha = 2 * M_PI * (float)i / (float)n;
00077 c1 = cos(alpha);
00078 s1 = sin(alpha);
00079 if (!inverse)
00080 s1 = -s1;
00081 exptab[i].re = c1;
00082 exptab[i].im = s1;
00083 }
00084 }
00085
00086 static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
00087 {
00088 int n, i, j, k, n2;
00089 double tmp_re, tmp_im, s, c;
00090 FFTComplex *q;
00091
00092 n = 1 << nbits;
00093 n2 = n >> 1;
00094 for (i = 0; i < n; i++) {
00095 tmp_re = 0;
00096 tmp_im = 0;
00097 q = tab;
00098 for (j = 0; j < n; j++) {
00099 k = (i * j) & (n - 1);
00100 if (k >= n2) {
00101 c = -exptab[k - n2].re;
00102 s = -exptab[k - n2].im;
00103 } else {
00104 c = exptab[k].re;
00105 s = exptab[k].im;
00106 }
00107 CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
00108 q++;
00109 }
00110 tabr[i].re = REF_SCALE(tmp_re, nbits);
00111 tabr[i].im = REF_SCALE(tmp_im, nbits);
00112 }
00113 }
00114
00115 static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
00116 {
00117 int n = 1<<nbits;
00118 int k, i, a;
00119 double sum, f;
00120
00121 for (i = 0; i < n; i++) {
00122 sum = 0;
00123 for (k = 0; k < n/2; k++) {
00124 a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
00125 f = cos(M_PI * a / (double)(2 * n));
00126 sum += f * in[k];
00127 }
00128 out[i] = REF_SCALE(-sum, nbits - 2);
00129 }
00130 }
00131
00132
00133 static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
00134 {
00135 int n = 1<<nbits;
00136 int k, i;
00137 double a, s;
00138
00139
00140 for (k = 0; k < n/2; k++) {
00141 s = 0;
00142 for (i = 0; i < n; i++) {
00143 a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
00144 s += input[i] * cos(a);
00145 }
00146 output[k] = REF_SCALE(s, nbits - 1);
00147 }
00148 }
00149
00150 #if CONFIG_FFT_FLOAT
00151 static void idct_ref(FFTSample *output, FFTSample *input, int nbits)
00152 {
00153 int n = 1<<nbits;
00154 int k, i;
00155 double a, s;
00156
00157
00158 for (i = 0; i < n; i++) {
00159 s = 0.5 * input[0];
00160 for (k = 1; k < n; k++) {
00161 a = M_PI*k*(i+0.5) / n;
00162 s += input[k] * cos(a);
00163 }
00164 output[i] = 2 * s / n;
00165 }
00166 }
00167 static void dct_ref(FFTSample *output, FFTSample *input, int nbits)
00168 {
00169 int n = 1<<nbits;
00170 int k, i;
00171 double a, s;
00172
00173
00174 for (k = 0; k < n; k++) {
00175 s = 0;
00176 for (i = 0; i < n; i++) {
00177 a = M_PI*k*(i+0.5) / n;
00178 s += input[i] * cos(a);
00179 }
00180 output[k] = s;
00181 }
00182 }
00183 #endif
00184
00185
00186 static FFTSample frandom(AVLFG *prng)
00187 {
00188 return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;
00189 }
00190
00191 static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
00192 {
00193 int i;
00194 double max= 0;
00195 double error= 0;
00196 int err = 0;
00197
00198 for (i = 0; i < n; i++) {
00199 double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;
00200 if (e >= 1e-3) {
00201 av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
00202 i, tab1[i], tab2[i]);
00203 err = 1;
00204 }
00205 error+= e*e;
00206 if(e>max) max= e;
00207 }
00208 av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error/n));
00209 return err;
00210 }
00211
00212
00213 static void help(void)
00214 {
00215 av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
00216 "-h print this help\n"
00217 "-s speed test\n"
00218 "-m (I)MDCT test\n"
00219 "-d (I)DCT test\n"
00220 "-r (I)RDFT test\n"
00221 "-i inverse transform test\n"
00222 "-n b set the transform size to 2^b\n"
00223 "-f x set scale factor for output data of (I)MDCT to x\n"
00224 );
00225 }
00226
00227 enum tf_transform {
00228 TRANSFORM_FFT,
00229 TRANSFORM_MDCT,
00230 TRANSFORM_RDFT,
00231 TRANSFORM_DCT,
00232 };
00233
00234 #if !HAVE_GETOPT
00235 #include "compat/getopt.c"
00236 #endif
00237
00238 int main(int argc, char **argv)
00239 {
00240 FFTComplex *tab, *tab1, *tab_ref;
00241 FFTSample *tab2;
00242 int it, i, c;
00243 int cpuflags;
00244 int do_speed = 0;
00245 int err = 1;
00246 enum tf_transform transform = TRANSFORM_FFT;
00247 int do_inverse = 0;
00248 FFTContext s1, *s = &s1;
00249 FFTContext m1, *m = &m1;
00250 #if CONFIG_FFT_FLOAT
00251 RDFTContext r1, *r = &r1;
00252 DCTContext d1, *d = &d1;
00253 int fft_size_2;
00254 #endif
00255 int fft_nbits, fft_size;
00256 double scale = 1.0;
00257 AVLFG prng;
00258 av_lfg_init(&prng, 1);
00259
00260 fft_nbits = 9;
00261 for(;;) {
00262 c = getopt(argc, argv, "hsimrdn:f:c:");
00263 if (c == -1)
00264 break;
00265 switch(c) {
00266 case 'h':
00267 help();
00268 return 1;
00269 case 's':
00270 do_speed = 1;
00271 break;
00272 case 'i':
00273 do_inverse = 1;
00274 break;
00275 case 'm':
00276 transform = TRANSFORM_MDCT;
00277 break;
00278 case 'r':
00279 transform = TRANSFORM_RDFT;
00280 break;
00281 case 'd':
00282 transform = TRANSFORM_DCT;
00283 break;
00284 case 'n':
00285 fft_nbits = atoi(optarg);
00286 break;
00287 case 'f':
00288 scale = atof(optarg);
00289 break;
00290 case 'c':
00291 cpuflags = av_get_cpu_flags();
00292
00293 if (av_parse_cpu_caps(&cpuflags, optarg) < 0)
00294 return 1;
00295
00296 av_force_cpu_flags(cpuflags);
00297 break;
00298 }
00299 }
00300
00301 fft_size = 1 << fft_nbits;
00302 tab = av_malloc(fft_size * sizeof(FFTComplex));
00303 tab1 = av_malloc(fft_size * sizeof(FFTComplex));
00304 tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
00305 tab2 = av_malloc(fft_size * sizeof(FFTSample));
00306
00307 switch (transform) {
00308 case TRANSFORM_MDCT:
00309 av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale);
00310 if (do_inverse)
00311 av_log(NULL, AV_LOG_INFO,"IMDCT");
00312 else
00313 av_log(NULL, AV_LOG_INFO,"MDCT");
00314 ff_mdct_init(m, fft_nbits, do_inverse, scale);
00315 break;
00316 case TRANSFORM_FFT:
00317 if (do_inverse)
00318 av_log(NULL, AV_LOG_INFO,"IFFT");
00319 else
00320 av_log(NULL, AV_LOG_INFO,"FFT");
00321 ff_fft_init(s, fft_nbits, do_inverse);
00322 fft_ref_init(fft_nbits, do_inverse);
00323 break;
00324 #if CONFIG_FFT_FLOAT
00325 case TRANSFORM_RDFT:
00326 if (do_inverse)
00327 av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
00328 else
00329 av_log(NULL, AV_LOG_INFO,"DFT_R2C");
00330 ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
00331 fft_ref_init(fft_nbits, do_inverse);
00332 break;
00333 case TRANSFORM_DCT:
00334 if (do_inverse)
00335 av_log(NULL, AV_LOG_INFO,"DCT_III");
00336 else
00337 av_log(NULL, AV_LOG_INFO,"DCT_II");
00338 ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);
00339 break;
00340 #endif
00341 default:
00342 av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
00343 return 1;
00344 }
00345 av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
00346
00347
00348
00349 for (i = 0; i < fft_size; i++) {
00350 tab1[i].re = frandom(&prng);
00351 tab1[i].im = frandom(&prng);
00352 }
00353
00354
00355 av_log(NULL, AV_LOG_INFO,"Checking...\n");
00356
00357 switch (transform) {
00358 case TRANSFORM_MDCT:
00359 if (do_inverse) {
00360 imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
00361 m->imdct_calc(m, tab2, (FFTSample *)tab1);
00362 err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);
00363 } else {
00364 mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
00365
00366 m->mdct_calc(m, tab2, (FFTSample *)tab1);
00367
00368 err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);
00369 }
00370 break;
00371 case TRANSFORM_FFT:
00372 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
00373 s->fft_permute(s, tab);
00374 s->fft_calc(s, tab);
00375
00376 fft_ref(tab_ref, tab1, fft_nbits);
00377 err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);
00378 break;
00379 #if CONFIG_FFT_FLOAT
00380 case TRANSFORM_RDFT:
00381 fft_size_2 = fft_size >> 1;
00382 if (do_inverse) {
00383 tab1[ 0].im = 0;
00384 tab1[fft_size_2].im = 0;
00385 for (i = 1; i < fft_size_2; i++) {
00386 tab1[fft_size_2+i].re = tab1[fft_size_2-i].re;
00387 tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;
00388 }
00389
00390 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
00391 tab2[1] = tab1[fft_size_2].re;
00392
00393 r->rdft_calc(r, tab2);
00394 fft_ref(tab_ref, tab1, fft_nbits);
00395 for (i = 0; i < fft_size; i++) {
00396 tab[i].re = tab2[i];
00397 tab[i].im = 0;
00398 }
00399 err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5);
00400 } else {
00401 for (i = 0; i < fft_size; i++) {
00402 tab2[i] = tab1[i].re;
00403 tab1[i].im = 0;
00404 }
00405 r->rdft_calc(r, tab2);
00406 fft_ref(tab_ref, tab1, fft_nbits);
00407 tab_ref[0].im = tab_ref[fft_size_2].re;
00408 err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);
00409 }
00410 break;
00411 case TRANSFORM_DCT:
00412 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
00413 d->dct_calc(d, (FFTSample *)tab);
00414 if (do_inverse) {
00415 idct_ref((FFTSample*)tab_ref, (FFTSample *)tab1, fft_nbits);
00416 } else {
00417 dct_ref((FFTSample*)tab_ref, (FFTSample *)tab1, fft_nbits);
00418 }
00419 err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);
00420 break;
00421 #endif
00422 }
00423
00424
00425
00426 if (do_speed) {
00427 int64_t time_start, duration;
00428 int nb_its;
00429
00430 av_log(NULL, AV_LOG_INFO,"Speed test...\n");
00431
00432 nb_its = 1;
00433 for(;;) {
00434 time_start = av_gettime();
00435 for (it = 0; it < nb_its; it++) {
00436 switch (transform) {
00437 case TRANSFORM_MDCT:
00438 if (do_inverse) {
00439 m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
00440 } else {
00441 m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
00442 }
00443 break;
00444 case TRANSFORM_FFT:
00445 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
00446 s->fft_calc(s, tab);
00447 break;
00448 #if CONFIG_FFT_FLOAT
00449 case TRANSFORM_RDFT:
00450 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
00451 r->rdft_calc(r, tab2);
00452 break;
00453 case TRANSFORM_DCT:
00454 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
00455 d->dct_calc(d, tab2);
00456 break;
00457 #endif
00458 }
00459 }
00460 duration = av_gettime() - time_start;
00461 if (duration >= 1000000)
00462 break;
00463 nb_its *= 2;
00464 }
00465 av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
00466 (double)duration / nb_its,
00467 (double)duration / 1000000.0,
00468 nb_its);
00469 }
00470
00471 switch (transform) {
00472 case TRANSFORM_MDCT:
00473 ff_mdct_end(m);
00474 break;
00475 case TRANSFORM_FFT:
00476 ff_fft_end(s);
00477 break;
00478 #if CONFIG_FFT_FLOAT
00479 case TRANSFORM_RDFT:
00480 ff_rdft_end(r);
00481 break;
00482 case TRANSFORM_DCT:
00483 ff_dct_end(d);
00484 break;
00485 #endif
00486 }
00487
00488 av_free(tab);
00489 av_free(tab1);
00490 av_free(tab2);
00491 av_free(tab_ref);
00492 av_free(exptab);
00493
00494 return err;
00495 }
