27 |
27 |
#include "libavutil/lfg.h"
|
28 |
28 |
#include "libavutil/log.h"
|
29 |
29 |
#include "fft.h"
|
|
30 |
#if CONFIG_FFT_FLOAT
|
30 |
31 |
#include "dct.h"
|
31 |
32 |
#include "rdft.h"
|
|
33 |
#endif
|
32 |
34 |
#include <math.h>
|
33 |
35 |
#include <unistd.h>
|
34 |
36 |
#include <sys/time.h>
|
... | ... | |
47 |
49 |
pim += (MUL16(are, bim) + MUL16(bre, aim));\
|
48 |
50 |
}
|
49 |
51 |
|
50 |
|
FFTComplex *exptab;
|
|
52 |
#if CONFIG_FFT_FLOAT
|
|
53 |
# define RANGE 1.0
|
|
54 |
# define REF_SCALE(x, bits) (x)
|
|
55 |
# define FMT "%10.6f"
|
|
56 |
#else
|
|
57 |
# define RANGE 16384
|
|
58 |
# define REF_SCALE(x, bits) ((x) / (1<<(bits)))
|
|
59 |
# define FMT "%6d"
|
|
60 |
#endif
|
|
61 |
|
|
62 |
struct {
|
|
63 |
float re, im;
|
|
64 |
} *exptab;
|
51 |
65 |
|
52 |
66 |
static void fft_ref_init(int nbits, int inverse)
|
53 |
67 |
{
|
... | ... | |
55 |
69 |
double c1, s1, alpha;
|
56 |
70 |
|
57 |
71 |
n = 1 << nbits;
|
58 |
|
exptab = av_malloc((n / 2) * sizeof(FFTComplex));
|
|
72 |
exptab = av_malloc((n / 2) * sizeof(*exptab));
|
59 |
73 |
|
60 |
74 |
for (i = 0; i < (n/2); i++) {
|
61 |
75 |
alpha = 2 * M_PI * (float)i / (float)n;
|
... | ... | |
92 |
106 |
CMAC(tmp_re, tmp_im, c, s, q->re, q->im);
|
93 |
107 |
q++;
|
94 |
108 |
}
|
95 |
|
tabr[i].re = tmp_re;
|
96 |
|
tabr[i].im = tmp_im;
|
|
109 |
tabr[i].re = REF_SCALE(tmp_re, nbits);
|
|
110 |
tabr[i].im = REF_SCALE(tmp_im, nbits);
|
97 |
111 |
}
|
98 |
112 |
}
|
99 |
113 |
|
100 |
|
static void imdct_ref(float *out, float *in, int nbits)
|
|
114 |
static void imdct_ref(FFTSample *out, FFTSample *in, int nbits)
|
101 |
115 |
{
|
102 |
116 |
int n = 1<<nbits;
|
103 |
117 |
int k, i, a;
|
... | ... | |
110 |
124 |
f = cos(M_PI * a / (double)(2 * n));
|
111 |
125 |
sum += f * in[k];
|
112 |
126 |
}
|
113 |
|
out[i] = -sum;
|
|
127 |
out[i] = REF_SCALE(-sum, nbits - 2);
|
114 |
128 |
}
|
115 |
129 |
}
|
116 |
130 |
|
117 |
131 |
/* NOTE: no normalisation by 1 / N is done */
|
118 |
|
static void mdct_ref(float *output, float *input, int nbits)
|
|
132 |
static void mdct_ref(FFTSample *output, FFTSample *input, int nbits)
|
119 |
133 |
{
|
120 |
134 |
int n = 1<<nbits;
|
121 |
135 |
int k, i;
|
... | ... | |
128 |
142 |
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n));
|
129 |
143 |
s += input[i] * cos(a);
|
130 |
144 |
}
|
131 |
|
output[k] = s;
|
|
145 |
output[k] = REF_SCALE(s, nbits - 1);
|
132 |
146 |
}
|
133 |
147 |
}
|
134 |
148 |
|
|
149 |
#if CONFIG_FFT_FLOAT
|
135 |
150 |
static void idct_ref(float *output, float *input, int nbits)
|
136 |
151 |
{
|
137 |
152 |
int n = 1<<nbits;
|
... | ... | |
164 |
179 |
output[k] = s;
|
165 |
180 |
}
|
166 |
181 |
}
|
|
182 |
#endif
|
167 |
183 |
|
168 |
184 |
|
169 |
|
static float frandom(AVLFG *prng)
|
|
185 |
static FFTSample frandom(AVLFG *prng)
|
170 |
186 |
{
|
171 |
|
return (int16_t)av_lfg_get(prng) / 32768.0;
|
|
187 |
return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE;
|
172 |
188 |
}
|
173 |
189 |
|
174 |
190 |
static int64_t gettime(void)
|
... | ... | |
178 |
194 |
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
|
179 |
195 |
}
|
180 |
196 |
|
181 |
|
static int check_diff(float *tab1, float *tab2, int n, double scale)
|
|
197 |
static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale)
|
182 |
198 |
{
|
183 |
199 |
int i;
|
184 |
200 |
double max= 0;
|
... | ... | |
186 |
202 |
int err = 0;
|
187 |
203 |
|
188 |
204 |
for (i = 0; i < n; i++) {
|
189 |
|
double e= fabsf(tab1[i] - (tab2[i] / scale));
|
|
205 |
double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE;
|
190 |
206 |
if (e >= 1e-3) {
|
191 |
|
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: %10.6f %10.6f\n",
|
|
207 |
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n",
|
192 |
208 |
i, tab1[i], tab2[i]);
|
193 |
209 |
err = 1;
|
194 |
210 |
}
|
... | ... | |
233 |
249 |
int do_inverse = 0;
|
234 |
250 |
FFTContext s1, *s = &s1;
|
235 |
251 |
FFTContext m1, *m = &m1;
|
|
252 |
#if CONFIG_FFT_FLOAT
|
236 |
253 |
RDFTContext r1, *r = &r1;
|
237 |
254 |
DCTContext d1, *d = &d1;
|
|
255 |
#endif
|
238 |
256 |
int fft_nbits, fft_size, fft_size_2;
|
239 |
257 |
double scale = 1.0;
|
240 |
258 |
AVLFG prng;
|
... | ... | |
297 |
315 |
ff_fft_init(s, fft_nbits, do_inverse);
|
298 |
316 |
fft_ref_init(fft_nbits, do_inverse);
|
299 |
317 |
break;
|
|
318 |
#if CONFIG_FFT_FLOAT
|
300 |
319 |
case TRANSFORM_RDFT:
|
301 |
320 |
if (do_inverse)
|
302 |
321 |
av_log(NULL, AV_LOG_INFO,"IDFT_C2R");
|
... | ... | |
312 |
331 |
av_log(NULL, AV_LOG_INFO,"DCT_II");
|
313 |
332 |
ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II);
|
314 |
333 |
break;
|
|
334 |
#endif
|
|
335 |
default:
|
|
336 |
av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n");
|
|
337 |
return 1;
|
315 |
338 |
}
|
316 |
339 |
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);
|
317 |
340 |
|
... | ... | |
328 |
351 |
switch (transform) {
|
329 |
352 |
case TRANSFORM_MDCT:
|
330 |
353 |
if (do_inverse) {
|
331 |
|
imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits);
|
332 |
|
m->imdct_calc(m, tab2, (float *)tab1);
|
333 |
|
err = check_diff((float *)tab_ref, tab2, fft_size, scale);
|
|
354 |
imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
|
|
355 |
m->imdct_calc(m, tab2, (FFTSample *)tab1);
|
|
356 |
err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale);
|
334 |
357 |
} else {
|
335 |
|
mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits);
|
|
358 |
mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits);
|
336 |
359 |
|
337 |
|
m->mdct_calc(m, tab2, (float *)tab1);
|
|
360 |
m->mdct_calc(m, tab2, (FFTSample *)tab1);
|
338 |
361 |
|
339 |
|
err = check_diff((float *)tab_ref, tab2, fft_size / 2, scale);
|
|
362 |
err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale);
|
340 |
363 |
}
|
341 |
364 |
break;
|
342 |
365 |
case TRANSFORM_FFT:
|
... | ... | |
345 |
368 |
s->fft_calc(s, tab);
|
346 |
369 |
|
347 |
370 |
fft_ref(tab_ref, tab1, fft_nbits);
|
348 |
|
err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 1.0);
|
|
371 |
err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0);
|
349 |
372 |
break;
|
|
373 |
#if CONFIG_FFT_FLOAT
|
350 |
374 |
case TRANSFORM_RDFT:
|
351 |
375 |
if (do_inverse) {
|
352 |
376 |
tab1[ 0].im = 0;
|
... | ... | |
387 |
411 |
}
|
388 |
412 |
err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);
|
389 |
413 |
break;
|
|
414 |
#endif
|
390 |
415 |
}
|
391 |
416 |
|
392 |
417 |
/* do a speed test */
|
... | ... | |
404 |
429 |
switch (transform) {
|
405 |
430 |
case TRANSFORM_MDCT:
|
406 |
431 |
if (do_inverse) {
|
407 |
|
m->imdct_calc(m, (float *)tab, (float *)tab1);
|
|
432 |
m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
|
408 |
433 |
} else {
|
409 |
|
m->mdct_calc(m, (float *)tab, (float *)tab1);
|
|
434 |
m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1);
|
410 |
435 |
}
|
411 |
436 |
break;
|
412 |
437 |
case TRANSFORM_FFT:
|
413 |
438 |
memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
|
414 |
439 |
s->fft_calc(s, tab);
|
415 |
440 |
break;
|
|
441 |
#if CONFIG_FFT_FLOAT
|
416 |
442 |
case TRANSFORM_RDFT:
|
417 |
443 |
memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
|
418 |
444 |
r->rdft_calc(r, tab2);
|
... | ... | |
421 |
447 |
memcpy(tab2, tab1, fft_size * sizeof(FFTSample));
|
422 |
448 |
d->dct_calc(d, tab2);
|
423 |
449 |
break;
|
|
450 |
#endif
|
424 |
451 |
}
|
425 |
452 |
}
|
426 |
453 |
duration = gettime() - time_start;
|
... | ... | |
441 |
468 |
case TRANSFORM_FFT:
|
442 |
469 |
ff_fft_end(s);
|
443 |
470 |
break;
|
|
471 |
#if CONFIG_FFT_FLOAT
|
444 |
472 |
case TRANSFORM_RDFT:
|
445 |
473 |
ff_rdft_end(r);
|
446 |
474 |
break;
|
447 |
475 |
case TRANSFORM_DCT:
|
448 |
476 |
ff_dct_end(d);
|
449 |
477 |
break;
|
|
478 |
#endif
|
450 |
479 |
}
|
451 |
480 |
|
452 |
481 |
av_free(tab);
|