• Main Page
  • Related Pages
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

libavcodec/fft-test.c

Go to the documentation of this file.
00001 /*
00002  * (c) 2002 Fabrice Bellard
00003  *
00004  * This file is part of FFmpeg.
00005  *
00006  * FFmpeg is free software; you can redistribute it and/or
00007  * modify it under the terms of the GNU Lesser General Public
00008  * License as published by the Free Software Foundation; either
00009  * version 2.1 of the License, or (at your option) any later version.
00010  *
00011  * FFmpeg is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with FFmpeg; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
00019  */
00020 
00026 #include "dsputil.h"
00027 #include <math.h>
00028 #include <unistd.h>
00029 #include <sys/time.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 
00033 #undef exit
00034 #undef random
00035 
00036 /* reference fft */
00037 
00038 #define MUL16(a,b) ((a) * (b))
00039 
00040 #define CMAC(pre, pim, are, aim, bre, bim) \
00041 {\
00042    pre += (MUL16(are, bre) - MUL16(aim, bim));\
00043    pim += (MUL16(are, bim) + MUL16(bre, aim));\
00044 }
00045 
00046 FFTComplex *exptab;
00047 
00048 void fft_ref_init(int nbits, int inverse)
00049 {
00050     int n, i;
00051     double c1, s1, alpha;
00052 
00053     n = 1 << nbits;
00054     exptab = av_malloc((n / 2) * sizeof(FFTComplex));
00055 
00056     for(i=0;i<(n/2);i++) {
00057         alpha = 2 * M_PI * (float)i / (float)n;
00058         c1 = cos(alpha);
00059         s1 = sin(alpha);
00060         if (!inverse)
00061             s1 = -s1;
00062         exptab[i].re = c1;
00063         exptab[i].im = s1;
00064     }
00065 }
00066 
00067 void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)
00068 {
00069     int n, i, j, k, n2;
00070     double tmp_re, tmp_im, s, c;
00071     FFTComplex *q;
00072 
00073     n = 1 << nbits;
00074     n2 = n >> 1;
00075     for(i=0;i<n;i++) {
00076         tmp_re = 0;
00077         tmp_im = 0;
00078         q = tab;
00079         for(j=0;j<n;j++) {
00080             k = (i * j) & (n - 1);
00081             if (k >= n2) {
00082                 c = -exptab[k - n2].re;
00083                 s = -exptab[k - n2].im;
00084             } else {
00085                 c = exptab[k].re;
00086                 s = exptab[k].im;
00087             }
00088             CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
00089             q++;
00090         }
00091         tabr[i].re = tmp_re;
00092         tabr[i].im = tmp_im;
00093     }
00094 }
00095 
00096 void imdct_ref(float *out, float *in, int nbits)
00097 {
00098     int n = 1<<nbits;
00099     int k, i, a;
00100     double sum, f;
00101 
00102     for(i=0;i<n;i++) {
00103         sum = 0;
00104         for(k=0;k<n/2;k++) {
00105             a = (2 * i + 1 + (n / 2)) * (2 * k + 1);
00106             f = cos(M_PI * a / (double)(2 * n));
00107             sum += f * in[k];
00108         }
00109         out[i] = -sum;
00110     }
00111 }
00112 
00113 /* NOTE: no normalisation by 1 / N is done */
00114 void mdct_ref(float *output, float *input, int nbits)
00115 {
00116     int n = 1<<nbits;
00117     int k, i;
00118     double a, s;
00119 
00120     /* do it by hand */
00121     for(k=0;k<n/2;k++) {
00122         s = 0;
00123         for(i=0;i<n;i++) {
00124             a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
00125             s += input[i] * cos(a);
00126         }
00127         output[k] = s;
00128     }
00129 }
00130 
00131 
00132 float frandom(void)
00133 {
00134     return (float)((random() & 0xffff) - 32768) / 32768.0;
00135 }
00136 
00137 int64_t gettime(void)
00138 {
00139     struct timeval tv;
00140     gettimeofday(&tv,NULL);
00141     return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
00142 }
00143 
00144 void check_diff(float *tab1, float *tab2, int n)
00145 {
00146     int i;
00147     double max= 0;
00148     double error= 0;
00149 
00150     for(i=0;i<n;i++) {
00151         double e= fabsf(tab1[i] - tab2[i]);
00152         if (e >= 1e-3) {
00153             av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n",
00154                    i, tab1[i], tab2[i]);
00155         }
00156         error+= e*e;
00157         if(e>max) max= e;
00158     }
00159     av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n);
00160 }
00161 
00162 
00163 void help(void)
00164 {
00165     av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
00166            "-h     print this help\n"
00167            "-s     speed test\n"
00168            "-m     (I)MDCT test\n"
00169            "-i     inverse transform test\n"
00170            "-n b   set the transform size to 2^b\n"
00171            );
00172     exit(1);
00173 }
00174 
00175 
00176 
00177 int main(int argc, char **argv)
00178 {
00179     FFTComplex *tab, *tab1, *tab_ref;
00180     FFTSample *tab2;
00181     int it, i, c;
00182     int do_speed = 0;
00183     int do_mdct = 0;
00184     int do_inverse = 0;
00185     FFTContext s1, *s = &s1;
00186     MDCTContext m1, *m = &m1;
00187     int fft_nbits, fft_size;
00188 
00189     fft_nbits = 9;
00190     for(;;) {
00191         c = getopt(argc, argv, "hsimn:");
00192         if (c == -1)
00193             break;
00194         switch(c) {
00195         case 'h':
00196             help();
00197             break;
00198         case 's':
00199             do_speed = 1;
00200             break;
00201         case 'i':
00202             do_inverse = 1;
00203             break;
00204         case 'm':
00205             do_mdct = 1;
00206             break;
00207         case 'n':
00208             fft_nbits = atoi(optarg);
00209             break;
00210         }
00211     }
00212 
00213     fft_size = 1 << fft_nbits;
00214     tab = av_malloc(fft_size * sizeof(FFTComplex));
00215     tab1 = av_malloc(fft_size * sizeof(FFTComplex));
00216     tab_ref = av_malloc(fft_size * sizeof(FFTComplex));
00217     tab2 = av_malloc(fft_size * sizeof(FFTSample));
00218 
00219     if (do_mdct) {
00220         if (do_inverse)
00221             av_log(NULL, AV_LOG_INFO,"IMDCT");
00222         else
00223             av_log(NULL, AV_LOG_INFO,"MDCT");
00224         ff_mdct_init(m, fft_nbits, do_inverse);
00225     } else {
00226         if (do_inverse)
00227             av_log(NULL, AV_LOG_INFO,"IFFT");
00228         else
00229             av_log(NULL, AV_LOG_INFO,"FFT");
00230         ff_fft_init(s, fft_nbits, do_inverse);
00231         fft_ref_init(fft_nbits, do_inverse);
00232     }
00233     av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
00234 
00235     /* generate random data */
00236 
00237     for(i=0;i<fft_size;i++) {
00238         tab1[i].re = frandom();
00239         tab1[i].im = frandom();
00240     }
00241 
00242     /* checking result */
00243     av_log(NULL, AV_LOG_INFO,"Checking...\n");
00244 
00245     if (do_mdct) {
00246         if (do_inverse) {
00247             imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits);
00248             ff_imdct_calc(m, tab2, (float *)tab1);
00249             check_diff((float *)tab_ref, tab2, fft_size);
00250         } else {
00251             mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits);
00252 
00253             ff_mdct_calc(m, tab2, (float *)tab1);
00254 
00255             check_diff((float *)tab_ref, tab2, fft_size / 2);
00256         }
00257     } else {
00258         memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
00259         ff_fft_permute(s, tab);
00260         ff_fft_calc(s, tab);
00261 
00262         fft_ref(tab_ref, tab1, fft_nbits);
00263         check_diff((float *)tab_ref, (float *)tab, fft_size * 2);
00264     }
00265 
00266     /* do a speed test */
00267 
00268     if (do_speed) {
00269         int64_t time_start, duration;
00270         int nb_its;
00271 
00272         av_log(NULL, AV_LOG_INFO,"Speed test...\n");
00273         /* we measure during about 1 seconds */
00274         nb_its = 1;
00275         for(;;) {
00276             time_start = gettime();
00277             for(it=0;it<nb_its;it++) {
00278                 if (do_mdct) {
00279                     if (do_inverse) {
00280                         ff_imdct_calc(m, (float *)tab, (float *)tab1);
00281                     } else {
00282                         ff_mdct_calc(m, (float *)tab, (float *)tab1);
00283                     }
00284                 } else {
00285                     memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
00286                     ff_fft_calc(s, tab);
00287                 }
00288             }
00289             duration = gettime() - time_start;
00290             if (duration >= 1000000)
00291                 break;
00292             nb_its *= 2;
00293         }
00294         av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
00295                (double)duration / nb_its,
00296                (double)duration / 1000000.0,
00297                nb_its);
00298     }
00299 
00300     if (do_mdct) {
00301         ff_mdct_end(m);
00302     } else {
00303         ff_fft_end(s);
00304     }
00305     return 0;
00306 }

Generated on Tue Nov 4 2014 12:59:21 for ffmpeg by  doxygen 1.7.1