001/* JOrbis 002 * Copyright (C) 2000 ymnk, JCraft,Inc. 003 * 004 * Written by: 2000 ymnk<ymnk@jcaft.com> 005 * 006 * Many thanks to 007 * Monty <monty@xiph.org> and 008 * The XIPHOPHORUS Company http://www.xiph.org/ . 009 * JOrbis has been based on their awesome works, Vorbis codec. 010 * 011 * This program is free software; you can redistribute it and/or 012 * modify it under the terms of the GNU Library General Public License 013 * as published by the Free Software Foundation; either version 2 of 014 * the License, or (at your option) any later version. 015 016 * This program is distributed in the hope that it will be useful, 017 * but WITHOUT ANY WARRANTY; without even the implied warranty of 018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 019 * GNU Library General Public License for more details. 020 * 021 * You should have received a copy of the GNU Library General Public 022 * License along with this program; if not, write to the Free Software 023 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 024 */ 025 026package com.jcraft.jorbis; 027 028class Lpc{ 029 // en/decode lookups 030 Drft fft=new Drft();; 031 032 int ln; 033 int m; 034 035 // Autocorrelation LPC coeff generation algorithm invented by 036 // N. Levinson in 1947, modified by J. Durbin in 1959. 037 038 // Input : n elements of time doamin data 039 // Output: m lpc coefficients, excitation energy 040 041 static float lpc_from_data(float[] data, float[] lpc,int n,int m){ 042 float[] aut=new float[m+1]; 043 float error; 044 int i,j; 045 046 // autocorrelation, p+1 lag coefficients 047 048 j=m+1; 049 while(j--!=0){ 050 float d=0; 051 for(i=j;i<n;i++)d+=data[i]*data[i-j]; 052 aut[j]=d; 053 } 054 055 // Generate lpc coefficients from autocorr values 056 057 error=aut[0]; 058 /* 059 if(error==0){ 060 for(int k=0; k<m; k++) lpc[k]=0.0f; 061 return 0; 062 } 063 */ 064 065 for(i=0;i<m;i++){ 066 float r=-aut[i+1]; 067 068 if(error==0){ 069 for(int k=0; k<m; k++) lpc[k]=0.0f; 070 return 0; 071 } 072 073 // Sum up this iteration's reflection coefficient; note that in 074 // Vorbis we don't save it. If anyone wants to recycle this code 075 // and needs reflection coefficients, save the results of 'r' from 076 // each iteration. 077 078 for(j=0;j<i;j++)r-=lpc[j]*aut[i-j]; 079 r/=error; 080 081 // Update LPC coefficients and total error 082 083 lpc[i]=r; 084 for(j=0;j<i/2;j++){ 085 float tmp=lpc[j]; 086 lpc[j]+=r*lpc[i-1-j]; 087 lpc[i-1-j]+=r*tmp; 088 } 089 if(i%2!=0)lpc[j]+=lpc[j]*r; 090 091 error*=1.0-r*r; 092 } 093 094 // we need the error value to know how big an impulse to hit the 095 // filter with later 096 097 return error; 098 } 099 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}