1 double LPCBurg(const double x[]/*in wave data*/, size_t N, double a[]/*out LPC*/, size_t M)
 2 {
 3     std::vector<double> b1 (N, 0);
 4     std::vector<double> b2 (N, 0);
 5     std::vector<double> aa (M, 0);
 6
 7     // (3)
 8
 9     double gain = 0.0;
10     for(size_t j = 0; j < N; j++)
11     {
12         gain += x[j] * x[j];
13     }
14
15     gain /= N;
16     if(gain <= 0)
17     {
18         return 0.0;    // warning empty
19     }
20
21     // (9)
22
23     b1[0] = x[0];
24     b2[N - 2] = x[N - 1];
25     for(size_t j = 1; j < N - 1; j++)
26     {
27         b1[j] = b2[j - 1] = x[j];
28     }
29
30     for(size_t i = 0; ; ++i)
31     {
32         // (7)
33
34         double num = 0.0, denum = 0.0;
35         for(size_t j = 0; j < N - i - 1; j++)
36         {
37             num += b1[j] * b2[j];
38             denum += b1[j] * b1[j] + b2[j] * b2[j];
39         }
40
41         if(denum <= 0)
42         {
43             return 0.0;    // warning ill-conditioned
44         }
45
46         a[i] = 2.0 * num / denum;
47
48         // (10)
49
50         gain *= 1.0 - a[i] * a[i];
51
52         // (5)
53
54         for(size_t j = 0; j+1 <= i; j++)
55         {
56             a[j] = aa[j] - a[i] * aa[i - j - 1];
57         }
58
59         if(i == M-1) break;
60
61         // (8)  Watch out: i -> i+1
62
63         for(size_t j = 0; j <= i; j++)
64         {
65             aa[j] = a[j];
66         }
67
68         for(size_t j = 0; j <= N - i - 2; j++)
69         {
70             b1[j] -= aa[i] * b2[j];
71             b2[j] = b2[j + 1] - aa[i] * b1[j + 1];
72         }
73     }
74
75     return gain;
76 }