eaglercraft-1.8/sources/teavm/java/com/jcraft/jorbis/Mdct.java

/* -*-mode:java; c-basic-offset:2; indent-tabs-mode:nil -*- */
/* JOrbis
 * Copyright (C) 2000 ymnk, JCraft,Inc.
 *  
 * Written by: 2000 ymnk<ymnk@jcraft.com>
 *   
 * Many thanks to 
 *   Monty <monty@xiph.org> and 
 *   The XIPHOPHORUS Company http://www.xiph.org/ .
 * JOrbis has been based on their awesome works, Vorbis codec.
 *   
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public License
 * as published by the Free Software Foundation; either version 2 of
 * the License, or (at your option) any later version.
   
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

package com.jcraft.jorbis;

class Mdct {

	int n;
	int log2n;

	float[] trig;
	int[] bitrev;

	float scale;

	void init(int n) {
		bitrev = new int[n / 4];
		trig = new float[n + n / 4];

		log2n = (int) Math.rint(Math.log(n) / Math.log(2));
		this.n = n;

		int AE = 0;
		int AO = 1;
		int BE = AE + n / 2;
		int BO = BE + 1;
		int CE = BE + n / 2;
		int CO = CE + 1;
		// trig lookups...
		for (int i = 0; i < n / 4; i++) {
			trig[AE + i * 2] = (float) Math.cos((Math.PI / n) * (4 * i));
			trig[AO + i * 2] = (float) -Math.sin((Math.PI / n) * (4 * i));
			trig[BE + i * 2] = (float) Math.cos((Math.PI / (2 * n)) * (2 * i + 1));
			trig[BO + i * 2] = (float) Math.sin((Math.PI / (2 * n)) * (2 * i + 1));
		}
		for (int i = 0; i < n / 8; i++) {
			trig[CE + i * 2] = (float) Math.cos((Math.PI / n) * (4 * i + 2));
			trig[CO + i * 2] = (float) -Math.sin((Math.PI / n) * (4 * i + 2));
		}

		{
			int mask = (1 << (log2n - 1)) - 1;
			int msb = 1 << (log2n - 2);
			for (int i = 0; i < n / 8; i++) {
				int acc = 0;
				for (int j = 0; msb >>> j != 0; j++)
					if (((msb >>> j) & i) != 0)
						acc |= 1 << j;
				bitrev[i * 2] = ((~acc) & mask);
				// bitrev[i*2]=((~acc)&mask)-1;
				bitrev[i * 2 + 1] = acc;
			}
		}
		scale = 4.f / n;
	}

	void clear() {
	}

	void forward(float[] in, float[] out) {
	}

	float[] _x = new float[1024];
	float[] _w = new float[1024];

	synchronized void backward(float[] in, float[] out) {
		if (_x.length < n / 2) {
			_x = new float[n / 2];
		}
		if (_w.length < n / 2) {
			_w = new float[n / 2];
		}
		float[] x = _x;
		float[] w = _w;
		int n2 = n >>> 1;
		int n4 = n >>> 2;
		int n8 = n >>> 3;

		// rotate + step 1
		{
			int inO = 1;
			int xO = 0;
			int A = n2;

			int i;
			for (i = 0; i < n8; i++) {
				A -= 2;
				x[xO++] = -in[inO + 2] * trig[A + 1] - in[inO] * trig[A];
				x[xO++] = in[inO] * trig[A + 1] - in[inO + 2] * trig[A];
				inO += 4;
			}

			inO = n2 - 4;

			for (i = 0; i < n8; i++) {
				A -= 2;
				x[xO++] = in[inO] * trig[A + 1] + in[inO + 2] * trig[A];
				x[xO++] = in[inO] * trig[A] - in[inO + 2] * trig[A + 1];
				inO -= 4;
			}
		}

		float[] xxx = mdct_kernel(x, w, n, n2, n4, n8);
		int xx = 0;

		// step 8

		{
			int B = n2;
			int o1 = n4, o2 = o1 - 1;
			int o3 = n4 + n2, o4 = o3 - 1;

			for (int i = 0; i < n4; i++) {
				float temp1 = (xxx[xx] * trig[B + 1] - xxx[xx + 1] * trig[B]);
				float temp2 = -(xxx[xx] * trig[B] + xxx[xx + 1] * trig[B + 1]);

				out[o1] = -temp1;
				out[o2] = temp1;
				out[o3] = temp2;
				out[o4] = temp2;

				o1++;
				o2--;
				o3++;
				o4--;
				xx += 2;
				B += 2;
			}
		}
	}

	private float[] mdct_kernel(float[] x, float[] w, int n, int n2, int n4, int n8) {
		// step 2

		int xA = n4;
		int xB = 0;
		int w2 = n4;
		int A = n2;

		for (int i = 0; i < n4;) {
			float x0 = x[xA] - x[xB];
			float x1;
			w[w2 + i] = x[xA++] + x[xB++];

			x1 = x[xA] - x[xB];
			A -= 4;

			w[i++] = x0 * trig[A] + x1 * trig[A + 1];
			w[i] = x1 * trig[A] - x0 * trig[A + 1];

			w[w2 + i] = x[xA++] + x[xB++];
			i++;
		}

		// step 3

		{
			for (int i = 0; i < log2n - 3; i++) {
				int k0 = n >>> (i + 2);
				int k1 = 1 << (i + 3);
				int wbase = n2 - 2;

				A = 0;
				float[] temp;

				for (int r = 0; r < (k0 >>> 2); r++) {
					int w1 = wbase;
					w2 = w1 - (k0 >> 1);
					float AEv = trig[A], wA;
					float AOv = trig[A + 1], wB;
					wbase -= 2;

					k0++;
					for (int s = 0; s < (2 << i); s++) {
						wB = w[w1] - w[w2];
						x[w1] = w[w1] + w[w2];

						wA = w[++w1] - w[++w2];
						x[w1] = w[w1] + w[w2];

						x[w2] = wA * AEv - wB * AOv;
						x[w2 - 1] = wB * AEv + wA * AOv;

						w1 -= k0;
						w2 -= k0;
					}
					k0--;
					A += k1;
				}

				temp = w;
				w = x;
				x = temp;
			}
		}

		// step 4, 5, 6, 7
		{
			int C = n;
			int bit = 0;
			int x1 = 0;
			int x2 = n2 - 1;

			for (int i = 0; i < n8; i++) {
				int t1 = bitrev[bit++];
				int t2 = bitrev[bit++];

				float wA = w[t1] - w[t2 + 1];
				float wB = w[t1 - 1] + w[t2];
				float wC = w[t1] + w[t2 + 1];
				float wD = w[t1 - 1] - w[t2];

				float wACE = wA * trig[C];
				float wBCE = wB * trig[C++];
				float wACO = wA * trig[C];
				float wBCO = wB * trig[C++];

				x[x1++] = (wC + wACO + wBCE) * .5f;
				x[x2--] = (-wD + wBCO - wACE) * .5f;
				x[x1++] = (wD + wBCO - wACE) * .5f;
				x[x2--] = (wC - wACO - wBCE) * .5f;
			}
		}
		return (x);
	}
}