|
|
|
@ -27,8 +27,10 @@ |
|
|
|
|
#include "libavutil/lfg.h" |
|
|
|
|
#include "libavutil/log.h" |
|
|
|
|
#include "fft.h" |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
#include "dct.h" |
|
|
|
|
#include "rdft.h" |
|
|
|
|
#endif |
|
|
|
|
#include <math.h> |
|
|
|
|
#include <unistd.h> |
|
|
|
|
#include <sys/time.h> |
|
|
|
@ -47,7 +49,19 @@ |
|
|
|
|
pim += (MUL16(are, bim) + MUL16(bre, aim));\
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
FFTComplex *exptab; |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
# define RANGE 1.0 |
|
|
|
|
# define REF_SCALE(x, bits) (x) |
|
|
|
|
# define FMT "%10.6f" |
|
|
|
|
#else |
|
|
|
|
# define RANGE 16384 |
|
|
|
|
# define REF_SCALE(x, bits) ((x) / (1<<(bits))) |
|
|
|
|
# define FMT "%6d" |
|
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
struct { |
|
|
|
|
float re, im; |
|
|
|
|
} *exptab; |
|
|
|
|
|
|
|
|
|
static void fft_ref_init(int nbits, int inverse) |
|
|
|
|
{ |
|
|
|
@ -55,7 +69,7 @@ static void fft_ref_init(int nbits, int inverse) |
|
|
|
|
double c1, s1, alpha; |
|
|
|
|
|
|
|
|
|
n = 1 << nbits; |
|
|
|
|
exptab = av_malloc((n / 2) * sizeof(FFTComplex)); |
|
|
|
|
exptab = av_malloc((n / 2) * sizeof(*exptab)); |
|
|
|
|
|
|
|
|
|
for (i = 0; i < (n/2); i++) { |
|
|
|
|
alpha = 2 * M_PI * (float)i / (float)n; |
|
|
|
@ -92,12 +106,12 @@ static void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) |
|
|
|
|
CMAC(tmp_re, tmp_im, c, s, q->re, q->im); |
|
|
|
|
q++; |
|
|
|
|
} |
|
|
|
|
tabr[i].re = tmp_re; |
|
|
|
|
tabr[i].im = tmp_im; |
|
|
|
|
tabr[i].re = REF_SCALE(tmp_re, nbits); |
|
|
|
|
tabr[i].im = REF_SCALE(tmp_im, nbits); |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
static void imdct_ref(float *out, float *in, int nbits) |
|
|
|
|
static void imdct_ref(FFTSample *out, FFTSample *in, int nbits) |
|
|
|
|
{ |
|
|
|
|
int n = 1<<nbits; |
|
|
|
|
int k, i, a; |
|
|
|
@ -110,12 +124,12 @@ static void imdct_ref(float *out, float *in, int nbits) |
|
|
|
|
f = cos(M_PI * a / (double)(2 * n)); |
|
|
|
|
sum += f * in[k]; |
|
|
|
|
} |
|
|
|
|
out[i] = -sum; |
|
|
|
|
out[i] = REF_SCALE(-sum, nbits - 2); |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
/* NOTE: no normalisation by 1 / N is done */ |
|
|
|
|
static void mdct_ref(float *output, float *input, int nbits) |
|
|
|
|
static void mdct_ref(FFTSample *output, FFTSample *input, int nbits) |
|
|
|
|
{ |
|
|
|
|
int n = 1<<nbits; |
|
|
|
|
int k, i; |
|
|
|
@ -128,10 +142,11 @@ static void mdct_ref(float *output, float *input, int nbits) |
|
|
|
|
a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); |
|
|
|
|
s += input[i] * cos(a); |
|
|
|
|
} |
|
|
|
|
output[k] = s; |
|
|
|
|
output[k] = REF_SCALE(s, nbits - 1); |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
static void idct_ref(float *output, float *input, int nbits) |
|
|
|
|
{ |
|
|
|
|
int n = 1<<nbits; |
|
|
|
@ -164,11 +179,12 @@ static void dct_ref(float *output, float *input, int nbits) |
|
|
|
|
output[k] = s; |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
#endif |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
static float frandom(AVLFG *prng) |
|
|
|
|
static FFTSample frandom(AVLFG *prng) |
|
|
|
|
{ |
|
|
|
|
return (int16_t)av_lfg_get(prng) / 32768.0; |
|
|
|
|
return (int16_t)av_lfg_get(prng) / 32768.0 * RANGE; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
static int64_t gettime(void) |
|
|
|
@ -178,7 +194,7 @@ static int64_t gettime(void) |
|
|
|
|
return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
static int check_diff(float *tab1, float *tab2, int n, double scale) |
|
|
|
|
static int check_diff(FFTSample *tab1, FFTSample *tab2, int n, double scale) |
|
|
|
|
{ |
|
|
|
|
int i; |
|
|
|
|
double max= 0; |
|
|
|
@ -186,9 +202,9 @@ static int check_diff(float *tab1, float *tab2, int n, double scale) |
|
|
|
|
int err = 0; |
|
|
|
|
|
|
|
|
|
for (i = 0; i < n; i++) { |
|
|
|
|
double e= fabsf(tab1[i] - (tab2[i] / scale)); |
|
|
|
|
double e = fabsf(tab1[i] - (tab2[i] / scale)) / RANGE; |
|
|
|
|
if (e >= 1e-3) { |
|
|
|
|
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: %10.6f %10.6f\n", |
|
|
|
|
av_log(NULL, AV_LOG_ERROR, "ERROR %5d: "FMT" "FMT"\n", |
|
|
|
|
i, tab1[i], tab2[i]); |
|
|
|
|
err = 1; |
|
|
|
|
} |
|
|
|
@ -233,8 +249,10 @@ int main(int argc, char **argv) |
|
|
|
|
int do_inverse = 0; |
|
|
|
|
FFTContext s1, *s = &s1; |
|
|
|
|
FFTContext m1, *m = &m1; |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
RDFTContext r1, *r = &r1; |
|
|
|
|
DCTContext d1, *d = &d1; |
|
|
|
|
#endif |
|
|
|
|
int fft_nbits, fft_size, fft_size_2; |
|
|
|
|
double scale = 1.0; |
|
|
|
|
AVLFG prng; |
|
|
|
@ -297,6 +315,7 @@ int main(int argc, char **argv) |
|
|
|
|
ff_fft_init(s, fft_nbits, do_inverse); |
|
|
|
|
fft_ref_init(fft_nbits, do_inverse); |
|
|
|
|
break; |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
case TRANSFORM_RDFT: |
|
|
|
|
if (do_inverse) |
|
|
|
|
av_log(NULL, AV_LOG_INFO,"IDFT_C2R"); |
|
|
|
@ -312,6 +331,10 @@ int main(int argc, char **argv) |
|
|
|
|
av_log(NULL, AV_LOG_INFO,"DCT_II"); |
|
|
|
|
ff_dct_init(d, fft_nbits, do_inverse ? DCT_III : DCT_II); |
|
|
|
|
break; |
|
|
|
|
#endif |
|
|
|
|
default: |
|
|
|
|
av_log(NULL, AV_LOG_ERROR, "Requested transform not supported\n"); |
|
|
|
|
return 1; |
|
|
|
|
} |
|
|
|
|
av_log(NULL, AV_LOG_INFO," %d test\n", fft_size); |
|
|
|
|
|
|
|
|
@ -328,15 +351,15 @@ int main(int argc, char **argv) |
|
|
|
|
switch (transform) { |
|
|
|
|
case TRANSFORM_MDCT: |
|
|
|
|
if (do_inverse) { |
|
|
|
|
imdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); |
|
|
|
|
m->imdct_calc(m, tab2, (float *)tab1); |
|
|
|
|
err = check_diff((float *)tab_ref, tab2, fft_size, scale); |
|
|
|
|
imdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits); |
|
|
|
|
m->imdct_calc(m, tab2, (FFTSample *)tab1); |
|
|
|
|
err = check_diff((FFTSample *)tab_ref, tab2, fft_size, scale); |
|
|
|
|
} else { |
|
|
|
|
mdct_ref((float *)tab_ref, (float *)tab1, fft_nbits); |
|
|
|
|
mdct_ref((FFTSample *)tab_ref, (FFTSample *)tab1, fft_nbits); |
|
|
|
|
|
|
|
|
|
m->mdct_calc(m, tab2, (float *)tab1); |
|
|
|
|
m->mdct_calc(m, tab2, (FFTSample *)tab1); |
|
|
|
|
|
|
|
|
|
err = check_diff((float *)tab_ref, tab2, fft_size / 2, scale); |
|
|
|
|
err = check_diff((FFTSample *)tab_ref, tab2, fft_size / 2, scale); |
|
|
|
|
} |
|
|
|
|
break; |
|
|
|
|
case TRANSFORM_FFT: |
|
|
|
@ -345,8 +368,9 @@ int main(int argc, char **argv) |
|
|
|
|
s->fft_calc(s, tab); |
|
|
|
|
|
|
|
|
|
fft_ref(tab_ref, tab1, fft_nbits); |
|
|
|
|
err = check_diff((float *)tab_ref, (float *)tab, fft_size * 2, 1.0); |
|
|
|
|
err = check_diff((FFTSample *)tab_ref, (FFTSample *)tab, fft_size * 2, 1.0); |
|
|
|
|
break; |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
case TRANSFORM_RDFT: |
|
|
|
|
if (do_inverse) { |
|
|
|
|
tab1[ 0].im = 0; |
|
|
|
@ -387,6 +411,7 @@ int main(int argc, char **argv) |
|
|
|
|
} |
|
|
|
|
err = check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0); |
|
|
|
|
break; |
|
|
|
|
#endif |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
/* do a speed test */ |
|
|
|
@ -404,15 +429,16 @@ int main(int argc, char **argv) |
|
|
|
|
switch (transform) { |
|
|
|
|
case TRANSFORM_MDCT: |
|
|
|
|
if (do_inverse) { |
|
|
|
|
m->imdct_calc(m, (float *)tab, (float *)tab1); |
|
|
|
|
m->imdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1); |
|
|
|
|
} else { |
|
|
|
|
m->mdct_calc(m, (float *)tab, (float *)tab1); |
|
|
|
|
m->mdct_calc(m, (FFTSample *)tab, (FFTSample *)tab1); |
|
|
|
|
} |
|
|
|
|
break; |
|
|
|
|
case TRANSFORM_FFT: |
|
|
|
|
memcpy(tab, tab1, fft_size * sizeof(FFTComplex)); |
|
|
|
|
s->fft_calc(s, tab); |
|
|
|
|
break; |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
case TRANSFORM_RDFT: |
|
|
|
|
memcpy(tab2, tab1, fft_size * sizeof(FFTSample)); |
|
|
|
|
r->rdft_calc(r, tab2); |
|
|
|
@ -421,6 +447,7 @@ int main(int argc, char **argv) |
|
|
|
|
memcpy(tab2, tab1, fft_size * sizeof(FFTSample)); |
|
|
|
|
d->dct_calc(d, tab2); |
|
|
|
|
break; |
|
|
|
|
#endif |
|
|
|
|
} |
|
|
|
|
} |
|
|
|
|
duration = gettime() - time_start; |
|
|
|
@ -441,12 +468,14 @@ int main(int argc, char **argv) |
|
|
|
|
case TRANSFORM_FFT: |
|
|
|
|
ff_fft_end(s); |
|
|
|
|
break; |
|
|
|
|
#if CONFIG_FFT_FLOAT |
|
|
|
|
case TRANSFORM_RDFT: |
|
|
|
|
ff_rdft_end(r); |
|
|
|
|
break; |
|
|
|
|
case TRANSFORM_DCT: |
|
|
|
|
ff_dct_end(d); |
|
|
|
|
break; |
|
|
|
|
#endif |
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
av_free(tab); |
|
|
|
|