zenilib  0.5.3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
lpc.c
Go to the documentation of this file.
1 /********************************************************************
2  * *
3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
7  * *
8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 *
9  * by the Xiph.Org Foundation http://www.xiph.org/ *
10  * *
11  ********************************************************************
12 
13  function: LPC low level routines
14  last mod: $Id: lpc.c 16227 2009-07-08 06:58:46Z xiphmont $
15 
16  ********************************************************************/
17 
18 /* Some of these routines (autocorrelator, LPC coefficient estimator)
19  are derived from code written by Jutta Degener and Carsten Bormann;
20  thus we include their copyright below. The entirety of this file
21  is freely redistributable on the condition that both of these
22  copyright notices are preserved without modification. */
23 
24 /* Preserved Copyright: *********************************************/
25 
26 /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
27 Technische Universita"t Berlin
28 
29 Any use of this software is permitted provided that this notice is not
30 removed and that neither the authors nor the Technische Universita"t
31 Berlin are deemed to have made any representations as to the
32 suitability of this software for any purpose nor are held responsible
33 for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
34 THIS SOFTWARE.
35 
36 As a matter of courtesy, the authors request to be informed about uses
37 this software has found, about bugs in this software, and about any
38 improvements that may be of general interest.
39 
40 Berlin, 28.11.1994
41 Jutta Degener
42 Carsten Bormann
43 
44 *********************************************************************/
45 
46 #include <stdlib.h>
47 #include <string.h>
48 #include <math.h>
49 #include "os.h"
50 #include "smallft.h"
51 #include "lpc.h"
52 #include "scales.h"
53 #include "misc.h"
54 
55 /* Autocorrelation LPC coeff generation algorithm invented by
56  N. Levinson in 1947, modified by J. Durbin in 1959. */
57 
58 /* Input : n elements of time doamin data
59  Output: m lpc coefficients, excitation energy */
60 
61 float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
62  double *aut=alloca(sizeof(*aut)*(m+1));
63  double *lpc=alloca(sizeof(*lpc)*(m));
64  double error;
65  double epsilon;
66  int i,j;
67 
68  /* autocorrelation, p+1 lag coefficients */
69  j=m+1;
70  while(j--){
71  double d=0; /* double needed for accumulator depth */
72  for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
73  aut[j]=d;
74  }
75 
76  /* Generate lpc coefficients from autocorr values */
77 
78  /* set our noise floor to about -100dB */
79  error=aut[0] * (1. + 1e-10);
80  epsilon=1e-9*aut[0]+1e-10;
81 
82  for(i=0;i<m;i++){
83  double r= -aut[i+1];
84 
85  if(error<epsilon){
86  memset(lpc+i,0,(m-i)*sizeof(*lpc));
87  goto done;
88  }
89 
90  /* Sum up this iteration's reflection coefficient; note that in
91  Vorbis we don't save it. If anyone wants to recycle this code
92  and needs reflection coefficients, save the results of 'r' from
93  each iteration. */
94 
95  for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
96  r/=error;
97 
98  /* Update LPC coefficients and total error */
99 
100  lpc[i]=r;
101  for(j=0;j<i/2;j++){
102  double tmp=lpc[j];
103 
104  lpc[j]+=r*lpc[i-1-j];
105  lpc[i-1-j]+=r*tmp;
106  }
107  if(i&1)lpc[j]+=lpc[j]*r;
108 
109  error*=1.-r*r;
110 
111  }
112 
113  done:
114 
115  /* slightly damp the filter */
116  {
117  double g = .99;
118  double damp = g;
119  for(j=0;j<m;j++){
120  lpc[j]*=damp;
121  damp*=g;
122  }
123  }
124 
125  for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
126 
127  /* we need the error value to know how big an impulse to hit the
128  filter with later */
129 
130  return error;
131 }
132 
133 void vorbis_lpc_predict(float *coeff,float *prime,int m,
134  float *data,long n){
135 
136  /* in: coeff[0...m-1] LPC coefficients
137  prime[0...m-1] initial values (allocated size of n+m-1)
138  out: data[0...n-1] data samples */
139 
140  long i,j,o,p;
141  float y;
142  float *work=alloca(sizeof(*work)*(m+n));
143 
144  if(!prime)
145  for(i=0;i<m;i++)
146  work[i]=0.f;
147  else
148  for(i=0;i<m;i++)
149  work[i]=prime[i];
150 
151  for(i=0;i<n;i++){
152  y=0;
153  o=i;
154  p=m;
155  for(j=0;j<m;j++)
156  y-=work[o++]*coeff[--p];
157 
158  data[i]=work[o]=y;
159  }
160 }
float vorbis_lpc_from_data(float *data, float *lpci, int n, int m)
Definition: lpc.c:61
GLboolean GLboolean g
Definition: glew.h:8736
GLclampf f
Definition: glew.h:3390
void vorbis_lpc_predict(float *coeff, float *prime, int m, float *data, long n)
Definition: lpc.c:133
GLclampd n
Definition: glew.h:7287
return Display return Display Bool Bool int int e
Definition: SDL_x11sym.h:30
int32_t j
Definition: e_log.c:102
#define memset
Definition: SDL_malloc.c:633
return Display return Display Bool Bool int d
Definition: SDL_x11sym.h:30
GLint GLenum GLsizei GLsizei GLsizei GLint GLsizei const GLvoid * data
Definition: gl2ext.h:848
FT_Error error
Definition: cffdrivr.c:407
GLfloat GLfloat p
Definition: glew.h:14938
EGLSurface EGLint EGLint y
Definition: eglext.h:293
GLdouble GLdouble GLdouble r
Definition: glew.h:1392
int i
Definition: pngrutil.c:1377
#define m(i, j)
const GLdouble * m
Definition: glew.h:8385