1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
#include "LPC.h"
// http://kbs.cs.tu-berlin.de/~jutta/gsm/lpc.html
#define P_MAX 8 /* order p of LPC analysis, typically 8..14 */
/* Invented by N. Levinson in 1947, modified by J. Durbin in 1959.
*/
double /* returns minimum mean square error */
levinson_durbin(
double const * ac, /* in: [0...p] autocorrelation values */
double * ref, /* out: [0...p-1] reflection coef's */
double * lpc) /* [0...p-1] LPC coefficients */
{
int i, j; double r, error = ac[0];
if (ac[0] == 0) {
for (i = 0; i < P_MAX; i++) {ref[i] = 0;} return 0; }
for (i = 0; i < P_MAX; i++) {
/* Sum up this iteration's reflection coefficient.
*/
r = -ac[i + 1];
for (j = 0; j < i; j++) r -= lpc[j] * ac[i - j];
ref[i] = r /= error;
/* Update LPC coefficients and total error.
*/
lpc[i] = r;
for (j = 0; j < i/2; j++) {
double tmp = lpc[j];
lpc[j] += r * lpc[i-1-j];
lpc[i-1-j] += r * tmp;
}
if (i % 2) lpc[j] += lpc[j] * r;
error *= 1.0 - r * r;
}
return error;
}
/* I. Sch"ur's recursion from 1917 is related to the Levinson-Durbin
* method, but faster on parallel architectures; where Levinson-Durbin
* would take time proportional to p * log(p), Sch"ur only requires
* time proportional to p. The GSM coder uses an integer version of
* the Sch"ur recursion.
*/
double /* returns the minimum mean square error */
schur( double const * ac, /* in: [0...p] autocorrelation c's */
double * ref) /* out: [0...p-1] reflection c's */
{
int i, m; double r, error = ac[0], G[2][P_MAX];
if (ac[0] == 0.0) {
for (i = 0; i < P_MAX; i++) ref[i] = 0;
return 0;
}
/* Initialize rows of the generator matrix G to ac[1...P]
*/
for (i = 0; i < P_MAX; i++) G[0][i] = G[1][i] = ac[i+1];
for (i = 0;;) {
/* Calculate this iteration's reflection
* coefficient and error.
*/
ref[i] = r = -G[1][0] / error;
error += G[1][0] * r;
if (++i >= P_MAX) return error;
/* Update the generator matrix.
*
* Unlike the Levinson-Durbin summing of reflection
* coefficients, this loop could be distributed to
* p processors who each take only constant time.
*/
for (m = 0; m < P_MAX - i; m++) {
G[1][m] = G[1][m+1] + r * G[0][m];
G[0][m] = G[1][m+1] * r + G[0][m];
}
}
}
|
