1 /*
2 * lpc.c --
3 *
4 * FIXME: This file needs a description here.
5 *
6 * Copyright (c) 1991-2002 The Regents of the University of California.
7 * All rights reserved.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions are met:
11 *
12 * A. Redistributions of source code must retain the above copyright notice,
13 * this list of conditions and the following disclaimer.
14 * B. Redistributions in binary form must reproduce the above copyright notice,
15 * this list of conditions and the following disclaimer in the documentation
16 * and/or other materials provided with the distribution.
17 * C. Neither the names of the copyright holders nor the names of its
18 * contributors may be used to endorse or promote products derived from this
19 * software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS
22 * IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
23 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
24 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 static const char rcsid[] =
35 "@(#) $Header: /usr/mash/src/repository/mash/mash-1/audio/lpc.c,v 1.7 2002/02/03 03:10:46 lim Exp $";
36
37 #include "config.h"
38 #ifdef WIN32
39 #include <winsock.h>
40 #define inline
41 #endif
42 #include <math.h>
43 #include "lpc.h"
44 #include "mulaw.h"
45
46 #define MAXWINDOW 1000 /* Max analysis window length */
47 #define FS 8000.0F /* Sampling rate */
48
49 #define DOWN 5 /* Decimation for pitch analyzer */
50 #define PITCHORDER 4 /* Model order for pitch analyzer */
51 #define FC 600.0F /* Pitch analyzer filter cutoff */
52 #define MINPIT 50.0 /* Minimum pitch */
53 #define MAXPIT 300.0 /* Maximum pitch */
54
55 #define MINPER (int)(FS/(DOWN*MAXPIT)+.5) /* Minimum period */
56 #define MAXPER (int)(FS/(DOWN*MINPIT)+.5) /* Maximum period */
57
58 #define WSCALE 1.5863F /* Energy loss due to windowing */
59
60 #define BUFLEN ((AUDIO_FRAMESIZE * 3) / 2)
61
62 static float s[MAXWINDOW], y[MAXWINDOW], h[MAXWINDOW];
63 static float fa[6], u, u1, yp1, yp2;
64
65 #define BIAS 0x84 /* define the add-in bias for 16 bit samples */
66 #define CLIP 32635
67
68
69 inline static void
70 auto_correl(float *w, int n, int p, float *r)
71 {
72 int i, k, nk;
73
74 for (k = 0; k <= p; k++) {
75 nk = n - k;
76 r[k] = 0.0f;
77 for (i = 0; i < nk; i++)
78 r[k] += w[i] * w[i + k];
79 }
80 }
81
82 inline static void
83 durbin(float *r, int p, float *k, float *g)
84 {
85 int i, j;
86 float a[LPC_FILTORDER + 1], at[LPC_FILTORDER + 1], e;
87
88 for (i = 0; i <= p; i++)
89 a[i] = at[i] = 0.0f;
90
91 e = r[0];
92 for (i = 1; i <= p; i++) {
93 k[i] = -r[i];
94 for (j = 1; j < i; j++) {
95 at[j] = a[j];
96 k[i] -= a[j] * r[i - j];
97 }
98 k[i] /= e;
99 a[i] = k[i];
100 for (j = 1; j < i; j++)
101 a[j] = at[j] + k[i] * at[i - j];
102 e *= 1.0f - k[i] * k[i];
103 }
104
105 *g = (float)sqrt(e);
106 }
107
108 inline static void
109 inverse_filter(float *w, float *k)
110 {
111 int i, j;
112 float b[PITCHORDER + 1], bp[PITCHORDER + 1], f[PITCHORDER + 1];
113
114 for (i = 0; i <= PITCHORDER; i++)
115 b[i] = f[i] = bp[i] = 0.0;
116
117 for (i = 0; i < BUFLEN / DOWN; i++) {
118 f[0] = b[0] = w[i];
119 for (j = 1; j <= PITCHORDER; j++) {
120 f[j] = f[j - 1] + k[j] * bp[j - 1];
121 b[j] = k[j] * f[j - 1] + bp[j - 1];
122 bp[j - 1] = b[j - 1];
123 }
124 w[i] = f[PITCHORDER];
125 }
126 }
127
128 inline static void
129 calc_pitch(float *w, float *per)
130 {
131 int i, j, rpos;
132 float d[MAXWINDOW / DOWN], k[PITCHORDER + 1], r[MAXPER + 1], g, rmax;
133 float rval, rm, rp;
134 float a, b, c, x, y;
135 static int vuv = 0;
136
137 for (i = 0, j = 0; i < BUFLEN; i += DOWN)
138 d[j++] = w[i];
139 auto_correl(d, BUFLEN / DOWN, PITCHORDER, r);
140 durbin(r, PITCHORDER, k, &g);
141 inverse_filter(d, k);
142 auto_correl(d, BUFLEN / DOWN, MAXPER + 1, r);
143 rpos = 0;
144 rmax = 0.0;
145 for (i = MINPER; i <= MAXPER; i++) {
146 if (r[i] > rmax) {
147 rmax = r[i];
148 rpos = i;
149 }
150 }
151
152 rm = r[rpos - 1];
153 rp = r[rpos + 1];
154 rval = rmax / r[0];
155
156 a = 0.5F * rm - rmax + 0.5F * rp;
157 b = -0.5F * rm * (2.0F * rpos + 1.0F) +
158 2.0F * rpos * rmax + 0.5F * rp * (1.0F - 2.0F * rpos);
159 c = 0.5F * rm * (rpos * rpos + rpos) +
160 rmax * (1.0F - rpos * rpos) + 0.5F * rp * (rpos * rpos - rpos);
161
162 x = -b / (2.0F * a);
163 y = a * x * x + b * x + c;
164 x *= DOWN;
165
166 rmax = y;
167 rval = rmax / r[0];
168 if (rval >= 0.4 || (vuv == 3 && rval >= 0.3)) {
169 *per = x;
170 vuv = (vuv & 1) * 2 + 1;
171 } else {
172 *per = 0.0;
173 vuv = (vuv & 1) * 2;
174 }
175 }
176
177 void
178 lpc_init(lpcstate_t* state)
179 {
180 int i;
181 float r, v, w, wcT;
182
183 for (i = 0; i < BUFLEN; i++) {
184 s[i] = 0.0;
185 h[i] = WSCALE * (0.54F - 0.46F *
186 (float) cos(2 * M_PI * i / (BUFLEN - 1.0)));
187 }
188 wcT = 2.0F * (float) M_PI * FC / FS;
189 r = 0.36891079f * wcT;
190 v = 0.18445539f * wcT;
191 w = 0.92307712f * wcT;
192 fa[1] = (float) -exp(-r);
193 fa[2] = 1.0F + fa[1];
194 fa[3] = -2.0F * (float) exp(-v) * (float) cos(w);
195 fa[4] = (float) exp(-2.0 * v);
196 fa[5] = 1.0F + fa[3] + fa[4];
197
198 u1 = 0.0f;
199 yp1 = 0.0f;
200 yp2 = 0.0f;
201
202 state->Oldper = 0.0f;
203 state->OldG = 0.0f;
204 for (i = 0; i <= LPC_FILTORDER; i++) {
205 state->Oldk[i] = 0.0f;
206 state->bp[i] = 0.0f;
207 }
208 state->pitchctr = 0;
209 }
210
211 void
212 lpc_analyze(const unsigned char *buf, lpcparams_t *params)
213 {
214 int i, j;
215 float w[MAXWINDOW], r[LPC_FILTORDER + 1];
216 float per, G, k[LPC_FILTORDER + 1];
217
218 for (i = 0, j = BUFLEN - AUDIO_FRAMESIZE; i < AUDIO_FRAMESIZE; i++, j++) {
219 s[j] = (float)(mulawtolin[*buf++]) / 32768.0F;
220 u = fa[2] * s[j] - fa[1] * u1;
221 y[j] = fa[5] * u1 - fa[3] * yp1 - fa[4] * yp2;
222 u1 = u;
223 yp2 = yp1;
224 yp1 = y[j];
225 }
226
227 calc_pitch(y, &per);
228
229 for (i = 0; i < BUFLEN; i++)
230 w[i] = s[i] * h[i];
231 auto_correl(w, BUFLEN, LPC_FILTORDER, r);
232 durbin(r, LPC_FILTORDER, k, &G);
233
234 params->period = (u_short)(per * 256.f);
235 params->gain = (u_char)(G * 256.f);
236 for (i = 0; i < LPC_FILTORDER; i++)
237 params->k[i] = (char)(k[i + 1] * 128.f);
238
239 memcpy(s, s + AUDIO_FRAMESIZE, (BUFLEN - AUDIO_FRAMESIZE) * sizeof(s[0]));
240 memcpy(y, y + AUDIO_FRAMESIZE, (BUFLEN - AUDIO_FRAMESIZE) * sizeof(y[0]));
241 }
242
243 void
244 lpc_synthesize(unsigned char *buf, lpcparams_t *params, lpcstate_t* state)
245 {
246 int i, j;
247 register double u, f, per, G, NewG, Ginc, Newper, perinc;
248 double k[LPC_FILTORDER + 1], Newk[LPC_FILTORDER + 1],
249 kinc[LPC_FILTORDER + 1];
250
251 per = (double) params->period / 256.;
252 G = (double) params->gain / 256.;
253 k[0] = 0.0;
254 for (i = 0; i < LPC_FILTORDER; i++)
255 k[i + 1] = (double) (params->k[i]) / 128.;
256
257 G /= sqrt(BUFLEN / (per == 0.0? 3.0 : per));
258 Newper = state->Oldper;
259 NewG = state->OldG;
260 for (i = 1; i <= LPC_FILTORDER; i++)
261 Newk[i] = state->Oldk[i];
262
263 if (state->Oldper != 0 && per != 0) {
264 perinc = (per - state->Oldper) / (double)AUDIO_FRAMESIZE;
265 Ginc = (G - state->OldG) / (double)AUDIO_FRAMESIZE;
266 for (i = 1; i <= LPC_FILTORDER; i++)
267 kinc[i] = (k[i] - state->Oldk[i]) / (double)AUDIO_FRAMESIZE;
268 } else {
269 perinc = 0.0;
270 Ginc = 0.0;
271 for (i = 1; i <= LPC_FILTORDER; i++)
272 kinc[i] = 0.0;
273 }
274
275 if (Newper == 0)
276 state->pitchctr = 0;
277
278 for (i = 0; i < AUDIO_FRAMESIZE; i++) {
279 if (Newper == 0) {
280 u = ((double)random() / 2147483648.0) * NewG;
281 } else {
282 if (state->pitchctr == 0) {
283 u = NewG;
284 state->pitchctr = (int) Newper;
285 } else {
286 u = 0.0;
287 state->pitchctr--;
288 }
289 }
290
291 f = u;
292 for (j = LPC_FILTORDER; j >= 1; j--) {
293 register double b = state->bp[j - 1];
294 register double kj = Newk[j];
295 Newk[j] = kj + kinc[j];
296 f -= b * kj;
297 b += f * kj;
298 state->bp[j] = b;
299 }
300 state->bp[0] = f;
301
302 *buf++ = lintomulaw[(int)(f * 32768.0) & 0xffff];
303
304 Newper += perinc;
305 NewG += Ginc;
306 }
307
308 state->Oldper = per;
309 state->OldG = G;
310 for (i = 1; i <= LPC_FILTORDER; i++)
311 state->Oldk[i] = k[i];
312 }
313
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.