/* 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){
    //memset(l,0,sizeof(lpc_lookup));

    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){

    //memset(curve,0,sizeof(float)*l->ln*2);
    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); // reappropriated ;-)

    {
      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));
      }
    }
  }  

/*
  // subtract or add an lpc filter to data.  Vorbis doesn't actually use this.

  static void lpc_residue(float[] coeff, float[] prime,int m,
			  float[] data, int n){

    // in: coeff[0...m-1] LPC coefficients 
    //     prime[0...m-1] initial values 
    //     data[0...n-1] data samples 
    // out: data[0...n-1] residuals from LPC prediction

    float[] work=new float[m+n];
    float y;

    if(prime==null){
      for(int i=0;i<m;i++){
	work[i]=0;
      }
    }
    else{
      for(int i=0;i<m;i++){
	work[i]=prime[i];
      }
    }

    for(int i=0;i<n;i++){
      y=0;
      for(int j=0;j<m;j++){
	y-=work[i+j]*coeff[m-j-1];
      }
      work[i+m]=data[i];
      data[i]-=y;
    }
  }

  static void lpc_predict(float[] coeff, float[] prime,int m,
			  float[] data, int n){

    // in: coeff[0...m-1] LPC coefficients 
    //     prime[0...m-1] initial values (allocated size of n+m-1)
    //     data[0...n-1] residuals from LPC prediction   
    // out: data[0...n-1] data samples

    int o,p;
    float y;
    float[] work=new float[m+n];

    if(prime==null){
      for(int i=0;i<m;i++){
	work[i]=0.f;
      }
    }
    else{
      for(int i=0;i<m;i++){
	work[i]=prime[i];
      }
    }

    for(int i=0;i<n;i++){
      y=data[i];
      o=i;
      p=m;
      for(int j=0;j<m;j++){
	y-=work[o++]*coeff[--p];
      }
      data[i]=work[o]=y;
    }
  }
*/
}