/*
* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.
* USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS
* GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE
* IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.
*
* THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002
* by the Xiph.Org Foundation http://www.xiph.org/
*/
package org.xiph.libvorbis;
import static org.xiph.libvorbis.vorbis_constants.integer_constants.*;
class mdct_lookup {
int n;
int log2n;
float[] trig;
int[] bitrev;
float scale;
public mdct_lookup() {
}
public mdct_lookup(int _n, int _log2n, float[] _trig, int[] _bitrev, float _scale) {
n = _n;
log2n = _log2n;
trig = (float[]) _trig.clone();
bitrev = (int[]) _bitrev.clone();
scale = _scale;
}
public mdct_lookup(mdct_lookup src) {
this(src.n, src.log2n, src.trig, src.bitrev, src.scale);
}
public void mdct_init(int _n) {
n = _n;
log2n = new Double(Math.rint(Math.log(new Integer(n).floatValue()) / Math.log(2.f))).intValue();
// DATA_TYPE *T=_ogg_malloc(sizeof(*T)*(n+n/4));
trig = new float[n + n / 4];
// int *bitrev=_ogg_malloc(sizeof(*bitrev)*(n/4));
bitrev = new int[n / 4];
int i, j;
int n2 = n >> 1;
// trig lookups...
for (i = 0; i < n / 4; i++) {
trig[i * 2] = new Double(Math.cos((M_PI / n) * (4 * i))).floatValue();
trig[i * 2 + 1] = new Double(-Math.sin((M_PI / n) * (4 * i))).floatValue();
trig[n2 + i * 2] = new Double(Math.cos((M_PI / (2 * n)) * (2 * i + 1))).floatValue();
trig[n2 + i * 2 + 1] = new Double(Math.sin((M_PI / (2 * n)) * (2 * i + 1))).floatValue();
}
for (i = 0; i < n / 8; i++) {
trig[n + i * 2] = new Double(Math.cos((M_PI / n) * (4 * i + 2)) * .5).floatValue();
trig[n + i * 2 + 1] = new Double(-Math.sin((M_PI / n) * (4 * i + 2)) * .5).floatValue();
}
// bitreverse lookup...
{
int mask = (1 << (log2n - 1)) - 1;
int msb = 1 << (log2n - 2);
for (i = 0; i < n / 8; i++) {
int acc = 0;
for (j = 0; (msb >> j) != 0; j++)
if (((msb >> j) & i) != 0)
acc |= 1 << j;
bitrev[i * 2] = ((~acc) & mask) - 1;
bitrev[i * 2 + 1] = acc;
}
}
scale = 4.f / n;
}
// 8 point butterfly (in place, 4 register)
private void mdct_butterfly_8(float[] x, int off) {
float r0 = x[off + 6] + x[off + 2];
float r1 = x[off + 6] - x[off + 2];
float r2 = x[off + 4] + x[off + 0];
float r3 = x[off + 4] - x[off + 0];
x[off + 6] = r0 + r2;
x[off + 4] = r0 - r2;
r0 = x[off + 5] - x[off + 1];
r2 = x[off + 7] - x[off + 3];
x[off + 0] = r1 + r0;
x[off + 2] = r1 - r0;
r0 = x[off + 5] + x[off + 1];
r1 = x[off + 7] + x[off + 3];
x[off + 3] = r2 + r3;
x[off + 1] = r2 - r3;
x[off + 7] = r1 + r0;
x[off + 5] = r1 - r0;
}
// 16 point butterfly (in place, 4 register)
private void mdct_butterfly_16(float[] x, int off) {
float r0 = x[off + 1] - x[off + 9];
float r1 = x[off + 0] - x[off + 8];
x[off + 8] += x[off + 0];
x[off + 9] += x[off + 1];
x[off + 0] = ((r0 + r1) * cPI2_8);
x[off + 1] = ((r0 - r1) * cPI2_8);
r0 = x[off + 3] - x[off + 11];
r1 = x[off + 10] - x[off + 2];
x[off + 10] += x[off + 2];
x[off + 11] += x[off + 3];
x[off + 2] = r0;
x[off + 3] = r1;
r0 = x[off + 12] - x[off + 4];
r1 = x[off + 13] - x[off + 5];
x[off + 12] += x[off + 4];
x[off + 13] += x[off + 5];
x[off + 4] = ((r0 - r1) * cPI2_8);
x[off + 5] = ((r0 + r1) * cPI2_8);
r0 = x[off + 14] - x[off + 6];
r1 = x[off + 15] - x[off + 7];
x[off + 14] += x[off + 6];
x[off + 15] += x[off + 7];
x[off + 6] = r0;
x[off + 7] = r1;
mdct_butterfly_8(x, 0 + off);
mdct_butterfly_8(x, 8 + off);
}
// 32 point butterfly (in place, 4 register)
private void mdct_butterfly_32(float[] x, int off) {
float r0 = x[off + 30] - x[off + 14];
float r1 = x[off + 31] - x[off + 15];
x[off + 30] += x[off + 14];
x[off + 31] += x[off + 15];
x[off + 14] = r0;
x[off + 15] = r1;
r0 = x[off + 28] - x[off + 12];
r1 = x[off + 29] - x[off + 13];
x[off + 28] += x[off + 12];
x[off + 29] += x[off + 13];
x[off + 12] = (r0 * cPI1_8 - r1 * cPI3_8);
x[off + 13] = (r0 * cPI3_8 + r1 * cPI1_8);
r0 = x[off + 26] - x[off + 10];
r1 = x[off + 27] - x[off + 11];
x[off + 26] += x[off + 10];
x[off + 27] += x[off + 11];
x[off + 10] = ((r0 - r1) * cPI2_8);
x[off + 11] = ((r0 + r1) * cPI2_8);
r0 = x[off + 24] - x[off + 8];
r1 = x[off + 25] - x[off + 9];
x[off + 24] += x[off + 8];
x[off + 25] += x[off + 9];
x[off + 8] = (r0 * cPI3_8 - r1 * cPI1_8);
x[off + 9] = (r1 * cPI3_8 + r0 * cPI1_8);
r0 = x[off + 22] - x[off + 6];
r1 = x[off + 7] - x[off + 23];
x[off + 22] += x[off + 6];
x[off + 23] += x[off + 7];
x[off + 6] = r1;
x[off + 7] = r0;
r0 = x[off + 4] - x[off + 20];
r1 = x[off + 5] - x[off + 21];
x[off + 20] += x[off + 4];
x[off + 21] += x[off + 5];
x[off + 4] = (r1 * cPI1_8 + r0 * cPI3_8);
x[off + 5] = (r1 * cPI3_8 - r0 * cPI1_8);
r0 = x[off + 2] - x[off + 18];
r1 = x[off + 3] - x[off + 19];
x[off + 18] += x[off + 2];
x[off + 19] += x[off + 3];
x[off + 2] = ((r1 + r0) * cPI2_8);
x[off + 3] = ((r1 - r0) * cPI2_8);
r0 = x[off + 0] - x[off + 16];
r1 = x[off + 1] - x[off + 17];
x[off + 16] += x[off + 0];
x[off + 17] += x[off + 1];
x[off + 0] = (r1 * cPI3_8 + r0 * cPI1_8);
x[off + 1] = (r1 * cPI1_8 - r0 * cPI3_8);
mdct_butterfly_16(x, 0 + off);
mdct_butterfly_16(x, 16 + off);
}
// N point first stage butterfly (in place, 2 register)
private void mdct_butterfly_first(float[] T, float[] x, int off, int points) {
// float *x1 = x + points - 8;
int x1 = off + points - 8;
// float *x2 = x + (points>>1) - 8;
int x2 = off + (points >> 1) - 8;
int t = 0;
float r0;
float r1;
do {
r0 = x[x1 + 6] - x[x2 + 6];
r1 = x[x1 + 7] - x[x2 + 7];
x[x1 + 6] += x[x2 + 6];
x[x1 + 7] += x[x2 + 7];
x[x2 + 6] = (r1 * T[t + 1] + r0 * T[t + 0]);
x[x2 + 7] = (r1 * T[t + 0] - r0 * T[t + 1]);
r0 = x[x1 + 4] - x[x2 + 4];
r1 = x[x1 + 5] - x[x2 + 5];
x[x1 + 4] += x[x2 + 4];
x[x1 + 5] += x[x2 + 5];
x[x2 + 4] = (r1 * T[t + 5] + r0 * T[t + 4]);
x[x2 + 5] = (r1 * T[t + 4] - r0 * T[t + 5]);
r0 = x[x1 + 2] - x[x2 + 2];
r1 = x[x1 + 3] - x[x2 + 3];
x[x1 + 2] += x[x2 + 2];
x[x1 + 3] += x[x2 + 3];
x[x2 + 2] = (r1 * T[t + 9] + r0 * T[t + 8]);
x[x2 + 3] = (r1 * T[t + 8] - r0 * T[t + 9]);
r0 = x[x1 + 0] - x[x2 + 0];
r1 = x[x1 + 1] - x[x2 + 1];
x[x1 + 0] += x[x2 + 0];
x[x1 + 1] += x[x2 + 1];
x[x2 + 0] = (r1 * T[t + 13] + r0 * T[t + 12]);
x[x2 + 1] = (r1 * T[t + 12] - r0 * T[t + 13]);
x1 -= 8;
x2 -= 8;
t += 16;
} while (x2 >= off);
// } while ( x2 >= x );
}
// N/stage point generic N stage butterfly (in place, 2 register)
private void mdct_butterfly_generic(float[] T, float[] x, int off, int points, int trigint) {
// float *x1 = x + points - 8;
int x1 = off + points - 8;
// float *x2 = x + (points>>1) - 8;
int x2 = off + (points >> 1) - 8;
int t = 0;
float r0;
float r1;
do {
r0 = x[x1 + 6] - x[x2 + 6];
r1 = x[x1 + 7] - x[x2 + 7];
x[x1 + 6] += x[x2 + 6];
x[x1 + 7] += x[x2 + 7];
x[x2 + 6] = (r1 * T[t + 1] + r0 * T[t + 0]);
x[x2 + 7] = (r1 * T[t + 0] - r0 * T[t + 1]);
t += trigint;
r0 = x[x1 + 4] - x[x2 + 4];
r1 = x[x1 + 5] - x[x2 + 5];
x[x1 + 4] += x[x2 + 4];
x[x1 + 5] += x[x2 + 5];
x[x2 + 4] = (r1 * T[t + 1] + r0 * T[t + 0]);
x[x2 + 5] = (r1 * T[t + 0] - r0 * T[t + 1]);
t += trigint;
r0 = x[x1 + 2] - x[x2 + 2];
r1 = x[x1 + 3] - x[x2 + 3];
x[x1 + 2] += x[x2 + 2];
x[x1 + 3] += x[x2 + 3];
x[x2 + 2] = (r1 * T[t + 1] + r0 * T[t + 0]);
x[x2 + 3] = (r1 * T[t + 0] - r0 * T[t + 1]);
t += trigint;
r0 = x[x1 + 0] - x[x2 + 0];
r1 = x[x1 + 1] - x[x2 + 1];
x[x1 + 0] += x[x2 + 0];
x[x1 + 1] += x[x2 + 1];
x[x2 + 0] = (r1 * T[t + 1] + r0 * T[t + 0]);
x[x2 + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
t += trigint;
x1 -= 8;
x2 -= 8;
} while (x2 >= off);
// } while ( x2 >= x );
}
private void mdct_butterflies(float[] x, int off, int points) {
float[] T = trig;
int stages = log2n - 5;
int i, j;
if (--stages > 0) {
mdct_butterfly_first(T, x, off, points);
}
for (i = 1; --stages > 0; i++) {
for (j = 0; j < (1 << i); j++)
mdct_butterfly_generic(T, x, off + ((points >> i) * j), points >> i, 4 << i);
}
for (j = 0; j < points; j += 32)
mdct_butterfly_32(x, off + j);
}
private void mdct_bitreverse(float[] x) {
// float *w0 = x;
// float *w1 = x = w0+(n>>1);
// float *T = init->trig+n;
int bit = 0;
int w0 = 0;
int w1 = (n >> 1);
int xoff = (n >> 1);
int t = n;
float[] T = trig;
do {
// float *x0 = x+bit[0];
// float *x1 = x+bit[1];
int x0 = xoff + bitrev[bit + 0];
int x1 = xoff + bitrev[bit + 1];
float r0 = x[x0 + 1] - x[x1 + 1];
float r1 = x[x0 + 0] + x[x1 + 0];
float r2 = (r1 * T[t + 0] + r0 * T[t + 1]);
float r3 = (r1 * T[t + 1] - r0 * T[t + 0]);
w1 -= 4;
r0 = (x[x0 + 1] + x[x1 + 1]) * .5f;
r1 = (x[x0 + 0] - x[x1 + 0]) * .5f;
x[w0 + 0] = r0 + r2;
x[w1 + 2] = r0 - r2;
x[w0 + 1] = r1 + r3;
x[w1 + 3] = r3 - r1;
x0 = xoff + bitrev[bit + 2];
x1 = xoff + bitrev[bit + 3];
r0 = x[x0 + 1] - x[x1 + 1];
r1 = x[x0 + 0] + x[x1 + 0];
r2 = (r1 * T[t + 2] + r0 * T[t + 3]);
r3 = (r1 * T[t + 3] - r0 * T[t + 2]);
r0 = (x[x0 + 1] + x[x1 + 1]) * .5f;
r1 = (x[x0 + 0] - x[x1 + 0]) * .5f;
x[w0 + 2] = r0 + r2;
x[w1 + 0] = r0 - r2;
x[w0 + 3] = r1 + r3;
x[w1 + 1] = r3 - r1;
t += 4;
bit += 4;
w0 += 4;
} while (w0 < w1);
}
public void mdct_backward(float[] in, float[] out) {
int n2 = n >> 1;
int n4 = n >> 2;
// rotate
// float *iX = in+n2-7;
// float *oX = out+n2+n4;
float[] T = trig;
int iX = n2 - 7;
int oX = n2 + n4;
int t = n4;
do {
oX -= 4;
out[oX + 0] = (-in[iX + 2] * T[t + 3] - in[iX + 0] * T[t + 2]);
out[oX + 1] = (in[iX + 0] * T[t + 3] - in[iX + 2] * T[t + 2]);
out[oX + 2] = (-in[iX + 6] * T[t + 1] - in[iX + 4] * T[t + 0]);
out[oX + 3] = (in[iX + 4] * T[t + 1] - in[iX + 6] * T[t + 0]);
iX -= 8;
t += 4;
} while (iX >= 0);
// float *iX = in+n2-8;
// float *oX = out+n2+n4;
iX = n2 - 8;
oX = n2 + n4;
t = n4;
do {
t -= 4;
out[oX + 0] = (in[iX + 4] * T[t + 3] + in[iX + 6] * T[t + 2]);
out[oX + 1] = (in[iX + 4] * T[t + 2] - in[iX + 6] * T[t + 3]);
out[oX + 2] = (in[iX + 0] * T[t + 1] + in[iX + 2] * T[t + 0]);
out[oX + 3] = (in[iX + 0] * T[t + 0] - in[iX + 2] * T[t + 1]);
iX -= 8;
oX += 4;
} while (iX >= 0);
mdct_butterflies(out, n2, n2);
mdct_bitreverse(out);
// roatate + window
// float *oX1=out+n2+n4;
// float *oX2=out+n2+n4;
// float *iX =out;
// T =init->trig+n2;
int oX1 = n2 + n4;
int oX2 = n2 + n4;
iX = 0;
t = n2;
do {
oX1 -= 4;
out[oX1 + 3] = (out[iX + 0] * T[t + 1] - out[iX + 1] * T[t + 0]);
out[oX2 + 0] = -(out[iX + 0] * T[t + 0] + out[iX + 1] * T[t + 1]);
out[oX1 + 2] = (out[iX + 2] * T[t + 3] - out[iX + 3] * T[t + 2]);
out[oX2 + 1] = -(out[iX + 2] * T[t + 2] + out[iX + 3] * T[t + 3]);
out[oX1 + 1] = (out[iX + 4] * T[t + 5] - out[iX + 5] * T[t + 4]);
out[oX2 + 2] = -(out[iX + 4] * T[t + 4] + out[iX + 5] * T[t + 5]);
out[oX1 + 0] = (out[iX + 6] * T[t + 7] - out[iX + 7] * T[t + 6]);
out[oX2 + 3] = -(out[iX + 6] * T[t + 6] + out[iX + 7] * T[t + 7]);
oX2 += 4;
iX += 8;
t += 8;
} while (iX < oX1);
// float *oX1=out+n4;
// float *oX2=oX1;
// float *iX =out+n2+n4;
iX = n2 + n4;
oX1 = n4;
oX2 = oX1;
do {
oX1 -= 4;
iX -= 4;
out[oX2 + 0] = -(out[oX1 + 3] = out[iX + 3]);
out[oX2 + 1] = -(out[oX1 + 2] = out[iX + 2]);
out[oX2 + 2] = -(out[oX1 + 1] = out[iX + 1]);
out[oX2 + 3] = -(out[oX1 + 0] = out[iX + 0]);
oX2 += 4;
} while (oX2 < iX);
// float *oX1=out+n2+n4;
// float *oX2=out+n2;
// float *iX =out+n2+n4;
iX = n2 + n4;
oX1 = n2 + n4;
oX2 = n2;
do {
oX1 -= 4;
out[oX1 + 0] = out[iX + 3];
out[oX1 + 1] = out[iX + 2];
out[oX1 + 2] = out[iX + 1];
out[oX1 + 3] = out[iX + 0];
iX += 4;
} while (oX1 > oX2);
}
public void mdct_forward(float[] in, float[] out) {
int n2 = n >> 1;
int n4 = n >> 2;
int n8 = n >> 3;
float[] T = trig;
// float *w=alloca(n*sizeof(*w)); // forward needs working space
float[] w = new float[n];
// float *w2=w+n2;
int w2 = n2;
// rotate
// window + rotate + step 1
float r0;
float r1;
// float *x0=in+n2+n4;
// float *x1=x0+1;
// float *T=init->trig+n2;
int x0 = n2 + n4;
int x1 = x0 + 1;
int t = n2;
int i = 0;
for (i = 0; i < n8; i += 2) {
x0 -= 4;
t -= 2;
r0 = in[x0 + 2] + in[x1 + 0];
r1 = in[x0 + 0] + in[x1 + 2];
w[w2 + i] = (r1 * T[t + 1] + r0 * T[t + 0]);
w[w2 + i + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
x1 += 4;
}
// x1=in+1;
x1 = 1;
for (; i < n2 - n8; i += 2) {
t -= 2;
x0 -= 4;
r0 = in[x0 + 2] - in[x1 + 0];
r1 = in[x0 + 0] - in[x1 + 2];
w[w2 + i] = (r1 * T[t + 1] + r0 * T[t + 0]);
w[w2 + i + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
x1 += 4;
}
// x0=in+n;
x0 = n;
for (; i < n2; i += 2) {
t -= 2;
x0 -= 4;
r0 = -in[x0 + 2] - in[x1 + 0];
r1 = -in[x0 + 0] - in[x1 + 2];
w[w2 + i] = (r1 * T[t + 1] + r0 * T[t + 0]);
w[w2 + i + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
x1 += 4;
}
mdct_butterflies(w, n2, n2);
mdct_bitreverse(w);
// roatate + window
// T=init->trig+n2;
// x0=out+n2;
t = n2;
x0 = n2;
int w1 = 0;
for (i = 0; i < n4; i++) {
x0--;
out[i] = ((w[w1 + 0] * T[t + 0] + w[w1 + 1] * T[t + 1]) * scale);
out[x0 + 0] = ((w[w1 + 0] * T[t + 1] - w[w1 + 1] * T[t + 0]) * scale);
w1 += 2;
t += 2;
}
}
public void mdct_forward_offset(float[] in, int offset, float[] out, int out_offset) {
int n2 = n >> 1;
int n4 = n >> 2;
int n8 = n >> 3;
float[] T = trig;
// float *w=alloca(n*sizeof(*w)); // forward needs working space
float[] w = new float[n];
// float *w2=w+n2;
int w2 = n2;
// rotate
// window + rotate + step 1
float r0;
float r1;
// float *x0=in+n2+n4;
// float *x1=x0+1;
// float *T=init->trig+n2;
int x0 = n2 + n4;
int x1 = x0 + 1;
int t = n2;
int i = 0;
for (i = 0; i < n8; i += 2) {
x0 -= 4;
t -= 2;
r0 = in[offset + x0 + 2] + in[offset + x1 + 0];
r1 = in[offset + x0 + 0] + in[offset + x1 + 2];
w[w2 + i] = (r1 * T[t + 1] + r0 * T[t + 0]);
w[w2 + i + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
x1 += 4;
}
// x1=in+1;
x1 = 1;
for (; i < n2 - n8; i += 2) {
t -= 2;
x0 -= 4;
r0 = in[offset + x0 + 2] - in[offset + x1 + 0];
r1 = in[offset + x0 + 0] - in[offset + x1 + 2];
w[w2 + i] = (r1 * T[t + 1] + r0 * T[t + 0]);
w[w2 + i + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
x1 += 4;
}
// x0=in+n;
x0 = n;
for (; i < n2; i += 2) {
t -= 2;
x0 -= 4;
r0 = -in[offset + x0 + 2] - in[offset + x1 + 0];
r1 = -in[offset + x0 + 0] - in[offset + x1 + 2];
w[w2 + i] = (r1 * T[t + 1] + r0 * T[t + 0]);
w[w2 + i + 1] = (r1 * T[t + 0] - r0 * T[t + 1]);
x1 += 4;
}
mdct_butterflies(w, n2, n2);
mdct_bitreverse(w);
// roatate + window
// T=init->trig+n2;
// x0=out+n2;
t = n2;
x0 = n2;
int w1 = 0;
for (i = 0; i < n4; i++) {
x0--;
out[out_offset + i] = ((w[w1 + 0] * T[t + 0] + w[w1 + 1] * T[t + 1]) * scale);
out[out_offset + x0 + 0] = ((w[w1 + 0] * T[t + 1] - w[w1 + 1] * T[t + 0]) * scale);
w1 += 2;
t += 2;
}
}
}