4d48b091 |
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 | } |