00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00030 #include <math.h>
00031 #include "libavutil/mathematics.h"
00032 #include "fft.h"
00033 #include "x86/fft.h"
00034
00035 #define DCT32_FLOAT
00036 #include "dct32.c"
00037
00038
00039 #define SIN(s,n,x) (s->costab[(n) - (x)])
00040
00041
00042 #define COS(s,n,x) (s->costab[x])
00043
00044 static void ff_dst_calc_I_c(DCTContext *ctx, FFTSample *data)
00045 {
00046 int n = 1 << ctx->nbits;
00047 int i;
00048
00049 data[0] = 0;
00050 for(i = 1; i < n/2; i++) {
00051 float tmp1 = data[i ];
00052 float tmp2 = data[n - i];
00053 float s = SIN(ctx, n, 2*i);
00054
00055 s *= tmp1 + tmp2;
00056 tmp1 = (tmp1 - tmp2) * 0.5f;
00057 data[i ] = s + tmp1;
00058 data[n - i] = s - tmp1;
00059 }
00060
00061 data[n/2] *= 2;
00062 ff_rdft_calc(&ctx->rdft, data);
00063
00064 data[0] *= 0.5f;
00065
00066 for(i = 1; i < n-2; i += 2) {
00067 data[i + 1] += data[i - 1];
00068 data[i ] = -data[i + 2];
00069 }
00070
00071 data[n-1] = 0;
00072 }
00073
00074 static void ff_dct_calc_I_c(DCTContext *ctx, FFTSample *data)
00075 {
00076 int n = 1 << ctx->nbits;
00077 int i;
00078 float next = -0.5f * (data[0] - data[n]);
00079
00080 for(i = 0; i < n/2; i++) {
00081 float tmp1 = data[i ];
00082 float tmp2 = data[n - i];
00083 float s = SIN(ctx, n, 2*i);
00084 float c = COS(ctx, n, 2*i);
00085
00086 c *= tmp1 - tmp2;
00087 s *= tmp1 - tmp2;
00088
00089 next += c;
00090
00091 tmp1 = (tmp1 + tmp2) * 0.5f;
00092 data[i ] = tmp1 - s;
00093 data[n - i] = tmp1 + s;
00094 }
00095
00096 ff_rdft_calc(&ctx->rdft, data);
00097 data[n] = data[1];
00098 data[1] = next;
00099
00100 for(i = 3; i <= n; i += 2)
00101 data[i] = data[i - 2] - data[i];
00102 }
00103
00104 static void ff_dct_calc_III_c(DCTContext *ctx, FFTSample *data)
00105 {
00106 int n = 1 << ctx->nbits;
00107 int i;
00108
00109 float next = data[n - 1];
00110 float inv_n = 1.0f / n;
00111
00112 for (i = n - 2; i >= 2; i -= 2) {
00113 float val1 = data[i ];
00114 float val2 = data[i - 1] - data[i + 1];
00115 float c = COS(ctx, n, i);
00116 float s = SIN(ctx, n, i);
00117
00118 data[i ] = c * val1 + s * val2;
00119 data[i + 1] = s * val1 - c * val2;
00120 }
00121
00122 data[1] = 2 * next;
00123
00124 ff_rdft_calc(&ctx->rdft, data);
00125
00126 for (i = 0; i < n / 2; i++) {
00127 float tmp1 = data[i ] * inv_n;
00128 float tmp2 = data[n - i - 1] * inv_n;
00129 float csc = ctx->csc2[i] * (tmp1 - tmp2);
00130
00131 tmp1 += tmp2;
00132 data[i ] = tmp1 + csc;
00133 data[n - i - 1] = tmp1 - csc;
00134 }
00135 }
00136
00137 static void ff_dct_calc_II_c(DCTContext *ctx, FFTSample *data)
00138 {
00139 int n = 1 << ctx->nbits;
00140 int i;
00141 float next;
00142
00143 for (i=0; i < n/2; i++) {
00144 float tmp1 = data[i ];
00145 float tmp2 = data[n - i - 1];
00146 float s = SIN(ctx, n, 2*i + 1);
00147
00148 s *= tmp1 - tmp2;
00149 tmp1 = (tmp1 + tmp2) * 0.5f;
00150
00151 data[i ] = tmp1 + s;
00152 data[n-i-1] = tmp1 - s;
00153 }
00154
00155 ff_rdft_calc(&ctx->rdft, data);
00156
00157 next = data[1] * 0.5;
00158 data[1] *= -1;
00159
00160 for (i = n - 2; i >= 0; i -= 2) {
00161 float inr = data[i ];
00162 float ini = data[i + 1];
00163 float c = COS(ctx, n, i);
00164 float s = SIN(ctx, n, i);
00165
00166 data[i ] = c * inr + s * ini;
00167
00168 data[i+1] = next;
00169
00170 next += s * inr - c * ini;
00171 }
00172 }
00173
00174 static void dct32_func(DCTContext *ctx, FFTSample *data)
00175 {
00176 ctx->dct32(data, data);
00177 }
00178
00179 void ff_dct_calc(DCTContext *s, FFTSample *data)
00180 {
00181 s->dct_calc(s, data);
00182 }
00183
00184 av_cold int ff_dct_init(DCTContext *s, int nbits, enum DCTTransformType inverse)
00185 {
00186 int n = 1 << nbits;
00187 int i;
00188
00189 s->nbits = nbits;
00190 s->inverse = inverse;
00191
00192 ff_init_ff_cos_tabs(nbits+2);
00193
00194 s->costab = ff_cos_tabs[nbits+2];
00195
00196 s->csc2 = av_malloc(n/2 * sizeof(FFTSample));
00197
00198 if (ff_rdft_init(&s->rdft, nbits, inverse == DCT_III) < 0) {
00199 av_free(s->csc2);
00200 return -1;
00201 }
00202
00203 for (i = 0; i < n/2; i++)
00204 s->csc2[i] = 0.5 / sin((M_PI / (2*n) * (2*i + 1)));
00205
00206 switch(inverse) {
00207 case DCT_I : s->dct_calc = ff_dct_calc_I_c; break;
00208 case DCT_II : s->dct_calc = ff_dct_calc_II_c ; break;
00209 case DCT_III: s->dct_calc = ff_dct_calc_III_c; break;
00210 case DST_I : s->dct_calc = ff_dst_calc_I_c; break;
00211 }
00212
00213 if (inverse == DCT_II && nbits == 5)
00214 s->dct_calc = dct32_func;
00215
00216 s->dct32 = dct32;
00217 if (HAVE_MMX) ff_dct_init_mmx(s);
00218
00219 return 0;
00220 }
00221
00222 av_cold void ff_dct_end(DCTContext *s)
00223 {
00224 ff_rdft_end(&s->rdft);
00225 av_free(s->csc2);
00226 }