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 https://xiph.org/ *
|
---|
10 | * *
|
---|
11 | ********************************************************************
|
---|
12 |
|
---|
13 | function: LPC low level routines
|
---|
14 |
|
---|
15 | ********************************************************************/
|
---|
16 |
|
---|
17 | /* Some of these routines (autocorrelator, LPC coefficient estimator)
|
---|
18 | are derived from code written by Jutta Degener and Carsten Bormann;
|
---|
19 | thus we include their copyright below. The entirety of this file
|
---|
20 | is freely redistributable on the condition that both of these
|
---|
21 | copyright notices are preserved without modification. */
|
---|
22 |
|
---|
23 | /* Preserved Copyright: *********************************************/
|
---|
24 |
|
---|
25 | /* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
|
---|
26 | Technische Universita"t Berlin
|
---|
27 |
|
---|
28 | Any use of this software is permitted provided that this notice is not
|
---|
29 | removed and that neither the authors nor the Technische Universita"t
|
---|
30 | Berlin are deemed to have made any representations as to the
|
---|
31 | suitability of this software for any purpose nor are held responsible
|
---|
32 | for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
|
---|
33 | THIS SOFTWARE.
|
---|
34 |
|
---|
35 | As a matter of courtesy, the authors request to be informed about uses
|
---|
36 | this software has found, about bugs in this software, and about any
|
---|
37 | improvements that may be of general interest.
|
---|
38 |
|
---|
39 | Berlin, 28.11.1994
|
---|
40 | Jutta Degener
|
---|
41 | Carsten Bormann
|
---|
42 |
|
---|
43 | *********************************************************************/
|
---|
44 |
|
---|
45 | #include <stdlib.h>
|
---|
46 | #include <string.h>
|
---|
47 | #include <math.h>
|
---|
48 | #include "os.h"
|
---|
49 | #include "smallft.h"
|
---|
50 | #include "lpc.h"
|
---|
51 | #include "scales.h"
|
---|
52 | #include "misc.h"
|
---|
53 |
|
---|
54 | /* Autocorrelation LPC coeff generation algorithm invented by
|
---|
55 | N. Levinson in 1947, modified by J. Durbin in 1959. */
|
---|
56 |
|
---|
57 | /* Input : n elements of time doamin data
|
---|
58 | Output: m lpc coefficients, excitation energy */
|
---|
59 |
|
---|
60 | float vorbis_lpc_from_data(float *data,float *lpci,int n,int m){
|
---|
61 | double *aut=alloca(sizeof(*aut)*(m+1));
|
---|
62 | double *lpc=alloca(sizeof(*lpc)*(m));
|
---|
63 | double error;
|
---|
64 | double epsilon;
|
---|
65 | int i,j;
|
---|
66 |
|
---|
67 | /* autocorrelation, p+1 lag coefficients */
|
---|
68 | j=m+1;
|
---|
69 | while(j--){
|
---|
70 | double d=0; /* double needed for accumulator depth */
|
---|
71 | for(i=j;i<n;i++)d+=(double)data[i]*data[i-j];
|
---|
72 | aut[j]=d;
|
---|
73 | }
|
---|
74 |
|
---|
75 | /* Generate lpc coefficients from autocorr values */
|
---|
76 |
|
---|
77 | /* set our noise floor to about -100dB */
|
---|
78 | error=aut[0] * (1. + 1e-10);
|
---|
79 | epsilon=1e-9*aut[0]+1e-10;
|
---|
80 |
|
---|
81 | for(i=0;i<m;i++){
|
---|
82 | double r= -aut[i+1];
|
---|
83 |
|
---|
84 | if(error<epsilon){
|
---|
85 | memset(lpc+i,0,(m-i)*sizeof(*lpc));
|
---|
86 | goto done;
|
---|
87 | }
|
---|
88 |
|
---|
89 | /* Sum up this iteration's reflection coefficient; note that in
|
---|
90 | Vorbis we don't save it. If anyone wants to recycle this code
|
---|
91 | and needs reflection coefficients, save the results of 'r' from
|
---|
92 | each iteration. */
|
---|
93 |
|
---|
94 | for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
|
---|
95 | r/=error;
|
---|
96 |
|
---|
97 | /* Update LPC coefficients and total error */
|
---|
98 |
|
---|
99 | lpc[i]=r;
|
---|
100 | for(j=0;j<i/2;j++){
|
---|
101 | double tmp=lpc[j];
|
---|
102 |
|
---|
103 | lpc[j]+=r*lpc[i-1-j];
|
---|
104 | lpc[i-1-j]+=r*tmp;
|
---|
105 | }
|
---|
106 | if(i&1)lpc[j]+=lpc[j]*r;
|
---|
107 |
|
---|
108 | error*=1.-r*r;
|
---|
109 |
|
---|
110 | }
|
---|
111 |
|
---|
112 | done:
|
---|
113 |
|
---|
114 | /* slightly damp the filter */
|
---|
115 | {
|
---|
116 | double g = .99;
|
---|
117 | double damp = g;
|
---|
118 | for(j=0;j<m;j++){
|
---|
119 | lpc[j]*=damp;
|
---|
120 | damp*=g;
|
---|
121 | }
|
---|
122 | }
|
---|
123 |
|
---|
124 | for(j=0;j<m;j++)lpci[j]=(float)lpc[j];
|
---|
125 |
|
---|
126 | /* we need the error value to know how big an impulse to hit the
|
---|
127 | filter with later */
|
---|
128 |
|
---|
129 | return error;
|
---|
130 | }
|
---|
131 |
|
---|
132 | void vorbis_lpc_predict(float *coeff,float *prime,int m,
|
---|
133 | float *data,long n){
|
---|
134 |
|
---|
135 | /* in: coeff[0...m-1] LPC coefficients
|
---|
136 | prime[0...m-1] initial values (allocated size of n+m-1)
|
---|
137 | out: data[0...n-1] data samples */
|
---|
138 |
|
---|
139 | long i,j,o,p;
|
---|
140 | float y;
|
---|
141 | float *work=alloca(sizeof(*work)*(m+n));
|
---|
142 |
|
---|
143 | if(!prime)
|
---|
144 | for(i=0;i<m;i++)
|
---|
145 | work[i]=0.f;
|
---|
146 | else
|
---|
147 | for(i=0;i<m;i++)
|
---|
148 | work[i]=prime[i];
|
---|
149 |
|
---|
150 | for(i=0;i<n;i++){
|
---|
151 | y=0;
|
---|
152 | o=i;
|
---|
153 | p=m;
|
---|
154 | for(j=0;j<m;j++)
|
---|
155 | y-=work[o++]*coeff[--p];
|
---|
156 |
|
---|
157 | data[i]=work[o]=y;
|
---|
158 | }
|
---|
159 | }
|
---|