FFmpeg: libavcodec/fft-test.c Source File

FFmpeg

[フレーム]

libavcodec

fft-test.c

Go to the documentation of this file.

1 /*

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 */

21 /**

22 * @file

23 * FFT and MDCT tests.

24 */

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>

44 /* reference fft */

46 #define MUL16(a,b) ((a) * (b))

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 }

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

68 struct {

69 float re, im;

70 } *exptab;

72 static void fft_ref_init(int nbits, int inverse)

73 {

74 int n, i;

75 double c1, s1, alpha;

77 n = 1 << nbits;

78 exptab = av_malloc((n / 2) * sizeof(*exptab));

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 }

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;

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 }

Generated on Sat Jan 25 2014 19:51:47 for FFmpeg by doxygen 1.8.2