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 "internal/stack_alloc.h"
50 #include "internal/lpc.h"
51
52 /* Autocorrelation LPC coeff generation algorithm invented by
53 N. Levinson in 1947, modified by J. Durbin in 1959. */
54
55 /* Input : n elements of time doamin data
56 Output: m lpc coefficients, excitation energy */
57
vorbis_lpc_from_data(float * data,float * lpci,int n,int m)58 float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
59 double *aut=alloca(sizeof(*aut)*(m+1));
60 double *lpc=alloca(sizeof(*lpc)*(m));
61 double error;
62 double epsilon;
63 int i,j;
64
65 /* autocorrelation, p+1 lag coefficients */
66 j=m+1;
67 while(j--){
68 double d=0; /* double needed for accumulator depth */
69 for(i=j;i<n;i++)d+=(double)data[i]*data[(i-j)];
70 aut[j]=d;
71 }
72
73 /* Generate lpc coefficients from autocorr values */
74
75 /* set our noise floor to about -100dB */
76 error=aut[0] * (1. + 1e-10);
77 epsilon=1e-9*aut[0]+1e-10;
78
79 for(i=0;i<m;i++){
80 double r= -aut[i+1];
81
82 if(error<epsilon){
83 memset(lpc+i,0,(m-i)*sizeof(*lpc));
84 goto done;
85 }
86
87 /* Sum up this iteration's reflection coefficient; note that in
88 Vorbis we don't save it. If anyone wants to recycle this code
89 and needs reflection coefficients, save the results of 'r' from
90 each iteration. */
91
92 for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
93 r/=error;
94
95 /* Update LPC coefficients and total error */
96
97 lpc[i]=r;
98 for(j=0;j<i/2;j++){
99 double tmp=lpc[j];
100
101 lpc[j]+=r*lpc[i-1-j];
102 lpc[i-1-j]+=r*tmp;
103 }
104 if(i&1)lpc[j]+=lpc[j]*r;
105
106 error*=1.-r*r;
107
108 }
109
110 done:
111
112 /* slightly damp the filter */
113 {
114 double g = .99;
115 double damp = g;
116 for(j=0;j<m;j++){
117 lpc[j]*=damp;
118 damp*=g;
119 }
120 }
121
122 for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
123
124 /* we need the error value to know how big an impulse to hit the
125 filter with later */
126
127 return error;
128 }
129
vorbis_lpc_predict(float * coeff,float * prime,int m,float * data,long n)130 void vorbis_lpc_predict(float *coeff,float *prime,int m,
131 float *data,long n){
132
133 /* in: coeff[0...m-1] LPC coefficients
134 prime[0...m-1] initial values (allocated size of n+m-1)
135 out: data[0...n-1] data samples */
136
137 long i,j,o,p;
138 float y;
139 float *work=alloca(sizeof(*work)*(m+n));
140
141 if(!prime)
142 for(i=0;i<m;i++)
143 work[i]=0.f;
144 else
145 for(i=0;i<m;i++)
146 work[i]=prime[i];
147
148 for(i=0;i<n;i++){
149 y=0;
150 o=i;
151 p=m;
152 for(j=0;j<m;j++)
153 y-=work[o++]*coeff[--p];
154
155 data[i]=work[o]=y;
156 }
157 }
158
159 #include "dumb.h"
160 #include "internal/dumb.h"
161 #include "internal/it.h"
162
163 enum { lpc_max = 256 }; /* Maximum number of input samples to train the function */
164 enum { lpc_order = 32 }; /* Order of the filter */
165 enum { lpc_extra = 64 }; /* How many samples of padding to predict or silence */
166
167
168 /* This extra sample padding is really only needed by the FIR resampler, but it helps the other resamplers as well. */
169
dumb_it_add_lpc(struct DUMB_IT_SIGDATA * sigdata)170 void dumb_it_add_lpc(struct DUMB_IT_SIGDATA *sigdata){
171 float lpc[lpc_order * 2];
172 float lpc_input[lpc_max * 2];
173 float lpc_output[lpc_extra * 2];
174
175 signed char * s8;
176 signed short * s16;
177
178 int n, o, offset, lpc_samples;
179
180 for ( n = 0; n < sigdata->n_samples; n++ ) {
181 IT_SAMPLE * sample = sigdata->sample + n;
182 if ( ( sample->flags & ( IT_SAMPLE_EXISTS | IT_SAMPLE_LOOP) ) == IT_SAMPLE_EXISTS ) {
183 /* If we have enough sample data to train the filter, use the filter to generate the padding */
184 if ( sample->length >= lpc_order ) {
185 lpc_samples = sample->length;
186 if (lpc_samples > lpc_max) lpc_samples = lpc_max;
187 offset = sample->length - lpc_samples;
188
189 if ( sample->flags & IT_SAMPLE_STEREO )
190 {
191 if ( sample->flags & IT_SAMPLE_16BIT )
192 {
193 s16 = ( signed short * ) sample->data;
194 s16 += offset * 2;
195 for ( o = 0; o < lpc_samples; o++ )
196 {
197 lpc_input[ o ] = s16[ o * 2 + 0 ];
198 lpc_input[ o + lpc_max ] = s16[ o * 2 + 1 ];
199 }
200 }
201 else
202 {
203 s8 = ( signed char * ) sample->data;
204 s8 += offset * 2;
205 for ( o = 0; o < lpc_samples; o++ )
206 {
207 lpc_input[ o ] = s8[ o * 2 + 0 ];
208 lpc_input[ o + lpc_max ] = s8[ o * 2 + 1 ];
209 }
210 }
211
212 vorbis_lpc_from_data( lpc_input, lpc, lpc_samples, lpc_order );
213 vorbis_lpc_from_data( lpc_input + lpc_max, lpc + lpc_order, lpc_samples, lpc_order );
214
215 vorbis_lpc_predict( lpc, lpc_input + lpc_samples - lpc_order, lpc_order, lpc_output, lpc_extra );
216 vorbis_lpc_predict( lpc + lpc_order, lpc_input + lpc_max + lpc_samples - lpc_order, lpc_order, lpc_output + lpc_extra, lpc_extra );
217
218 if ( sample->flags & IT_SAMPLE_16BIT )
219 {
220 s16 = ( signed short * ) realloc( sample->data, ( sample->length + lpc_extra ) * 2 * sizeof(short) );
221 sample->data = s16;
222
223 s16 += sample->length * 2;
224 sample->length += lpc_extra;
225
226 for ( o = 0; o < lpc_extra; o++ )
227 {
228 s16[ o * 2 + 0 ] = lpc_output[ o ];
229 s16[ o * 2 + 1 ] = lpc_output[ o + lpc_extra ];
230 }
231 }
232 else
233 {
234 s8 = ( signed char * ) realloc( sample->data, ( sample->length + lpc_extra ) * 2 );
235 sample->data = s8;
236
237 s8 += sample->length * 2;
238 sample->length += lpc_extra;
239
240 for ( o = 0; o < lpc_extra; o++ )
241 {
242 s8[ o * 2 + 0 ] = lpc_output[ o ];
243 s8[ o * 2 + 1 ] = lpc_output[ o + lpc_extra ];
244 }
245 }
246 }
247 else
248 {
249 if ( sample->flags & IT_SAMPLE_16BIT )
250 {
251 s16 = ( signed short * ) sample->data;
252 s16 += offset;
253 for ( o = 0; o < lpc_samples; o++ )
254 {
255 lpc_input[ o ] = s16[ o ];
256 }
257 }
258 else
259 {
260 s8 = ( signed char * ) sample->data;
261 s8 += offset;
262 for ( o = 0; o < lpc_samples; o++ )
263 {
264 lpc_input[ o ] = s8[ o ];
265 }
266 }
267
268 vorbis_lpc_from_data( lpc_input, lpc, lpc_samples, lpc_order );
269
270 vorbis_lpc_predict( lpc, lpc_input + lpc_samples - lpc_order, lpc_order, lpc_output, lpc_extra );
271
272 if ( sample->flags & IT_SAMPLE_16BIT )
273 {
274 s16 = ( signed short * ) realloc( sample->data, ( sample->length + lpc_extra ) * sizeof(short) );
275 sample->data = s16;
276
277 s16 += sample->length;
278 sample->length += lpc_extra;
279
280 for ( o = 0; o < lpc_extra; o++ )
281 {
282 s16[ o ] = lpc_output[ o ];
283 }
284 }
285 else
286 {
287 s8 = ( signed char * ) realloc( sample->data, sample->length + lpc_extra );
288 sample->data = s8;
289
290 s8 += sample->length;
291 sample->length += lpc_extra;
292
293 for ( o = 0; o < lpc_extra; o++ )
294 {
295 s8[ o ] = lpc_output[ o ];
296 }
297 }
298 }
299 }
300 else
301 /* Otherwise, pad with silence. */
302 {
303 offset = sample->length;
304 lpc_samples = lpc_extra;
305
306 sample->length += lpc_samples;
307
308 n = 1;
309 if ( sample->flags & IT_SAMPLE_STEREO ) n *= 2;
310 if ( sample->flags & IT_SAMPLE_16BIT ) n *= 2;
311
312 offset *= n;
313 lpc_samples *= n;
314
315 sample->data = realloc( sample->data, offset + lpc_samples );
316 memset( (char*)sample->data + offset, 0, lpc_samples );
317 }
318 }
319 }
320 }
321