Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

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