fft-test.c
Go to the documentation of this file.
1 /*
2  * (c) 2002 Fabrice Bellard
3  *
4  * This file is part of Libav.
5  *
6  * Libav 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.
10  *
11  * Libav 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.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with Libav; if not, write to the Free Software
18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19  */
20 
26 #include "libavutil/mathematics.h"
27 #include "libavutil/lfg.h"
28 #include "libavutil/log.h"
29 #include "fft.h"
30 #if CONFIG_FFT_FLOAT
31 #include "dct.h"
32 #include "rdft.h"
33 #endif
34 #include <math.h>
35 #include <unistd.h>
36 #include <sys/time.h>
37 #include <stdlib.h>
38 #include <string.h>
39 
40 /* reference fft */
41 
42 #define MUL16(a,b) ((a) * (b))
43 
44 #define CMAC(pre, pim, are, aim, bre, bim) \
45 {\
46  pre += (MUL16(are, bre) - MUL16(aim, bim));\
47  pim += (MUL16(are, bim) + MUL16(bre, aim));\
48 }
49 
50 #if CONFIG_FFT_FLOAT
51 # define RANGE 1.0
52 # define REF_SCALE(x, bits) (x)
53 # define FMT "%10.6f"
54 #else
55 # define RANGE 16384
56 # define REF_SCALE(x, bits) ((x) / (1<<(bits)))
57 # define FMT "%6d"
58 #endif
59 
60 struct {
61  float re, im;
62 } *exptab;
63 
64 static void fft_ref_init(int nbits, int inverse)
65 {
66  int n, i;
67  double c1, s1, alpha;
68 
69  n = 1 << nbits;
70  exptab = av_malloc((n / 2) * sizeof(*exptab));
71 
72  for (i = 0; i < (n/2); i++) {
73  alpha = 2 * M_PI * (float)i / (float)n;
74  c1 = cos(alpha);
75  s1 = sin(alpha);
76  if (!inverse)
77  s1 = -s1;
78  exptab[i].re = c1;
79  exptab[i].im = s1;
80  }
81 }
82 
83 static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
84 {
85  int n, i, j, k, n2;
86  double tmp_re, tmp_im, s, c;
87  FFTComplex *q;
88 
89  n = 1 << nbits;
90  n2 = n >> 1;
91  for (i = 0; i < n; i++) {
92  tmp_re = 0;
93  tmp_im = 0;
94  q = tab;
95  for (j = 0; j < n; j++) {
96  k = (i * j) & (n - 1);
97  if (k >= n2) {
98  c = -exptab[k - n2].re;
99  s = -exptab[k - n2].im;
100  } else {
101  c = exptab[k].re;
102  s = exptab[k].im;
103  }
104  CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
105  q++;
106  }
107  tabr[i].re = REF_SCALE(tmp_re, nbits);
108  tabr[i].im = REF_SCALE(tmp_im, nbits);
109  }
110 }
111 
112 static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
113 {
114  int n = 1<<nbits;
115  int k, i, a;
116  double sum, f;
117 
118  for (i = 0; i < n; i++) {
119  sum = 0;
120  for (k = 0; k < n/2; k++) {
121  a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
122  f = cos(M_PI * a / (double)(2 * n));
123  sum += f * in[k];
124  }
125  out[i] = REF_SCALE(-sum, nbits - 2);
126  }
127 }
128 
129 /* NOTE: no normalisation by 1 / N is done */
130 static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
131 {
132  int n = 1<<nbits;
133  int k, i;
134  double a, s;
135 
136  /* do it by hand */
137  for (k = 0; k < n/2; k++) {
138  s = 0;
139  for (i = 0; i < n; i++) {
140  a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
141  s += input[i] * cos(a);
142  }
143  output[k] = REF_SCALE(s, nbits - 1);
144  }
145 }
146 
147 #if CONFIG_FFT_FLOAT
148 static void idct_ref(float *output, float *input, int nbits)
149 {
150  int n = 1<<nbits;
151  int k, i;
152  double a, s;
153 
154  /* do it by hand */
155  for (i = 0; i < n; i++) {
156  s = 0.5 * input[0];
157  for (k = 1; k < n; k++) {
158  a = M_PI*k*(i+0.5) / n;
159  s += input[k] * cos(a);
160  }
161  output[i] = 2 * s / n;
162  }
163 }
164 static void dct_ref(float *output, float *input, int nbits)
165 {
166  int n = 1<<nbits;
167  int k, i;
168  double a, s;
169 
170  /* do it by hand */
171  for (k = 0; k < n; k++) {
172  s = 0;
173  for (i = 0; i < n; i++) {
174  a = M_PI*k*(i+0.5) / n;
175  s += input[i] * cos(a);
176  }
177  output[k] = s;
178  }
179 }
180 #endif
181 
182 
183 static FFTSample frandom(AVLFG *prng)
184 {
185  return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;
186 }
187 
188 static int64_t gettime(void)
189 {
190  struct timeval tv;
191  gettimeofday(&tv,NULL);
192  return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
193 }
194 
195 static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
196 {
197  int i;
198  double max= 0;
199  double error= 0;
200  int err = 0;
201 
202  for (i = 0; i < n; i++) {
203  double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;
204  if (e >= 1e-3) {
205  av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
206  i, tab1[i], tab2[i]);
207  err = 1;
208  }
209  error+= e*e;
210  if(e>max) max= e;
211  }
212  av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n);
213  return err;
214 }
215 
216 
217 static void help(void)
218 {
219  av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
220  "-h print this help\n"
221  "-s speed test\n"
222  "-m (I)MDCT test\n"
223  "-d (I)DCT test\n"
224  "-r (I)RDFT test\n"
225  "-i inverse transform test\n"
226  "-n b set the transform size to 2^b\n"
227  "-f x set scale factor for output data of (I)MDCT to x\n"
228  );
229 }
230 
236 };
237 
238 int main(int argc, char **argv)
239 {
240  FFTComplex *tab, *tab1, *tab_ref;
241  FFTSample *tab2;
242  int it, i, c;
243  int do_speed = 0;
244  int err = 1;
245  enum tf_transform transform = TRANSFORM_FFT;
246  int do_inverse = 0;
247  FFTContext s1, *s = &s1;
248  FFTContext m1, *m = &m1;
249 #if CONFIG_FFT_FLOAT
250  RDFTContext r1, *r = &r1;
251  DCTContext d1, *d = &d1;
252  int fft_size_2;
253 #endif
254  int fft_nbits, fft_size;
255  double scale = 1.0;
256  AVLFG prng;
257  av_lfg_init(&prng, 1);
258 
259  fft_nbits = 9;
260  for(;;) {
261  c = getopt(argc, argv, "hsimrdn:f:");
262  if (c == -1)
263  break;
264  switch(c) {
265  case 'h':
266  help();
267  return 1;
268  case 's':
269  do_speed = 1;
270  break;
271  case 'i':
272  do_inverse = 1;
273  break;
274  case 'm':
275  transform = TRANSFORM_MDCT;
276  break;
277  case 'r':
278  transform = TRANSFORM_RDFT;
279  break;
280  case 'd':
281  transform = TRANSFORM_DCT;
282  break;
283  case 'n':
284  fft_nbits = atoi(optarg);
285  break;
286  case 'f':
287  scale = atof(optarg);
288  break;
289  }
290  }
291 
292  fft_size = 1 << fft_nbits;
293  tab = av_malloc(fft_size * sizeof(FFTComplex));
294  tab1 = av_malloc(fft_size * sizeof(FFTComplex));
295  tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
296  tab2 = av_malloc(fft_size * sizeof(FFTSample));
297 
298  switch (transform) {
299  case TRANSFORM_MDCT:
300  av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale);
301  if (do_inverse)
302  av_log(NULL, AV_LOG_INFO,"IMDCT");
303  else
304  av_log(NULL, AV_LOG_INFO,"MDCT");
305  ff_mdct_init(m, fft_nbits, do_inverse, scale);
306  break;
307  case TRANSFORM_FFT:
308  if (do_inverse)
309  av_log(NULL, AV_LOG_INFO,"IFFT");
310  else
311  av_log(NULL, AV_LOG_INFO,"FFT");
312  ff_fft_init(s, fft_nbits, do_inverse);
313  fft_ref_init(fft_nbits, do_inverse);
314  break;
315 #if CONFIG_FFT_FLOAT
316  case TRANSFORM_RDFT:
317  if (do_inverse)
318  av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
319  else
320  av_log(NULL, AV_LOG_INFO,"DFT_R2C");
321  ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);
322  fft_ref_init(fft_nbits, do_inverse);
323  break;
324  case TRANSFORM_DCT:
325  if (do_inverse)
326  av_log(NULL, AV_LOG_INFO,"DCT_III");
327  else
328  av_log(NULL, AV_LOG_INFO,"DCT_II");
329  ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);
330  break;
331 #endif
332  default:
333  av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
334  return 1;
335  }
336  av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
337 
338  /* generate random data */
339 
340  for (i = 0; i < fft_size; i++) {
341  tab1[i].re = frandom(&prng);
342  tab1[i].im = frandom(&prng);
343  }
344 
345  /* checking result */
346  av_log(NULL, AV_LOG_INFO,"Checking...\n");
347 
348  switch (transform) {
349  case TRANSFORM_MDCT:
350  if (do_inverse) {
351  imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
352  m->imdct_calc(m, tab2, (FFTSample *)tab1);
353  err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);
354  } else {
355  mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
356 
357  m->mdct_calc(m, tab2, (FFTSample *)tab1);
358 
359  err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);
360  }
361  break;
362  case TRANSFORM_FFT:
363  memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
364  s->fft_permute(s, tab);
365  s->fft_calc(s, tab);
366 
367  fft_ref(tab_ref, tab1, fft_nbits);
368  err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);
369  break;
370 #if CONFIG_FFT_FLOAT
371  case TRANSFORM_RDFT:
372  fft_size_2 = fft_size >> 1;
373  if (do_inverse) {
374  tab1[ 0].im = 0;
375  tab1[fft_size_2].im = 0;
376  for (i = 1; i < fft_size_2; i++) {
377  tab1[fft_size_2+i].re = tab1[fft_size_2-i].re;
378  tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;
379  }
380 
381  memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
382  tab2[1] = tab1[fft_size_2].re;
383 
384  r->rdft_calc(r, tab2);
385  fft_ref(tab_ref, tab1, fft_nbits);
386  for (i = 0; i < fft_size; i++) {
387  tab[i].re = tab2[i];
388  tab[i].im = 0;
389  }
390  err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5);
391  } else {
392  for (i = 0; i < fft_size; i++) {
393  tab2[i] = tab1[i].re;
394  tab1[i].im = 0;
395  }
396  r->rdft_calc(r, tab2);
397  fft_ref(tab_ref, tab1, fft_nbits);
398  tab_ref[0].im = tab_ref[fft_size_2].re;
399  err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);
400  }
401  break;
402  case TRANSFORM_DCT:
403  memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
404  d->dct_calc(d, tab);
405  if (do_inverse) {
406  idct_ref(tab_ref, tab1, fft_nbits);
407  } else {
408  dct_ref(tab_ref, tab1, fft_nbits);
409  }
410  err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);
411  break;
412 #endif
413  }
414 
415  /* do a speed test */
416 
417  if (do_speed) {
418  int64_t time_start, duration;
419  int nb_its;
420 
421  av_log(NULL, AV_LOG_INFO,"Speed test...\n");
422  /* we measure during about 1 seconds */
423  nb_its = 1;
424  for(;;) {
425  time_start = gettime();
426  for (it = 0; it < nb_its; it++) {
427  switch (transform) {
428  case TRANSFORM_MDCT:
429  if (do_inverse) {
430  m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
431  } else {
432  m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
433  }
434  break;
435  case TRANSFORM_FFT:
436  memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
437  s->fft_calc(s, tab);
438  break;
439 #if CONFIG_FFT_FLOAT
440  case TRANSFORM_RDFT:
441  memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
442  r->rdft_calc(r, tab2);
443  break;
444  case TRANSFORM_DCT:
445  memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
446  d->dct_calc(d, tab2);
447  break;
448 #endif
449  }
450  }
451  duration = gettime() - time_start;
452  if (duration >= 1000000)
453  break;
454  nb_its *= 2;
455  }
456  av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
457  (double)duration / nb_its,
458  (double)duration / 1000000.0,
459  nb_its);
460  }
461 
462  switch (transform) {
463  case TRANSFORM_MDCT:
464  ff_mdct_end(m);
465  break;
466  case TRANSFORM_FFT:
467  ff_fft_end(s);
468  break;
469 #if CONFIG_FFT_FLOAT
470  case TRANSFORM_RDFT:
471  ff_rdft_end(r);
472  break;
473  case TRANSFORM_DCT:
474  ff_dct_end(d);
475  break;
476 #endif
477  }
478 
479  av_free(tab);
480  av_free(tab1);
481  av_free(tab2);
482  av_free(tab_ref);
483  av_free(exptab);
484 
485  return err;
486 }