2 * Copyright (C) 2000 ymnk, JCraft,Inc.
4 * Written by: 2000 ymnk<ymnk@jcraft.com>
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.
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.
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.
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.
26 package com.jcraft.jorbis;
35 // Autocorrelation LPC coeff generation algorithm invented by
36 // N. Levinson in 1947, modified by J. Durbin in 1959.
38 // Input : n elements of time doamin data
39 // Output: m lpc coefficients, excitation energy
41 static float lpc_from_data(float[] data, float[] lpc,int n,int m){
42 float[] aut=new float[m+1];
46 // autocorrelation, p+1 lag coefficients
51 for(i=j;i<n;i++)d+=data[i]*data[i-j];
55 // Generate lpc coefficients from autocorr values
60 for(int k=0; k<m; k++) lpc[k]=0.0f;
69 for(int k=0; k<m; k++) lpc[k]=0.0f;
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
78 for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
81 // Update LPC coefficients and total error
89 if(i%2!=0)lpc[j]+=lpc[j]*r;
94 // we need the error value to know how big an impulse to hit the
100 // Input : n element envelope spectral curve
101 // Output: m lpc coefficients, excitation energy
103 float lpc_from_curve(float[] curve, float[] lpc){
105 float[] work=new float[n+n];
106 float fscale=(float)(.5/n);
109 // input is a real curve. make it complex-real
110 // This mixes phase, but the LPC generation doesn't care.
112 work[i*2]=curve[i]*fscale;
115 work[n*2-1]=curve[n-1]*fscale;
120 // The autocorrelation will not be circular. Shift, else we lose
121 // most of the power in the edges.
123 for(i=0,j=n/2;i<n/2;){
129 return(lpc_from_data(work,lpc,n,m));
132 void init(int mapped, int m){
133 //memset(l,0,sizeof(lpc_lookup));
138 // we cheat decoding the LPC spectrum via FFTs
146 static float FAST_HYPOT(float a, float b){
147 return (float)Math.sqrt((a)*(a) + (b)*(b));
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.
154 // This version does a linear curve generation and then later
155 // interpolates the log curve from the linear curve.
157 void lpc_to_curve(float[] curve, float[] lpc, float amp){
159 //memset(curve,0,sizeof(float)*l->ln*2);
160 for(int i=0; i<ln*2; i++)curve[i]=0.0f;
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);
169 fft.backward(curve); // reappropriated ;-)
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]);
179 float a = real + unit;
180 curve[i] = (float)(1.0 / FAST_HYPOT(a, imag));
186 // subtract or add an lpc filter to data. Vorbis doesn't actually use this.
188 static void lpc_residue(float[] coeff, float[] prime,int m,
189 float[] data, int n){
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
196 float[] work=new float[m+n];
200 for(int i=0;i<m;i++){
205 for(int i=0;i<m;i++){
210 for(int i=0;i<n;i++){
212 for(int j=0;j<m;j++){
213 y-=work[i+j]*coeff[m-j-1];
220 static void lpc_predict(float[] coeff, float[] prime,int m,
221 float[] data, int n){
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
230 float[] work=new float[m+n];
233 for(int i=0;i<m;i++){
238 for(int i=0;i<m;i++){
243 for(int i=0;i<n;i++){
247 for(int j=0;j<m;j++){
248 y-=work[o++]*coeff[--p];