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 libavcodec/fft-test.c 23 * FFT and MDCT tests. 24 */ 25 26#include "dsputil.h" 27#include <math.h> 28#include <unistd.h> 29#include <sys/time.h> 30#include <stdlib.h> 31#include <string.h> 32 33#undef exit 34#undef random 35 36/* reference fft */ 37 38#define MUL16(a,b) ((a) * (b)) 39 40#define CMAC(pre, pim, are, aim, bre, bim) \ 41{\ 42 pre += (MUL16(are, bre) - MUL16(aim, bim));\ 43 pim += (MUL16(are, bim) + MUL16(bre, aim));\ 44} 45 46FFTComplex *exptab; 47 48void fft_ref_init(int nbits, int inverse) 49{ 50 int n, i; 51 double c1, s1, alpha; 52 53 n = 1 << nbits; 54 exptab = av_malloc((n / 2) * sizeof(FFTComplex)); 55 56 for(i=0;i<(n/2);i++) { 57 alpha = 2 * M_PI * (float)i / (float)n; 58 c1 = cos(alpha); 59 s1 = sin(alpha); 60 if (!inverse) 61 s1 = -s1; 62 exptab[i].re = c1; 63 exptab[i].im = s1; 64 } 65} 66 67void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) 68{ 69 int n, i, j, k, n2; 70 double tmp_re, tmp_im, s, c; 71 FFTComplex *q; 72 73 n = 1 << nbits; 74 n2 = n >> 1; 75 for(i=0;i<n;i++) { 76 tmp_re = 0; 77 tmp_im = 0; 78 q = tab; 79 for(j=0;j<n;j++) { 80 k = (i * j) & (n - 1); 81 if (k >= n2) { 82 c = -exptab[k - n2].re; 83 s = -exptab[k - n2].im; 84 } else { 85 c = exptab[k].re; 86 s = exptab[k].im; 87 } 88 CMAC(tmp_re, tmp_im, c, s, q->re, q->im); 89 q++; 90 } 91 tabr[i].re = tmp_re; 92 tabr[i].im = tmp_im; 93 } 94} 95 96void imdct_ref(float *out, float *in, int nbits) 97{ 98 int n = 1<<nbits; 99 int k, i, a; 100 double sum, f; 101 102 for(i=0;i<n;i++) { 103 sum = 0; 104 for(k=0;k<n/2;k++) { 105 a = (2 * i + 1 + (n / 2)) * (2 * k + 1); 106 f = cos(M_PI * a / (double)(2 * n)); 107 sum += f * in[k]; 108 } 109 out[i] = -sum; 110 } 111} 112 113/* NOTE: no normalisation by 1 / N is done */ 114void mdct_ref(float *output, float *input, int nbits) 115{ 116 int n = 1<<nbits; 117 int k, i; 118 double a, s; 119 120 /* do it by hand */ 121 for(k=0;k<n/2;k++) { 122 s = 0; 123 for(i=0;i<n;i++) { 124 a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); 125 s += input[i] * cos(a); 126 } 127 output[k] = s; 128 } 129} 130 131 132float frandom(void) 133{ 134 return (float)((random() & 0xffff) - 32768) / 32768.0; 135} 136 137int64_t gettime(void) 138{ 139 struct timeval tv; 140 gettimeofday(&tv,NULL); 141 return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; 142} 143 144void check_diff(float *tab1, float *tab2, int n) 145{ 146 int i; 147 double max= 0; 148 double error= 0; 149 150 for(i=0;i<n;i++) { 151 double e= fabsf(tab1[i] - tab2[i]); 152 if (e >= 1e-3) { 153 av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n", 154 i, tab1[i], tab2[i]); 155 } 156 error+= e*e; 157 if(e>max) max= e; 158 } 159 av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error)/n); 160} 161 162 163void help(void) 164{ 165 av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n" 166 "-h print this help\n" 167 "-s speed test\n" 168 "-m (I)MDCT test\n" 169 "-i inverse transform test\n" 170 "-n b set the transform size to 2^b\n" 171 ); 172 exit(1); 173} 174 175 176 177int main(int argc, char **argv) 178{ 179 FFTComplex *tab, *tab1, *tab_ref; 180 FFTSample *tab2; 181 int it, i, c; 182 int do_speed = 0; 183 int do_mdct = 0; 184 int do_inverse = 0; 185 FFTContext s1, *s = &s1; 186 MDCTContext m1, *m = &m1; 187 int fft_nbits, fft_size; 188 189 fft_nbits = 9; 190 for(;;) { 191 c = getopt(argc, argv, "hsimn:"); 192 if (c == -1) 193 break; 194 switch(c) { 195 case 'h': 196 help(); 197 break; 198 case 's': 199 do_speed = 1; 200 break; 201 case 'i': 202 do_inverse = 1; 203 break; 204 case 'm': 205 do_mdct = 1; 206 break; 207 case 'n': 208 fft_nbits = atoi(optarg); 209 break; 210 } 211 } 212 213 fft_size = 1 << fft_nbits; 214 tab = av_malloc(fft_size * sizeof(FFTComplex)); 215 tab1 = av_malloc(fft_size * sizeof(FFTComplex)); 216 tab_ref = av_malloc(fft_size * sizeof(FFTComplex)); 217 tab2 = av_malloc(fft_size * sizeof(FFTSample)); 218 219 if (do_mdct) { 220 if (do_inverse) 221 av_log(NULL, AV_LOG_INFO,"IMDCT"); 222 else 223 av_log(NULL, AV_LOG_INFO,"MDCT"); 224 ff_mdct_init(m, fft_nbits, do_inverse); 225 } else { 226 if (do_inverse) 227 av_log(NULL, AV_LOG_INFO,"IFFT"); 228 else 229 av_log(NULL, AV_LOG_INFO,"FFT"); 230 ff_fft_init(s, fft_nbits, do_inverse); 231 fft_ref_init(fft_nbits, do_inverse); 232 } 233 av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); 234 235 /* generate random data */ 236 237 for(i=0;i<fft_size;i++) { 238 tab1[i].re = frandom(); 239 tab1[i].im = frandom(); 240 } 241 242 /* checking result */ 243 av_log(NULL, AV_LOG_INFO,"Checking...\n"); 244 245 if (do_mdct) { 246 if (do_inverse) { 247 imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); 248 ff_imdct_calc(m, tab2, (float *)tab1); 249 check_diff((float *)tab_ref, tab2, fft_size); 250 } else { 251 mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); 252 253 ff_mdct_calc(m, tab2, (float *)tab1); 254 255 check_diff((float *)tab_ref, tab2, fft_size / 2); 256 } 257 } else { 258 memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); 259 ff_fft_permute(s, tab); 260 ff_fft_calc(s, tab); 261 262 fft_ref(tab_ref, tab1, fft_nbits); 263 check_diff((float *)tab_ref, (float *)tab, fft_size * 2); 264 } 265 266 /* do a speed test */ 267 268 if (do_speed) { 269 int64_t time_start, duration; 270 int nb_its; 271 272 av_log(NULL, AV_LOG_INFO,"Speed test...\n"); 273 /* we measure during about 1 seconds */ 274 nb_its = 1; 275 for(;;) { 276 time_start = gettime(); 277 for(it=0;it<nb_its;it++) { 278 if (do_mdct) { 279 if (do_inverse) { 280 ff_imdct_calc(m, (float *)tab, (float *)tab1); 281 } else { 282 ff_mdct_calc(m, (float *)tab, (float *)tab1); 283 } 284 } else { 285 memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); 286 ff_fft_calc(s, tab); 287 } 288 } 289 duration = gettime() - time_start; 290 if (duration >= 1000000) 291 break; 292 nb_its *= 2; 293 } 294 av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n", 295 (double)duration / nb_its, 296 (double)duration / 1000000.0, 297 nb_its); 298 } 299 300 if (do_mdct) { 301 ff_mdct_end(m); 302 } else { 303 ff_fft_end(s); 304 } 305 return 0; 306} 307