/* -*-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 Lpc {
	// en/decode lookups
	Drft fft = new Drft();;

	int ln;
	int m;

	// Autocorrelation LPC coeff generation algorithm invented by
	// N. Levinson in 1947, modified by J. Durbin in 1959.

	// Input : n elements of time doamin data
	// Output: m lpc coefficients, excitation energy

	static float lpc_from_data(float[] data, float[] lpc, int n, int m) {
		float[] aut = new float[m + 1];
		float error;
		int i, j;

		// autocorrelation, p+1 lag coefficients

		j = m + 1;
		while (j-- != 0) {
			float d = 0;
			for (i = j; i < n; i++)
				d += data[i] * data[i - j];
			aut[j] = d;
		}

		// Generate lpc coefficients from autocorr values

		error = aut[0];
		/*
		 * if(error==0){ for(int k=0; k<m; k++) lpc[k]=0.0f; return 0; }
		 */

		for (i = 0; i < m; i++) {
			float r = -aut[i + 1];

			if (error == 0) {
				for (int k = 0; k < m; k++)
					lpc[k] = 0.0f;
				return 0;
			}

			// Sum up this iteration's reflection coefficient; note that in
			// Vorbis we don't save it. If anyone wants to recycle this code
			// and needs reflection coefficients, save the results of 'r' from
			// each iteration.

			for (j = 0; j < i; j++)
				r -= lpc[j] * aut[i - j];
			r /= error;

			// Update LPC coefficients and total error

			lpc[i] = r;
			for (j = 0; j < i / 2; j++) {
				float tmp = lpc[j];
				lpc[j] += r * lpc[i - 1 - j];
				lpc[i - 1 - j] += r * tmp;
			}
			if (i % 2 != 0)
				lpc[j] += lpc[j] * r;

			error *= 1.0 - r * r;
		}

		// we need the error value to know how big an impulse to hit the
		// filter with later

		return error;
	}

	// Input : n element envelope spectral curve
	// Output: m lpc coefficients, excitation energy

	float lpc_from_curve(float[] curve, float[] lpc) {
		int n = ln;
		float[] work = new float[n + n];
		float fscale = (float) (.5 / n);
		int i, j;

		// input is a real curve. make it complex-real
		// This mixes phase, but the LPC generation doesn't care.
		for (i = 0; i < n; i++) {
			work[i * 2] = curve[i] * fscale;
			work[i * 2 + 1] = 0;
		}
		work[n * 2 - 1] = curve[n - 1] * fscale;

		n *= 2;
		fft.backward(work);

		// The autocorrelation will not be circular. Shift, else we lose
		// most of the power in the edges.

		for (i = 0, j = n / 2; i < n / 2;) {
			float temp = work[i];
			work[i++] = work[j];
			work[j++] = temp;
		}

		return (lpc_from_data(work, lpc, n, m));
	}

	void init(int mapped, int m) {
		ln = mapped;
		this.m = m;

		// we cheat decoding the LPC spectrum via FFTs
		fft.init(mapped * 2);
	}

	void clear() {
		fft.clear();
	}

	static float FAST_HYPOT(float a, float b) {
		return (float) Math.sqrt((a) * (a) + (b) * (b));
	}

	// One can do this the long way by generating the transfer function in
	// the time domain and taking the forward FFT of the result. The
	// results from direct calculation are cleaner and faster.
	//
	// This version does a linear curve generation and then later
	// interpolates the log curve from the linear curve.

	void lpc_to_curve(float[] curve, float[] lpc, float amp) {

		for (int i = 0; i < ln * 2; i++)
			curve[i] = 0.0f;

		if (amp == 0)
			return;

		for (int i = 0; i < m; i++) {
			curve[i * 2 + 1] = lpc[i] / (4 * amp);
			curve[i * 2 + 2] = -lpc[i] / (4 * amp);
		}

		fft.backward(curve);

		{
			int l2 = ln * 2;
			float unit = (float) (1. / amp);
			curve[0] = (float) (1. / (curve[0] * 2 + unit));
			for (int i = 1; i < ln; i++) {
				float real = (curve[i] + curve[l2 - i]);
				float imag = (curve[i] - curve[l2 - i]);

				float a = real + unit;
				curve[i] = (float) (1.0 / FAST_HYPOT(a, imag));
			}
		}
	}
}