Lpc.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 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));
      }
    }
  }
}