Line data Source code
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$
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 0 : float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
62 0 : double *aut=alloca(sizeof(*aut)*(m+1));
63 0 : double *lpc=alloca(sizeof(*lpc)*(m));
64 : double error;
65 : double epsilon;
66 : int i,j;
67 :
68 : /* autocorrelation, p+1 lag coefficients */
69 0 : j=m+1;
70 0 : while(j--){
71 0 : double d=0; /* double needed for accumulator depth */
72 0 : for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
73 0 : aut[j]=d;
74 : }
75 :
76 : /* Generate lpc coefficients from autocorr values */
77 :
78 : /* set our noise floor to about -100dB */
79 0 : error=aut[0] * (1. + 1e-10);
80 0 : epsilon=1e-9*aut[0]+1e-10;
81 :
82 0 : for(i=0;i<m;i++){
83 0 : double r= -aut[i+1];
84 :
85 0 : if(error<epsilon){
86 0 : memset(lpc+i,0,(m-i)*sizeof(*lpc));
87 0 : 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 0 : for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
96 0 : r/=error;
97 :
98 : /* Update LPC coefficients and total error */
99 :
100 0 : lpc[i]=r;
101 0 : for(j=0;j<i/2;j++){
102 0 : double tmp=lpc[j];
103 :
104 0 : lpc[j]+=r*lpc[i-1-j];
105 0 : lpc[i-1-j]+=r*tmp;
106 : }
107 0 : if(i&1)lpc[j]+=lpc[j]*r;
108 :
109 0 : error*=1.-r*r;
110 :
111 : }
112 :
113 : done:
114 :
115 : /* slightly damp the filter */
116 : {
117 0 : double g = .99;
118 0 : double damp = g;
119 0 : for(j=0;j<m;j++){
120 0 : lpc[j]*=damp;
121 0 : damp*=g;
122 : }
123 : }
124 :
125 0 : 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 0 : return error;
131 : }
132 :
133 0 : 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 0 : float *work=alloca(sizeof(*work)*(m+n));
143 :
144 0 : if(!prime)
145 0 : for(i=0;i<m;i++)
146 0 : work[i]=0.f;
147 : else
148 0 : for(i=0;i<m;i++)
149 0 : work[i]=prime[i];
150 :
151 0 : for(i=0;i<n;i++){
152 0 : y=0;
153 0 : o=i;
154 0 : p=m;
155 0 : for(j=0;j<m;j++)
156 0 : y-=work[o++]*coeff[--p];
157 :
158 0 : data[i]=work[o]=y;
159 : }
160 0 : }
|
