summaryrefslogtreecommitdiff
path: root/lpc.c
diff options
context:
space:
mode:
authorerdgeist <erdgeist@erdgeist.org>2025-08-15 12:42:40 +0200
committererdgeist <erdgeist@erdgeist.org>2025-08-15 12:42:40 +0200
commit30325d24d107dbf133da39f7c96d1510fd1c9449 (patch)
tree932baa5b2a4475821f16dccf9e3e05011daa6d92 /lpc.c
parent9022d768021bbe15c7815cc6f8b64218b46f0e10 (diff)
Bump to codec2 version 1.2.0erdgeist-bump-to-1.2.0
Diffstat (limited to 'lpc.c')
-rw-r--r--lpc.c223
1 files changed, 97 insertions, 126 deletions
diff --git a/lpc.c b/lpc.c
index 31442cd..3117a1a 100644
--- a/lpc.c
+++ b/lpc.c
@@ -25,16 +25,18 @@
25 along with this program; if not, see <http://www.gnu.org/licenses/>. 25 along with this program; if not, see <http://www.gnu.org/licenses/>.
26*/ 26*/
27 27
28#define LPC_MAX_N 512 /* maximum no. of samples in frame */ 28#define LPC_MAX_N 512 /* maximum no. of samples in frame */
29#define PI 3.141592654 /* mathematical constant */ 29#define PI 3.141592654 /* mathematical constant */
30 30
31#define ALPHA 1.0 31#define ALPHA 1.0
32#define BETA 0.94 32#define BETA 0.94
33
34#include "lpc.h"
33 35
34#include <assert.h> 36#include <assert.h>
35#include <math.h> 37#include <math.h>
38
36#include "defines.h" 39#include "defines.h"
37#include "lpc.h"
38 40
39/*---------------------------------------------------------------------------*\ 41/*---------------------------------------------------------------------------*\
40 42
@@ -42,29 +44,25 @@
42 44
43 Pre-emphasise (high pass filter with zero close to 0 Hz) a frame of 45 Pre-emphasise (high pass filter with zero close to 0 Hz) a frame of
44 speech samples. Helps reduce dynamic range of LPC spectrum, giving 46 speech samples. Helps reduce dynamic range of LPC spectrum, giving
45 greater weight and hense a better match to low energy formants. 47 greater weight and hence a better match to low energy formants.
46 48
47 Should be balanced by de-emphasis of the output speech. 49 Should be balanced by de-emphasis of the output speech.
48 50
49\*---------------------------------------------------------------------------*/ 51\*---------------------------------------------------------------------------*/
50 52
51void pre_emp( 53void pre_emp(float Sn_pre[], /* output frame of speech samples */
52 float Sn_pre[], /* output frame of speech samples */ 54 float Sn[], /* input frame of speech samples */
53 float Sn[], /* input frame of speech samples */ 55 float *mem, /* Sn[-1]single sample memory */
54 float *mem, /* Sn[-1]single sample memory */ 56 int Nsam /* number of speech samples to use */
55 int Nsam /* number of speech samples to use */ 57) {
56) 58 int i;
57{
58 int i;
59
60 for(i=0; i<Nsam; i++) {
61 Sn_pre[i] = Sn[i] - ALPHA * mem[0];
62 mem[0] = Sn[i];
63 }
64 59
60 for (i = 0; i < Nsam; i++) {
61 Sn_pre[i] = Sn[i] - ALPHA * mem[0];
62 mem[0] = Sn[i];
63 }
65} 64}
66 65
67
68/*---------------------------------------------------------------------------*\ 66/*---------------------------------------------------------------------------*\
69 67
70 de_emp() 68 de_emp()
@@ -73,23 +71,19 @@ void pre_emp(
73 71
74\*---------------------------------------------------------------------------*/ 72\*---------------------------------------------------------------------------*/
75 73
76void de_emp( 74void de_emp(float Sn_de[], /* output frame of speech samples */
77 float Sn_de[], /* output frame of speech samples */ 75 float Sn[], /* input frame of speech samples */
78 float Sn[], /* input frame of speech samples */ 76 float *mem, /* Sn[-1]single sample memory */
79 float *mem, /* Sn[-1]single sample memory */ 77 int Nsam /* number of speech samples to use */
80 int Nsam /* number of speech samples to use */ 78) {
81) 79 int i;
82{
83 int i;
84
85 for(i=0; i<Nsam; i++) {
86 Sn_de[i] = Sn[i] + BETA * mem[0];
87 mem[0] = Sn_de[i];
88 }
89 80
81 for (i = 0; i < Nsam; i++) {
82 Sn_de[i] = Sn[i] + BETA * mem[0];
83 mem[0] = Sn_de[i];
84 }
90} 85}
91 86
92
93/*---------------------------------------------------------------------------*\ 87/*---------------------------------------------------------------------------*\
94 88
95 hanning_window() 89 hanning_window()
@@ -98,16 +92,14 @@ void de_emp(
98 92
99\*---------------------------------------------------------------------------*/ 93\*---------------------------------------------------------------------------*/
100 94
101void hanning_window( 95void hanning_window(float Sn[], /* input frame of speech samples */
102 float Sn[], /* input frame of speech samples */ 96 float Wn[], /* output frame of windowed samples */
103 float Wn[], /* output frame of windowed samples */ 97 int Nsam /* number of samples */
104 int Nsam /* number of samples */ 98) {
105) 99 int i; /* loop variable */
106{
107 int i; /* loop variable */
108 100
109 for(i=0; i<Nsam; i++) 101 for (i = 0; i < Nsam; i++)
110 Wn[i] = Sn[i]*(0.5 - 0.5*cosf(2*PI*(float)i/(Nsam-1))); 102 Wn[i] = Sn[i] * (0.5 - 0.5 * cosf(2 * PI * (float)i / (Nsam - 1)));
111} 103}
112 104
113/*---------------------------------------------------------------------------*\ 105/*---------------------------------------------------------------------------*\
@@ -119,19 +111,16 @@ void hanning_window(
119 111
120\*---------------------------------------------------------------------------*/ 112\*---------------------------------------------------------------------------*/
121 113
122void autocorrelate( 114void autocorrelate(float Sn[], /* frame of Nsam windowed speech samples */
123 float Sn[], /* frame of Nsam windowed speech samples */ 115 float Rn[], /* array of P+1 autocorrelation coefficients */
124 float Rn[], /* array of P+1 autocorrelation coefficients */ 116 int Nsam, /* number of windowed samples to use */
125 int Nsam, /* number of windowed samples to use */ 117 int order /* order of LPC analysis */
126 int order /* order of LPC analysis */ 118) {
127) 119 int i, j; /* loop variables */
128{
129 int i,j; /* loop variables */
130 120
131 for(j=0; j<order+1; j++) { 121 for (j = 0; j < order + 1; j++) {
132 Rn[j] = 0.0; 122 Rn[j] = 0.0;
133 for(i=0; i<Nsam-j; i++) 123 for (i = 0; i < Nsam - j; i++) Rn[j] += Sn[i] * Sn[i + j];
134 Rn[j] += Sn[i]*Sn[i+j];
135 } 124 }
136} 125}
137 126
@@ -150,36 +139,31 @@ void autocorrelate(
150 139
151\*---------------------------------------------------------------------------*/ 140\*---------------------------------------------------------------------------*/
152 141
153void levinson_durbin( 142void levinson_durbin(float R[], /* order+1 autocorrelation coeff */
154 float R[], /* order+1 autocorrelation coeff */ 143 float lpcs[], /* order+1 LPC's */
155 float lpcs[], /* order+1 LPC's */ 144 int order /* order of the LPC analysis */
156 int order /* order of the LPC analysis */ 145) {
157) 146 float a[order + 1][order + 1];
158{
159 float a[order+1][order+1];
160 float sum, e, k; 147 float sum, e, k;
161 int i,j; /* loop variables */ 148 int i, j; /* loop variables */
162 149
163 e = R[0]; /* Equation 38a, Makhoul */ 150 e = R[0]; /* Equation 38a, Makhoul */
164 151
165 for(i=1; i<=order; i++) { 152 for (i = 1; i <= order; i++) {
166 sum = 0.0; 153 sum = 0.0;
167 for(j=1; j<=i-1; j++) 154 for (j = 1; j <= i - 1; j++) sum += a[i - 1][j] * R[i - j];
168 sum += a[i-1][j]*R[i-j]; 155 k = -1.0 * (R[i] + sum) / e; /* Equation 38b, Makhoul */
169 k = -1.0*(R[i] + sum)/e; /* Equation 38b, Makhoul */ 156 if (fabsf(k) > 1.0) k = 0.0;
170 if (fabsf(k) > 1.0)
171 k = 0.0;
172 157
173 a[i][i] = k; 158 a[i][i] = k;
174 159
175 for(j=1; j<=i-1; j++) 160 for (j = 1; j <= i - 1; j++)
176 a[i][j] = a[i-1][j] + k*a[i-1][i-j]; /* Equation 38c, Makhoul */ 161 a[i][j] = a[i - 1][j] + k * a[i - 1][i - j]; /* Equation 38c, Makhoul */
177 162
178 e *= (1-k*k); /* Equation 38d, Makhoul */ 163 e *= (1 - k * k); /* Equation 38d, Makhoul */
179 } 164 }
180 165
181 for(i=1; i<=order; i++) 166 for (i = 1; i <= order; i++) lpcs[i] = a[order][i];
182 lpcs[i] = a[order][i];
183 lpcs[0] = 1.0; 167 lpcs[0] = 1.0;
184} 168}
185 169
@@ -194,20 +178,17 @@ void levinson_durbin(
194 178
195\*---------------------------------------------------------------------------*/ 179\*---------------------------------------------------------------------------*/
196 180
197void inverse_filter( 181void inverse_filter(float Sn[], /* Nsam input samples */
198 float Sn[], /* Nsam input samples */ 182 float a[], /* LPCs for this frame of samples */
199 float a[], /* LPCs for this frame of samples */ 183 int Nsam, /* number of samples */
200 int Nsam, /* number of samples */ 184 float res[], /* Nsam residual samples */
201 float res[], /* Nsam residual samples */ 185 int order /* order of LPC */
202 int order /* order of LPC */ 186) {
203) 187 int i, j; /* loop variables */
204{ 188
205 int i,j; /* loop variables */ 189 for (i = 0; i < Nsam; i++) {
206
207 for(i=0; i<Nsam; i++) {
208 res[i] = 0.0; 190 res[i] = 0.0;
209 for(j=0; j<=order; j++) 191 for (j = 0; j <= order; j++) res[i] += Sn[i - j] * a[j];
210 res[i] += Sn[i-j]*a[j];
211 } 192 }
212} 193}
213 194
@@ -223,7 +204,7 @@ void inverse_filter(
223 The synthesis filter has memory as well, this is treated in the same way 204 The synthesis filter has memory as well, this is treated in the same way
224 as the memory for the inverse filter (see inverse_filter() notes above). 205 as the memory for the inverse filter (see inverse_filter() notes above).
225 The difference is that the memory for the synthesis filter is stored in 206 The difference is that the memory for the synthesis filter is stored in
226 the output array, wheras the memory of the inverse filter is stored in the 207 the output array, whereas the memory of the inverse filter is stored in the
227 input array. 208 input array.
228 209
229 Note: the calling function must update the filter memory. 210 Note: the calling function must update the filter memory.
@@ -231,21 +212,19 @@ void inverse_filter(
231\*---------------------------------------------------------------------------*/ 212\*---------------------------------------------------------------------------*/
232 213
233void synthesis_filter( 214void synthesis_filter(
234 float res[], /* Nsam input residual (excitation) samples */ 215 float res[], /* Nsam input residual (excitation) samples */
235 float a[], /* LPCs for this frame of speech samples */ 216 float a[], /* LPCs for this frame of speech samples */
236 int Nsam, /* number of speech samples */ 217 int Nsam, /* number of speech samples */
237 int order, /* LPC order */ 218 int order, /* LPC order */
238 float Sn_[] /* Nsam output synthesised speech samples */ 219 float Sn_[] /* Nsam output synthesised speech samples */
239) 220) {
240{ 221 int i, j; /* loop variables */
241 int i,j; /* loop variables */
242 222
243 /* Filter Nsam samples */ 223 /* Filter Nsam samples */
244 224
245 for(i=0; i<Nsam; i++) { 225 for (i = 0; i < Nsam; i++) {
246 Sn_[i] = res[i]*a[0]; 226 Sn_[i] = res[i] * a[0];
247 for(j=1; j<=order; j++) 227 for (j = 1; j <= order; j++) Sn_[i] -= Sn_[i - j] * a[j];
248 Sn_[i] -= Sn_[i-j]*a[j];
249 } 228 }
250} 229}
251 230
@@ -258,29 +237,25 @@ void synthesis_filter(
258 237
259\*---------------------------------------------------------------------------*/ 238\*---------------------------------------------------------------------------*/
260 239
261void find_aks( 240void find_aks(float Sn[], /* Nsam samples with order sample memory */
262 float Sn[], /* Nsam samples with order sample memory */ 241 float a[], /* order+1 LPCs with first coeff 1.0 */
263 float a[], /* order+1 LPCs with first coeff 1.0 */ 242 int Nsam, /* number of input speech samples */
264 int Nsam, /* number of input speech samples */ 243 int order, /* order of the LPC analysis */
265 int order, /* order of the LPC analysis */ 244 float *E /* residual energy */
266 float *E /* residual energy */ 245) {
267) 246 float Wn[LPC_MAX_N]; /* windowed frame of Nsam speech samples */
268{ 247 float R[order + 1]; /* order+1 autocorrelation values of Sn[] */
269 float Wn[LPC_MAX_N]; /* windowed frame of Nsam speech samples */
270 float R[order+1]; /* order+1 autocorrelation values of Sn[] */
271 int i; 248 int i;
272 249
273 assert(Nsam < LPC_MAX_N); 250 assert(Nsam < LPC_MAX_N);
274 251
275 hanning_window(Sn,Wn,Nsam); 252 hanning_window(Sn, Wn, Nsam);
276 autocorrelate(Wn,R,Nsam,order); 253 autocorrelate(Wn, R, Nsam, order);
277 levinson_durbin(R,a,order); 254 levinson_durbin(R, a, order);
278 255
279 *E = 0.0; 256 *E = 0.0;
280 for(i=0; i<=order; i++) 257 for (i = 0; i <= order; i++) *E += a[i] * R[i];
281 *E += a[i]*R[i]; 258 if (*E < 0.0) *E = 1E-12;
282 if (*E < 0.0)
283 *E = 1E-12;
284} 259}
285 260
286/*---------------------------------------------------------------------------*\ 261/*---------------------------------------------------------------------------*\
@@ -291,16 +266,12 @@ void find_aks(
291 266
292\*---------------------------------------------------------------------------*/ 267\*---------------------------------------------------------------------------*/
293 268
294void weight( 269void weight(float ak[], /* vector of order+1 LPCs */
295 float ak[], /* vector of order+1 LPCs */ 270 float gamma, /* weighting factor */
296 float gamma, /* weighting factor */ 271 int order, /* num LPCs (excluding leading 1.0) */
297 int order, /* num LPCs (excluding leading 1.0) */ 272 float akw[] /* weighted vector of order+1 LPCs */
298 float akw[] /* weighted vector of order+1 LPCs */ 273) {
299)
300{
301 int i; 274 int i;
302 275
303 for(i=1; i<=order; i++) 276 for (i = 1; i <= order; i++) akw[i] = ak[i] * powf(gamma, (float)i);
304 akw[i] = ak[i]*powf(gamma,(float)i);
305} 277}
306