jorbis/src/com/jcraft/jorbis/Lpc.java

   1 /* JOrbis
   2  * Copyright (C) 2000 ymnk, JCraft,Inc.
   3  *
   4  * Written by: 2000 ymnk<ymnk@jcraft.com>
   5  *
   6  * Many thanks to
   7  *   Monty <monty@xiph.org> and
   8  *   The XIPHOPHORUS Company http://www.xiph.org/ .
   9  * JOrbis has been based on their awesome works, Vorbis codec.
  10  *
  11  * This program is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Library General Public License
  13  * as published by the Free Software Foundation; either version 2 of
  14  * the License, or (at your option) any later version.
  15
  16  * This program is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19  * GNU Library General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Library General Public
  22  * License along with this program; if not, write to the Free Software
  23  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  24  */
  25
  26 package com.jcraft.jorbis;
  27
  28 class Lpc{
  29   // en/decode lookups
  30   Drft fft=new Drft();;
  31
  32   int ln;
  33   int m;
  34
  35   // Autocorrelation LPC coeff generation algorithm invented by
  36   // N. Levinson in 1947, modified by J. Durbin in 1959.
  37
  38   // Input : n elements of time doamin data
  39   // Output: m lpc coefficients, excitation energy
  40
  41   static float lpc_from_data(float[] data, float[] lpc,int n,int m){
  42     float[] aut=new float[m+1];
  43     float error;
  44     int i,j;
  45
  46     // autocorrelation, p+1 lag coefficients
  47
  48     j=m+1;
  49     while(j--!=0){
  50       float d=0;
  51       for(i=j;i<n;i++)d+=data[i]*data[i-j];
  52       aut[j]=d;
  53     }
  54
  55     // Generate lpc coefficients from autocorr values
  56
  57     error=aut[0];
  58     /*
  59     if(error==0){
  60       for(int k=0; k<m; k++) lpc[k]=0.0f;
  61       return 0;
  62     }
  63     */
  64
  65     for(i=0;i<m;i++){
  66       float r=-aut[i+1];
  67
  68       if(error==0){
  69         for(int k=0; k<m; k++) lpc[k]=0.0f;
  70         return 0;
  71       }
  72
  73       // Sum up this iteration's reflection coefficient; note that in
  74       // Vorbis we don't save it.  If anyone wants to recycle this code
  75       // and needs reflection coefficients, save the results of 'r' from
  76       // each iteration.
  77
  78       for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
  79       r/=error;
  80
  81       // Update LPC coefficients and total error
  82
  83       lpc[i]=r;
  84       for(j=0;j<i/2;j++){
  85         float tmp=lpc[j];
  86         lpc[j]+=r*lpc[i-1-j];
  87         lpc[i-1-j]+=r*tmp;
  88       }
  89       if(i%2!=0)lpc[j]+=lpc[j]*r;
  90
  91       error*=1.0-r*r;
  92     }
  93
  94     // we need the error value to know how big an impulse to hit the
  95     // filter with later
  96
  97     return error;
  98   }
  99
 100   // Input : n element envelope spectral curve
 101   // Output: m lpc coefficients, excitation energy
 102
 103   float lpc_from_curve(float[] curve, float[] lpc){
 104     int n=ln;
 105     float[] work=new float[n+n];
 106     float fscale=(float)(.5/n);
 107     int i,j;
 108
 109     // input is a real curve. make it complex-real
 110     // This mixes phase, but the LPC generation doesn't care.
 111     for(i=0;i<n;i++){
 112       work[i*2]=curve[i]*fscale;
 113       work[i*2+1]=0;
 114     }
 115     work[n*2-1]=curve[n-1]*fscale;
 116
 117     n*=2;
 118     fft.backward(work);
 119
 120     // The autocorrelation will not be circular.  Shift, else we lose
 121     // most of the power in the edges.
 122
 123     for(i=0,j=n/2;i<n/2;){
 124       float temp=work[i];
 125       work[i++]=work[j];
 126       work[j++]=temp;
 127     }
 128
 129     return(lpc_from_data(work,lpc,n,m));
 130   }
 131
 132   void init(int mapped, int m){
 133     //memset(l,0,sizeof(lpc_lookup));
 134
 135     ln=mapped;
 136     this.m=m;
 137
 138     // we cheat decoding the LPC spectrum via FFTs
 139     fft.init(mapped*2);
 140   }
 141
 142   void clear(){
 143     fft.clear();
 144   }
 145
 146   static float FAST_HYPOT(float a, float b){
 147     return (float)Math.sqrt((a)*(a) + (b)*(b));
 148   }
 149
 150   // One can do this the long way by generating the transfer function in
 151   // the time domain and taking the forward FFT of the result.  The
 152   // results from direct calculation are cleaner and faster.
 153   //
 154   // This version does a linear curve generation and then later
 155   // interpolates the log curve from the linear curve.
 156
 157   void lpc_to_curve(float[] curve, float[] lpc, float amp){
 158
 159     //memset(curve,0,sizeof(float)*l->ln*2);
 160     for(int i=0; i<ln*2; i++)curve[i]=0.0f;
 161
 162     if(amp==0)return;
 163
 164     for(int i=0;i<m;i++){
 165       curve[i*2+1]=lpc[i]/(4*amp);
 166       curve[i*2+2]=-lpc[i]/(4*amp);
 167     }
 168
 169     fft.backward(curve); // reappropriated ;-)
 170
 171     {
 172       int l2=ln*2;
 173       float unit=(float)(1./amp);
 174       curve[0]=(float)(1./(curve[0]*2+unit));
 175       for(int i=1;i<ln;i++){
 176         float real=(curve[i]+curve[l2-i]);
 177         float imag=(curve[i]-curve[l2-i]);
 178
 179         float a = real + unit;
 180         curve[i] = (float)(1.0 / FAST_HYPOT(a, imag));
 181       }
 182     }
 183   }
 184
 185 /*
 186   // subtract or add an lpc filter to data.  Vorbis doesn't actually use this.
 187
 188   static void lpc_residue(float[] coeff, float[] prime,int m,
 189                           float[] data, int n){
 190
 191     // in: coeff[0...m-1] LPC coefficients
 192     //     prime[0...m-1] initial values
 193     //     data[0...n-1] data samples
 194     // out: data[0...n-1] residuals from LPC prediction
 195
 196     float[] work=new float[m+n];
 197     float y;
 198
 199     if(prime==null){
 200       for(int i=0;i<m;i++){
 201         work[i]=0;
 202       }
 203     }
 204     else{
 205       for(int i=0;i<m;i++){
 206         work[i]=prime[i];
 207       }
 208     }
 209
 210     for(int i=0;i<n;i++){
 211       y=0;
 212       for(int j=0;j<m;j++){
 213         y-=work[i+j]*coeff[m-j-1];
 214       }
 215       work[i+m]=data[i];
 216       data[i]-=y;
 217     }
 218   }
 219
 220   static void lpc_predict(float[] coeff, float[] prime,int m,
 221                           float[] data, int n){
 222
 223     // in: coeff[0...m-1] LPC coefficients
 224     //     prime[0...m-1] initial values (allocated size of n+m-1)
 225     //     data[0...n-1] residuals from LPC prediction
 226     // out: data[0...n-1] data samples
 227
 228     int o,p;
 229     float y;
 230     float[] work=new float[m+n];
 231
 232     if(prime==null){
 233       for(int i=0;i<m;i++){
 234         work[i]=0.f;
 235       }
 236     }
 237     else{
 238       for(int i=0;i<m;i++){
 239         work[i]=prime[i];
 240       }
 241     }
 242
 243     for(int i=0;i<n;i++){
 244       y=data[i];
 245       o=i;
 246       p=m;
 247       for(int j=0;j<m;j++){
 248         y-=work[o++]*coeff[--p];
 249       }
 250       data[i]=work[o]=y;
 251     }
 252   }
 253 */
 254 }