00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "dsputil.h"
00022
00028
00029 #define BESSEL_I0_ITER 50 // default: 50 iterations of Bessel I0 approximation
00030 void ff_kbd_window_init(float *window, float alpha, int n)
00031 {
00032 int i, j;
00033 double sum = 0.0, bessel, tmp;
00034 double local_window[n];
00035 double alpha2 = (alpha * M_PI / n) * (alpha * M_PI / n);
00036
00037 for (i = 0; i < n; i++) {
00038 tmp = i * (n - i) * alpha2;
00039 bessel = 1.0;
00040 for (j = BESSEL_I0_ITER; j > 0; j--)
00041 bessel = bessel * tmp / (j * j) + 1;
00042 sum += bessel;
00043 local_window[i] = sum;
00044 }
00045
00046 sum++;
00047 for (i = 0; i < n; i++)
00048 window[i] = sqrt(local_window[i] / sum);
00049 }
00050
00054 int ff_mdct_init(MDCTContext *s, int nbits, int inverse)
00055 {
00056 int n, n4, i;
00057 float alpha;
00058
00059 memset(s, 0, sizeof(*s));
00060 n = 1 << nbits;
00061 s->nbits = nbits;
00062 s->n = n;
00063 n4 = n >> 2;
00064 s->tcos = av_malloc(n4 * sizeof(FFTSample));
00065 if (!s->tcos)
00066 goto fail;
00067 s->tsin = av_malloc(n4 * sizeof(FFTSample));
00068 if (!s->tsin)
00069 goto fail;
00070
00071 for(i=0;i<n4;i++) {
00072 alpha = 2 * M_PI * (i + 1.0 / 8.0) / n;
00073 s->tcos[i] = -cos(alpha);
00074 s->tsin[i] = -sin(alpha);
00075 }
00076 if (ff_fft_init(&s->fft, s->nbits - 2, inverse) < 0)
00077 goto fail;
00078 return 0;
00079 fail:
00080 av_freep(&s->tcos);
00081 av_freep(&s->tsin);
00082 return -1;
00083 }
00084
00085
00086 #define CMUL(pre, pim, are, aim, bre, bim) \
00087 {\
00088 float _are = (are);\
00089 float _aim = (aim);\
00090 float _bre = (bre);\
00091 float _bim = (bim);\
00092 (pre) = _are * _bre - _aim * _bim;\
00093 (pim) = _are * _bim + _aim * _bre;\
00094 }
00095
00102 void ff_imdct_calc(MDCTContext *s, FFTSample *output,
00103 const FFTSample *input, FFTSample *tmp)
00104 {
00105 int k, n8, n4, n2, n, j;
00106 const uint16_t *revtab = s->fft.revtab;
00107 const FFTSample *tcos = s->tcos;
00108 const FFTSample *tsin = s->tsin;
00109 const FFTSample *in1, *in2;
00110 FFTComplex *z = (FFTComplex *)tmp;
00111
00112 n = 1 << s->nbits;
00113 n2 = n >> 1;
00114 n4 = n >> 2;
00115 n8 = n >> 3;
00116
00117
00118 in1 = input;
00119 in2 = input + n2 - 1;
00120 for(k = 0; k < n4; k++) {
00121 j=revtab[k];
00122 CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
00123 in1 += 2;
00124 in2 -= 2;
00125 }
00126 ff_fft_calc(&s->fft, z);
00127
00128
00129
00130 for(k = 0; k < n4; k++) {
00131 CMUL(z[k].re, z[k].im, z[k].re, z[k].im, tcos[k], tsin[k]);
00132 }
00133 for(k = 0; k < n8; k++) {
00134 output[2*k] = -z[n8 + k].im;
00135 output[n2-1-2*k] = z[n8 + k].im;
00136
00137 output[2*k+1] = z[n8-1-k].re;
00138 output[n2-1-2*k-1] = -z[n8-1-k].re;
00139
00140 output[n2 + 2*k]=-z[k+n8].re;
00141 output[n-1- 2*k]=-z[k+n8].re;
00142
00143 output[n2 + 2*k+1]=z[n8-k-1].im;
00144 output[n-2 - 2 * k] = z[n8-k-1].im;
00145 }
00146 }
00147
00154 void ff_mdct_calc(MDCTContext *s, FFTSample *out,
00155 const FFTSample *input, FFTSample *tmp)
00156 {
00157 int i, j, n, n8, n4, n2, n3;
00158 FFTSample re, im, re1, im1;
00159 const uint16_t *revtab = s->fft.revtab;
00160 const FFTSample *tcos = s->tcos;
00161 const FFTSample *tsin = s->tsin;
00162 FFTComplex *x = (FFTComplex *)tmp;
00163
00164 n = 1 << s->nbits;
00165 n2 = n >> 1;
00166 n4 = n >> 2;
00167 n8 = n >> 3;
00168 n3 = 3 * n4;
00169
00170
00171 for(i=0;i<n8;i++) {
00172 re = -input[2*i+3*n4] - input[n3-1-2*i];
00173 im = -input[n4+2*i] + input[n4-1-2*i];
00174 j = revtab[i];
00175 CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
00176
00177 re = input[2*i] - input[n2-1-2*i];
00178 im = -(input[n2+2*i] + input[n-1-2*i]);
00179 j = revtab[n8 + i];
00180 CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
00181 }
00182
00183 ff_fft_calc(&s->fft, x);
00184
00185
00186 for(i=0;i<n4;i++) {
00187 re = x[i].re;
00188 im = x[i].im;
00189 CMUL(re1, im1, re, im, -tsin[i], -tcos[i]);
00190 out[2*i] = im1;
00191 out[n2-1-2*i] = re1;
00192 }
00193 }
00194
00195 void ff_mdct_end(MDCTContext *s)
00196 {
00197 av_freep(&s->tcos);
00198 av_freep(&s->tsin);
00199 ff_fft_end(&s->fft);
00200 }