]>
Commit | Line | Data |
---|---|---|
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 | #ifdef HAVE_CONFIG_H | |
47 | # include "config.h" | |
48 | #endif | |
49 | ||
50 | #if HAVE_STDINT_H | |
51 | # include <stdint.h> | |
52 | #endif | |
53 | ||
54 | #include <stdlib.h> | |
55 | #include <string.h> | |
56 | #include <math.h> | |
57 | #include "lpc.h" | |
58 | #include "lpcm.h" | |
59 | ||
60 | /* Autocorrelation LPC coeff generation algorithm invented by | |
61 | N. Levinson in 1947, modified by J. Durbin in 1959. */ | |
62 | ||
63 | /* Input : n elements of time doamin data | |
64 | Output: m lpc coefficients, excitation energy */ | |
65 | ||
66 | float vorbis_lpc_from_data(short *data,float *lpci,int n,int m,int stride){ | |
67 | double *aut=malloc(sizeof(*aut)*(m+1)); | |
68 | double *lpc=malloc(sizeof(*lpc)*(m)); | |
69 | double error; | |
70 | double epsilon; | |
71 | int i,j; | |
72 | ||
73 | /* autocorrelation, p+1 lag coefficients */ | |
74 | j=m+1; | |
75 | while(j--){ | |
76 | double d=0; /* double needed for accumulator depth */ | |
77 | for(i=j;i<n;i++)d+=(double)data[i*stride]*data[(i-j)*stride]/1073741824.0; | |
78 | aut[j]=d; | |
79 | } | |
80 | ||
81 | /* Generate lpc coefficients from autocorr values */ | |
82 | ||
83 | /* set our noise floor to about -100dB */ | |
84 | error=aut[0] * (1. + 1e-10); | |
85 | epsilon=1e-9*aut[0]+1e-10; | |
86 | ||
87 | for(i=0;i<m;i++){ | |
88 | double r= -aut[i+1]; | |
89 | ||
90 | if(error<epsilon){ | |
91 | memset(lpc+i,0,(m-i)*sizeof(*lpc)); | |
92 | goto done; | |
93 | } | |
94 | ||
95 | /* Sum up this iteration's reflection coefficient; note that in | |
96 | Vorbis we don't save it. If anyone wants to recycle this code | |
97 | and needs reflection coefficients, save the results of 'r' from | |
98 | each iteration. */ | |
99 | ||
100 | for(j=0;j<i;j++)r-=lpc[j]*aut[i-j]; | |
101 | r/=error; | |
102 | ||
103 | /* Update LPC coefficients and total error */ | |
104 | ||
105 | lpc[i]=r; | |
106 | for(j=0;j<i/2;j++){ | |
107 | double tmp=lpc[j]; | |
108 | ||
109 | lpc[j]+=r*lpc[i-1-j]; | |
110 | lpc[i-1-j]+=r*tmp; | |
111 | } | |
112 | if(i&1)lpc[j]+=lpc[j]*r; | |
113 | ||
114 | error*=1.-r*r; | |
115 | ||
116 | } | |
117 | ||
118 | done: | |
119 | ||
120 | /* slightly damp the filter */ | |
121 | { | |
122 | double g = .99; | |
123 | double damp = g; | |
124 | for(j=0;j<m;j++){ | |
125 | lpc[j]*=damp; | |
126 | damp*=g; | |
127 | } | |
128 | } | |
129 | ||
130 | for(j=0;j<m;j++)lpci[j]=(float)lpc[j]; | |
131 | ||
132 | /* we need the error value to know how big an impulse to hit the | |
133 | filter with later */ | |
134 | ||
135 | free(aut); | |
136 | free(lpc); | |
137 | return error; | |
138 | } | |
139 | ||
140 | void vorbis_lpc_predict(float *coeff,short *prime,int m, | |
141 | short *data,long n,int stride){ | |
142 | ||
143 | /* in: coeff[0...m-1] LPC coefficients | |
144 | prime[0...m-1] initial values (allocated size of n+m-1) | |
145 | out: data[0...n-1] data samples */ | |
146 | ||
147 | long i,j,o,p; | |
148 | float y; | |
149 | float *work=malloc(sizeof(*work)*(m+n)); | |
150 | ||
151 | if(!prime) | |
152 | for(i=0;i<m;i++) | |
153 | work[i]=0.f; | |
154 | else | |
155 | for(i=0;i<m;i++) | |
156 | work[i]=prime[i*stride]/32768.0f; | |
157 | ||
158 | for(i=0;i<n;i++){ | |
159 | y=0; | |
160 | o=i; | |
161 | p=m; | |
162 | for(j=0;j<m;j++) | |
163 | y-=work[o++]*coeff[--p]; | |
164 | ||
165 | work[o]=y; | |
166 | data[i*stride]=lrint(pcm_clip(y*32768.0,-32768.0,32767.0)); | |
167 | } | |
168 | free(work); | |
169 | } |