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 <stdlib.h>

41 #include <string.h>

43 /* reference fft */

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

47 #define CMAC(pre, pim, are, aim, bre, bim) \

48 {\

49 pre += (MUL16(are, bre) - MUL16(aim, bim));\

50 pim += (MUL16(are, bim) + MUL16(bre, aim));\

51 }

53 #if CONFIG_FFT_FLOAT

54 # define RANGE 1.0

55 # define REF_SCALE(x, bits) (x)

56 # define FMT "%10.6f"

57 #else

58 # define RANGE 16384

59 # define REF_SCALE(x, bits) ((x) / (1<<(bits)))

60 # define FMT "%6d"

61 #endif

63 struct {

64 float re, im;

65 } *exptab;

67 static void fft_ref_init(int nbits, int inverse)

68 {

69 int n, i;

70 double c1, s1, alpha;

72 n = 1 << nbits;

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

75 for (i = 0; i < (n/2); i++) {

76 alpha = 2 * M_PI * (float)i / (float)n;

77 c1 = cos(alpha);

78 s1 = sin(alpha);

79 if (!inverse)

80 s1 = -s1;

81 exptab[i].re = c1;

82 exptab[i].im = s1;

83 }

84 }

86 static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits)

87 {

88 int n, i, j, k, n2;

89 double tmp_re, tmp_im, s, c;

90 FFTComplex *q;

92 n = 1 << nbits;

93 n2 = n >> 1;

94 for (i = 0; i < n; i++) {

95 tmp_re = 0;

96 tmp_im = 0;

97 q = tab;

98 for (j = 0; j < n; j++) {

99 k = (i * j) & (n - 1);

100 if (k >= n2) {

101 c = -exptab[k - n2].re;

102 s = -exptab[k - n2].im;

103 } else {

104 c = exptab[k].re;

105 s = exptab[k].im;

106 }

107 CMAC(tmp_re, tmp_im, c, s, q->re, q->im);

108 q++;

109 }

110 tabr[i].re = REF_SCALE(tmp_re, nbits);

111 tabr[i].im = REF_SCALE(tmp_im, nbits);

112 }

113 }

114

115 static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)

116 {

117 int n = 1<<nbits;

118 int k, i, a;

119 double sum, f;

120

121 for (i = 0; i < n; i++) {

122 sum = 0;

123 for (k = 0; k < n/2; k++) {

124 a = (2 * i + 1 + (n / 2)) * (2 * k + 1);

125 f = cos(M_PI * a / (double)(2 * n));

126 sum += f * in[k];

127 }

128 out[i] = REF_SCALE(-sum, nbits - 2);

129 }

130 }

131

132 /* NOTE: no normalisation by 1 / N is done */

133 static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)

134 {

135 int n = 1<<nbits;

136 int k, i;

137 double a, s;

138

139 /* do it by hand */

140 for (k = 0; k < n/2; k++) {

141 s = 0;

142 for (i = 0; i < n; i++) {

143 a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));

144 s += input[i] * cos(a);

145 }

146 output[k] = REF_SCALE(s, nbits - 1);

147 }

148 }

149

150 #if CONFIG_FFT_FLOAT

151 static void idct_ref(FFTSample *output, FFTSample *input, int nbits)

152 {

153 int n = 1<<nbits;

154 int k, i;

155 double a, s;

156

157 /* do it by hand */

158 for (i = 0; i < n; i++) {

159 s = 0.5 * input[0];

160 for (k = 1; k < n; k++) {

161 a = M_PI*k*(i+0.5) / n;

162 s += input[k] * cos(a);

163 }

164 output[i] = 2 * s / n;

165 }

166 }

167 static void dct_ref(FFTSample *output, FFTSample *input, int nbits)

168 {

169 int n = 1<<nbits;

170 int k, i;

171 double a, s;

172

173 /* do it by hand */

174 for (k = 0; k < n; k++) {

175 s = 0;

176 for (i = 0; i < n; i++) {

177 a = M_PI*k*(i+0.5) / n;

178 s += input[i] * cos(a);

179 }

180 output[k] = s;

181 }

182 }

183 #endif

184

185

186 static FFTSample frandom(AVLFG *prng)

187 {

188 return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;

189 }

190

191 static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)

192 {

193 int i;

194 double max= 0;

195 double error= 0;

196 int err = 0;

197

198 for (i = 0; i < n; i++) {

199 double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;

200 if (e >= 1e-3) {

201 av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT " "FMT "\n",

202 i, tab1[i], tab2[i]);

203 err = 1;

204 }

205 error+= e*e;

206 if(e>max) max= e;

207 }

208 av_log(NULL, AV_LOG_INFO, "max:%f e:%g\n", max, sqrt(error/n));

209 return err;

210 }

211

212

213 static void help(void)

214 {

215 av_log(NULL, AV_LOG_INFO,"usage: fft-test [-h] [-s] [-i] [-n b]\n"

216 "-h print this help\n"

217 "-s speed test\n"

218 "-m (I)MDCT test\n"

219 "-d (I)DCT test\n"

220 "-r (I)RDFT test\n"

221 "-i inverse transform test\n"

222 "-n b set the transform size to 2^b\n"

223 "-f x set scale factor for output data of (I)MDCT to x\n"

224 );

225 }

226

227 enum tf_transform {

228 TRANSFORM_FFT,

229 TRANSFORM_MDCT,

230 TRANSFORM_RDFT,

231 TRANSFORM_DCT,

232 };

233

234 #if !HAVE_GETOPT

235 #include "compat/getopt.c"

236 #endif

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 cpuflags;

244 int do_speed = 0;

245 int err = 1;

246 enum tf_transform transform = TRANSFORM_FFT;

247 int do_inverse = 0;

248 FFTContext s1, *s = &s1;

249 FFTContext m1, *m = &m1;

250 #if CONFIG_FFT_FLOAT

251 RDFTContext r1, *r = &r1;

252 DCTContext d1, *d = &d1;

253 int fft_size_2;

254 #endif

255 int fft_nbits, fft_size;

256 double scale = 1.0;

257 AVLFG prng;

258 av_lfg_init(&prng, 1);

259

260 fft_nbits = 9;

261 for(;;) {

262 c = getopt(argc, argv, "hsimrdn:f:c:");

263 if (c == -1)

264 break;

265 switch(c) {

266 case 'h':

267 help();

268 return 1;

269 case 's':

270 do_speed = 1;

271 break;

272 case 'i':

273 do_inverse = 1;

274 break;

275 case 'm':

276 transform = TRANSFORM_MDCT;

277 break;

278 case 'r':

279 transform = TRANSFORM_RDFT;

280 break;

281 case 'd':

282 transform = TRANSFORM_DCT;

283 break;

284 case 'n':

285 fft_nbits = atoi(optarg);

286 break;

287 case 'f':

288 scale = atof(optarg);

289 break;

290 case 'c':

291 cpuflags = av_get_cpu_flags();

292

293 if (av_parse_cpu_caps(&cpuflags, optarg) < 0)

294 return 1;

295

296 av_force_cpu_flags(cpuflags);

297 break;

298 }

299 }

300

301 fft_size = 1 << fft_nbits;

302 tab = av_malloc(fft_size * sizeof(FFTComplex));

303 tab1 = av_malloc(fft_size * sizeof(FFTComplex));

304 tab_ref = av_malloc(fft_size * sizeof(FFTComplex));

305 tab2 = av_malloc(fft_size * sizeof(FFTSample));

306

307 switch (transform) {

308 case TRANSFORM_MDCT:

309 av_log(NULL, AV_LOG_INFO,"Scale factor is set to %f\n", scale);

310 if (do_inverse)

311 av_log(NULL, AV_LOG_INFO,"IMDCT");

312 else

313 av_log(NULL, AV_LOG_INFO,"MDCT");

314 ff_mdct_init(m, fft_nbits, do_inverse, scale);

315 break;

316 case TRANSFORM_FFT:

317 if (do_inverse)

318 av_log(NULL, AV_LOG_INFO,"IFFT");

319 else

320 av_log(NULL, AV_LOG_INFO,"FFT");

321 ff_fft_init(s, fft_nbits, do_inverse);

322 fft_ref_init(fft_nbits, do_inverse);

323 break;

324 #if CONFIG_FFT_FLOAT

325 case TRANSFORM_RDFT:

326 if (do_inverse)

327 av_log(NULL, AV_LOG_INFO,"IDFT_C2R");

328 else

329 av_log(NULL, AV_LOG_INFO,"DFT_R2C");

330 ff_rdft_init(r, fft_nbits, do_inverse ? IDFT_C2R : DFT_R2C);

331 fft_ref_init(fft_nbits, do_inverse);

332 break;

333 case TRANSFORM_DCT:

334 if (do_inverse)

335 av_log(NULL, AV_LOG_INFO,"DCT_III");

336 else

337 av_log(NULL, AV_LOG_INFO,"DCT_II");

338 ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);

339 break;

340 #endif

341 default:

342 av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");

343 return 1;

344 }

345 av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);

346

347 /* generate random data */

348

349 for (i = 0; i < fft_size; i++) {

350 tab1[i].re = frandom(&prng);

351 tab1[i].im = frandom(&prng);

352 }

353

354 /* checking result */

355 av_log(NULL, AV_LOG_INFO,"Checking...\n");

356

357 switch (transform) {

358 case TRANSFORM_MDCT:

359 if (do_inverse) {

360 imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);

361 m->imdct_calc(m, tab2, (FFTSample *)tab1);

362 err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);

363 } else {

364 mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);

365

366 m->mdct_calc(m, tab2, (FFTSample *)tab1);

367

368 err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);

369 }

370 break;

371 case TRANSFORM_FFT:

372 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

373 s->fft_permute(s, tab);

374 s->fft_calc(s, tab);

375

376 fft_ref(tab_ref, tab1, fft_nbits);

377 err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);

378 break;

379 #if CONFIG_FFT_FLOAT

380 case TRANSFORM_RDFT:

381 fft_size_2 = fft_size >> 1;

382 if (do_inverse) {

383 tab1[ 0].im = 0;

384 tab1[fft_size_2].im = 0;

385 for (i = 1; i < fft_size_2; i++) {

386 tab1[fft_size_2+i].re = tab1[fft_size_2-i].re;

387 tab1[fft_size_2+i].im = -tab1[fft_size_2-i].im;

388 }

389

390 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

391 tab2[1] = tab1[fft_size_2].re;

392

393 r->rdft_calc(r, tab2);

394 fft_ref(tab_ref, tab1, fft_nbits);

395 for (i = 0; i < fft_size; i++) {

396 tab[i].re = tab2[i];

397 tab[i].im = 0;

398 }

399 err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 0.5);

400 } else {

401 for (i = 0; i < fft_size; i++) {

402 tab2[i] = tab1[i].re;

403 tab1[i].im = 0;

404 }

405 r->rdft_calc(r, tab2);

406 fft_ref(tab_ref, tab1, fft_nbits);

407 tab_ref[0].im = tab_ref[fft_size_2].re;

408 err = check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);

409 }

410 break;

411 case TRANSFORM_DCT:

412 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

413 d->dct_calc(d, (FFTSample *)tab);

414 if (do_inverse) {

415 idct_ref((FFTSample*)tab_ref, (FFTSample *)tab1, fft_nbits);

416 } else {

417 dct_ref((FFTSample*)tab_ref, (FFTSample *)tab1, fft_nbits);

418 }

419 err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);

420 break;

421 #endif

422 }

423

424 /* do a speed test */

425

426 if (do_speed) {

427 int64_t time_start, duration;

428 int nb_its;

429

430 av_log(NULL, AV_LOG_INFO,"Speed test...\n");

431 /* we measure during about 1 seconds */

432 nb_its = 1;

433 for(;;) {

434 time_start = av_gettime();

435 for (it = 0; it < nb_its; it++) {

436 switch (transform) {

437 case TRANSFORM_MDCT:

438 if (do_inverse) {

439 m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);

440 } else {

441 m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);

442 }

443 break;

444 case TRANSFORM_FFT:

445 memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

446 s->fft_calc(s, tab);

447 break;

448 #if CONFIG_FFT_FLOAT

449 case TRANSFORM_RDFT:

450 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

451 r->rdft_calc(r, tab2);

452 break;

453 case TRANSFORM_DCT:

454 memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

455 d->dct_calc(d, tab2);

456 break;

457 #endif

458 }

459 }

460 duration = av_gettime() - time_start;

461 if (duration >= 1000000)

462 break;

463 nb_its *= 2;

464 }

465 av_log(NULL, AV_LOG_INFO,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",

466 (double)duration / nb_its,

467 (double)duration / 1000000.0,

468 nb_its);

469 }

470

471 switch (transform) {

472 case TRANSFORM_MDCT:

473 ff_mdct_end(m);

474 break;

475 case TRANSFORM_FFT:

476 ff_fft_end(s);

477 break;

478 #if CONFIG_FFT_FLOAT

479 case TRANSFORM_RDFT:

480 ff_rdft_end(r);

481 break;

482 case TRANSFORM_DCT:

483 ff_dct_end(d);

484 break;

485 #endif

486 }

487

488 av_free(tab);

489 av_free(tab1);

490 av_free(tab2);

491 av_free(tab_ref);

492 av_free(exptab);

493

494 return err;

495 }

Generated on Wed Jul 10 2013 23:47:57 for FFmpeg by doxygen 1.8.2