summaryrefslogtreecommitdiff
path: root/quantise.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 /quantise.c
parent9022d768021bbe15c7815cc6f8b64218b46f0e10 (diff)
Bump to codec2 version 1.2.0erdgeist-bump-to-1.2.0
Diffstat (limited to 'quantise.c')
-rw-r--r--quantise.c1992
1 files changed, 536 insertions, 1456 deletions
diff --git a/quantise.c b/quantise.c
index 37bf8be..da1d821 100644
--- a/quantise.c
+++ b/quantise.c
@@ -24,26 +24,27 @@
24 24
25*/ 25*/
26 26
27#include "quantise.h"
28
27#include <assert.h> 29#include <assert.h>
28#include <ctype.h> 30#include <ctype.h>
31#include <math.h>
29#include <stdio.h> 32#include <stdio.h>
30#include <stdlib.h> 33#include <stdlib.h>
31#include <string.h> 34#include <string.h>
32#include <math.h>
33 35
36#include "codec2_fft.h"
34#include "defines.h" 37#include "defines.h"
35#include "dump.h" 38#include "dump.h"
36#include "quantise.h"
37#include "lpc.h" 39#include "lpc.h"
38#include "lsp.h" 40#include "lsp.h"
39#include "codec2_fft.h"
40#include "phase.h"
41#include "mbest.h" 41#include "mbest.h"
42#include "phase.h"
42 43
43#undef PROFILE 44#undef PROFILE
44#include "machdep.h" 45#include "machdep.h"
45 46
46#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */ 47#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */
47 48
48/*---------------------------------------------------------------------------*\ 49/*---------------------------------------------------------------------------*\
49 50
@@ -52,7 +53,7 @@
52\*---------------------------------------------------------------------------*/ 53\*---------------------------------------------------------------------------*/
53 54
54float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[], 55float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[],
55 int m_pitch, int order); 56 int m_pitch, int order);
56 57
57/*---------------------------------------------------------------------------*\ 58/*---------------------------------------------------------------------------*\
58 59
@@ -60,46 +61,11 @@ float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[],
60 61
61\*---------------------------------------------------------------------------*/ 62\*---------------------------------------------------------------------------*/
62 63
63int lsp_bits(int i) { 64int lsp_bits(int i) { return lsp_cb[i].log2m; }
64 return lsp_cb[i].log2m;
65}
66 65
67int lspd_bits(int i) { 66int lspd_bits(int i) { return lsp_cbd[i].log2m; }
68 return lsp_cbd[i].log2m;
69}
70 67
71#ifndef CORTEX_M4 68int lsp_pred_vq_bits(int i) { return lsp_cbjmv[i].log2m; }
72int mel_bits(int i) {
73 return mel_cb[i].log2m;
74}
75
76int lspmelvq_cb_bits(int i) {
77 return lspmelvq_cb[i].log2m;
78}
79#endif
80
81#ifdef __EXPERIMENTAL__
82int lspdt_bits(int i) {
83 return lsp_cbdt[i].log2m;
84}
85#endif
86
87int lsp_pred_vq_bits(int i) {
88 return lsp_cbjvm[i].log2m;
89}
90
91/*---------------------------------------------------------------------------*\
92
93 quantise_init
94
95 Loads the entire LSP quantiser comprised of several vector quantisers
96 (codebooks).
97
98\*---------------------------------------------------------------------------*/
99
100void quantise_init()
101{
102}
103 69
104/*---------------------------------------------------------------------------*\ 70/*---------------------------------------------------------------------------*\
105 71
@@ -111,7 +77,7 @@ void quantise_init()
111 77
112\*---------------------------------------------------------------------------*/ 78\*---------------------------------------------------------------------------*/
113 79
114long quantise(const float * cb, float vec[], float w[], int k, int m, float *se) 80long quantise(const float *cb, float vec[], float w[], int k, int m, float *se)
115/* float cb[][K]; current VQ codebook */ 81/* float cb[][K]; current VQ codebook */
116/* float vec[]; vector to quantise */ 82/* float vec[]; vector to quantise */
117/* float w[]; weighting vector */ 83/* float w[]; weighting vector */
@@ -119,363 +85,123 @@ long quantise(const float * cb, float vec[], float w[], int k, int m, float *se)
119/* int m; size of codebook */ 85/* int m; size of codebook */
120/* float *se; accumulated squared error */ 86/* float *se; accumulated squared error */
121{ 87{
122 float e; /* current error */ 88 float e; /* current error */
123 long besti; /* best index so far */ 89 long besti; /* best index so far */
124 float beste; /* best error so far */ 90 float beste; /* best error so far */
125 long j; 91 long j;
126 int i; 92 int i;
127 float diff; 93 float diff;
128
129 besti = 0;
130 beste = 1E32;
131 for(j=0; j<m; j++) {
132 e = 0.0;
133 for(i=0; i<k; i++) {
134 diff = cb[j*k+i]-vec[i];
135 e += (diff*w[i] * diff*w[i]);
136 }
137 if (e < beste) {
138 beste = e;
139 besti = j;
140 }
141 }
142
143 *se += beste;
144
145 return(besti);
146}
147
148
149
150/*---------------------------------------------------------------------------*\
151
152 encode_lspds_scalar()
153
154 Scalar/VQ LSP difference quantiser.
155
156\*---------------------------------------------------------------------------*/
157 94
158void encode_lspds_scalar( 95 besti = 0;
159 int indexes[], 96 beste = 1E32;
160 float lsp[], 97 for (j = 0; j < m; j++) {
161 int order 98 e = 0.0;
162) 99 for (i = 0; i < k; i++) {
163{ 100 diff = cb[j * k + i] - vec[i];
164 int i,k,m; 101 e += (diff * w[i] * diff * w[i]);
165 float lsp_hz[order];
166 float lsp__hz[order];
167 float dlsp[order];
168 float dlsp_[order];
169 float wt[order];
170 const float *cb;
171 float se;
172
173 for(i=0; i<order; i++) {
174 wt[i] = 1.0;
175 } 102 }
176 103 if (e < beste) {
177 /* convert from radians to Hz so we can use human readable 104 beste = e;
178 frequencies */ 105 besti = j;
179
180 for(i=0; i<order; i++)
181 lsp_hz[i] = (4000.0/PI)*lsp[i];
182
183 //printf("\n");
184
185 wt[0] = 1.0;
186 for(i=0; i<order; i++) {
187
188 /* find difference from previous qunatised lsp */
189
190 if (i)
191 dlsp[i] = lsp_hz[i] - lsp__hz[i-1];
192 else
193 dlsp[0] = lsp_hz[0];
194
195 k = lsp_cbd[i].k;
196 m = lsp_cbd[i].m;
197 cb = lsp_cbd[i].cb;
198 indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se);
199 dlsp_[i] = cb[indexes[i]*k];
200
201
202 if (i)
203 lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
204 else
205 lsp__hz[0] = dlsp_[0];
206
207 //printf("%d lsp %3.2f dlsp %3.2f dlsp_ %3.2f lsp_ %3.2f\n", i, lsp_hz[i], dlsp[i], dlsp_[i], lsp__hz[i]);
208 } 106 }
107 }
209 108
210} 109 *se += beste;
211
212
213void decode_lspds_scalar(
214 float lsp_[],
215 int indexes[],
216 int order
217)
218{
219 int i,k;
220 float lsp__hz[order];
221 float dlsp_[order];
222 const float *cb;
223
224 for(i=0; i<order; i++) {
225
226 k = lsp_cbd[i].k;
227 cb = lsp_cbd[i].cb;
228 dlsp_[i] = cb[indexes[i]*k];
229
230 if (i)
231 lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
232 else
233 lsp__hz[0] = dlsp_[0];
234
235 lsp_[i] = (PI/4000.0)*lsp__hz[i];
236
237 //printf("%d dlsp_ %3.2f lsp_ %3.2f\n", i, dlsp_[i], lsp__hz[i]);
238 }
239 110
111 return (besti);
240} 112}
241 113
242#ifdef __EXPERIMENTAL__
243/*---------------------------------------------------------------------------*\ 114/*---------------------------------------------------------------------------*\
244 115
245 lspvq_quantise 116 encode_lspds_scalar()
246 117
247 Vector LSP quantiser. 118 Scalar/VQ LSP difference-in-frequency quantiser.
248 119
249\*---------------------------------------------------------------------------*/ 120\*---------------------------------------------------------------------------*/
250 121
251void lspvq_quantise( 122void encode_lspds_scalar(int indexes[], float lsp[], int order) {
252 float lsp[], 123 int i, k, m;
253 float lsp_[], 124 float lsp_hz[order];
254 int order 125 float lsp__hz[order];
255) 126 float dlsp[order];
256{ 127 float dlsp_[order];
257 int i,k,m,ncb, nlsp; 128 float wt[order];
258 float wt[order], lsp_hz[order]; 129 const float *cb;
259 const float *cb; 130 float se;
260 float se; 131
261 int index; 132 for (i = 0; i < order; i++) {
262 133 wt[i] = 1.0;
263 for(i=0; i<order; i++) { 134 }
264 wt[i] = 1.0;
265 lsp_hz[i] = 4000.0*lsp[i]/PI;
266 }
267
268 /* scalar quantise LSPs 1,2,3,4 */
269
270 /* simple uniform scalar quantisers */
271
272 for(i=0; i<4; i++) {
273 k = lsp_cb[i].k;
274 m = lsp_cb[i].m;
275 cb = lsp_cb[i].cb;
276 index = quantise(cb, &lsp_hz[i], wt, k, m, &se);
277 lsp_[i] = cb[index*k]*PI/4000.0;
278 }
279
280 //#define WGHT
281#ifdef WGHT
282 for(i=4; i<9; i++) {
283 wt[i] = 1.0/(lsp[i]-lsp[i-1]) + 1.0/(lsp[i+1]-lsp[i]);
284 //printf("wt[%d] = %f\n", i, wt[i]);
285 }
286 wt[9] = 1.0/(lsp[i]-lsp[i-1]);
287#endif
288
289 /* VQ LSPs 5,6,7,8,9,10 */
290
291 ncb = 4;
292 nlsp = 4;
293 k = lsp_cbjnd[ncb].k;
294 m = lsp_cbjnd[ncb].m;
295 cb = lsp_cbjnd[ncb].cb;
296 index = quantise(cb, &lsp_hz[nlsp], &wt[nlsp], k, m, &se);
297 for(i=4; i<order; i++) {
298 lsp_[i] = cb[index*k+i-4]*(PI/4000.0);
299 //printf("%4.f (%4.f) ", lsp_hz[i], cb[index*k+i-4]);
300 }
301}
302
303/*---------------------------------------------------------------------------*\
304
305 lspjnd_quantise
306 135
307 Experimental JND LSP quantiser. 136 /* convert from radians to Hz so we can use human readable
137 frequencies */
308 138
309\*---------------------------------------------------------------------------*/ 139 for (i = 0; i < order; i++) lsp_hz[i] = (4000.0 / PI) * lsp[i];
310 140
311void lspjnd_quantise(float lsps[], float lsps_[], int order) 141 wt[0] = 1.0;
312{ 142 for (i = 0; i < order; i++) {
313 int i,k,m; 143 /* find difference from previous quantised lsp */
314 float wt[order], lsps_hz[order];
315 const float *cb;
316 float se = 0.0;
317 int index;
318
319 for(i=0; i<order; i++) {
320 wt[i] = 1.0;
321 }
322
323 /* convert to Hz */
324 144
325 for(i=0; i<order; i++) { 145 if (i)
326 lsps_hz[i] = lsps[i]*(4000.0/PI); 146 dlsp[i] = lsp_hz[i] - lsp__hz[i - 1];
327 lsps_[i] = lsps[i]; 147 else
328 } 148 dlsp[0] = lsp_hz[0];
329 149
330 /* simple uniform scalar quantisers */ 150 k = lsp_cbd[i].k;
151 m = lsp_cbd[i].m;
152 cb = lsp_cbd[i].cb;
153 indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se);
154 dlsp_[i] = cb[indexes[i] * k];
331 155
332 for(i=0; i<4; i++) { 156 if (i)
333 k = lsp_cbjnd[i].k; 157 lsp__hz[i] = lsp__hz[i - 1] + dlsp_[i];
334 m = lsp_cbjnd[i].m; 158 else
335 cb = lsp_cbjnd[i].cb; 159 lsp__hz[0] = dlsp_[0];
336 index = quantise(cb, &lsps_hz[i], wt, k, m, &se); 160 }
337 lsps_[i] = cb[index*k]*(PI/4000.0);
338 }
339
340 /* VQ LSPs 5,6,7,8,9,10 */
341
342 k = lsp_cbjnd[4].k;
343 m = lsp_cbjnd[4].m;
344 cb = lsp_cbjnd[4].cb;
345 index = quantise(cb, &lsps_hz[4], &wt[4], k, m, &se);
346 //printf("k = %d m = %d c[0] %f cb[k] %f\n", k,m,cb[0],cb[k]);
347 //printf("index = %4d: ", index);
348 for(i=4; i<order; i++) {
349 lsps_[i] = cb[index*k+i-4]*(PI/4000.0);
350 //printf("%4.f (%4.f) ", lsps_hz[i], cb[index*k+i-4]);
351 }
352 //printf("\n");
353} 161}
354 162
355void compute_weights(const float *x, float *w, int ndim); 163void decode_lspds_scalar(float lsp_[], int indexes[], int order) {
356 164 int i, k;
357/*---------------------------------------------------------------------------*\ 165 float lsp__hz[order];
358 166 float dlsp_[order];
359 lspdt_quantise 167 const float *cb;
360
361 LSP difference in time quantiser. Split VQ, encoding LSPs 1-4 with
362 one VQ, and LSPs 5-10 with a second. Update of previous lsp memory
363 is done outside of this function to handle dT between 10 or 20ms
364 frames.
365
366 mode action
367 ------------------
368
369 LSPDT_ALL VQ LSPs 1-4 and 5-10
370 LSPDT_LOW Just VQ LSPs 1-4, for LSPs 5-10 just copy previous
371 LSPDT_HIGH Just VQ LSPs 5-10, for LSPs 1-4 just copy previous
372
373\*---------------------------------------------------------------------------*/
374
375void lspdt_quantise(float lsps[], float lsps_[], float lsps__prev[], int mode)
376{
377 int i;
378 float wt[LPC_ORD];
379 float lsps_dt[LPC_ORD];
380#ifdef TRY_LSPDT_VQ
381 int k,m;
382 int index;
383 const float *cb;
384 float se = 0.0;
385#endif // TRY_LSPDT_VQ
386
387 //compute_weights(lsps, wt, LPC_ORD);
388 for(i=0; i<LPC_ORD; i++) {
389 wt[i] = 1.0;
390 }
391
392 //compute_weights(lsps, wt, LPC_ORD );
393 168
394 for(i=0; i<LPC_ORD; i++) { 169 for (i = 0; i < order; i++) {
395 lsps_dt[i] = lsps[i] - lsps__prev[i]; 170 k = lsp_cbd[i].k;
396 lsps_[i] = lsps__prev[i]; 171 cb = lsp_cbd[i].cb;
397 } 172 dlsp_[i] = cb[indexes[i] * k];
398 173
399 //#define TRY_LSPDT_VQ 174 if (i)
400#ifdef TRY_LSPDT_VQ 175 lsp__hz[i] = lsp__hz[i - 1] + dlsp_[i];
401 /* this actually improves speech a bit, but 40ms updates works surprsingly well.... */ 176 else
402 k = lsp_cbdt[0].k; 177 lsp__hz[0] = dlsp_[0];
403 m = lsp_cbdt[0].m;
404 cb = lsp_cbdt[0].cb;
405 index = quantise(cb, lsps_dt, wt, k, m, &se);
406 for(i=0; i<LPC_ORD; i++) {
407 lsps_[i] += cb[index*k + i];
408 }
409#endif
410 178
179 lsp_[i] = (PI / 4000.0) * lsp__hz[i];
180 }
411} 181}
412#endif
413 182
414#define MIN(a,b) ((a)<(b)?(a):(b)) 183#define MIN(a, b) ((a) < (b) ? (a) : (b))
415#define MAX_ENTRIES 16384 184#define MAX_ENTRIES 16384
416 185
417void compute_weights(const float *x, float *w, int ndim) 186void compute_weights(const float *x, float *w, int ndim) {
418{
419 int i; 187 int i;
420 w[0] = MIN(x[0], x[1]-x[0]); 188 w[0] = MIN(x[0], x[1] - x[0]);
421 for (i=1;i<ndim-1;i++) 189 for (i = 1; i < ndim - 1; i++) w[i] = MIN(x[i] - x[i - 1], x[i + 1] - x[i]);
422 w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]); 190 w[ndim - 1] = MIN(x[ndim - 1] - x[ndim - 2], PI - x[ndim - 1]);
423 w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], PI-x[ndim-1]);
424
425 for (i=0;i<ndim;i++)
426 w[i] = 1./(.01+w[i]);
427 //w[0]*=3;
428 //w[1]*=2;
429}
430
431#ifdef __EXPERIMENTAL__
432/* LSP weight calculation ported from m-file function kindly submitted
433 by Anssi, OH3GDD */
434 191
435void compute_weights_anssi_mode2(const float *x, float *w, int ndim) 192 for (i = 0; i < ndim; i++) w[i] = 1. / (.01 + w[i]);
436{
437 int i;
438 float d[LPC_ORD];
439
440 assert(ndim == LPC_ORD);
441
442 for(i=0; i<LPC_ORD; i++)
443 d[i] = 1.0;
444
445 d[0] = x[1];
446 for (i=1; i<LPC_ORD-1; i++)
447 d[i] = x[i+1] - x[i-1];
448 d[LPC_ORD-1] = PI - x[8];
449 for (i=0; i<LPC_ORD; i++) {
450 if (x[i]<((400.0/4000.0)*PI))
451 w[i]=5.0/(0.01+d[i]);
452 else if (x[i]<((700.0/4000.0)*PI))
453 w[i]=4.0/(0.01+d[i]);
454 else if (x[i]<((1200.0/4000.0)*PI))
455 w[i]=3.0/(0.01+d[i]);
456 else if (x[i]<((2000.0/4000.0)*PI))
457 w[i]=2.0/(0.01+d[i]);
458 else
459 w[i]=1.0/(0.01+d[i]);
460
461 w[i]=powf(w[i]+0.3, 0.66);
462 }
463} 193}
464#endif
465 194
466int find_nearest(const float *codebook, int nb_entries, float *x, int ndim) 195int find_nearest(const float *codebook, int nb_entries, float *x, int ndim) {
467{
468 int i, j; 196 int i, j;
469 float min_dist = 1e15; 197 float min_dist = 1e15;
470 int nearest = 0; 198 int nearest = 0;
471 199
472 for (i=0;i<nb_entries;i++) 200 for (i = 0; i < nb_entries; i++) {
473 { 201 float dist = 0;
474 float dist=0; 202 for (j = 0; j < ndim; j++)
475 for (j=0;j<ndim;j++) 203 dist += (x[j] - codebook[i * ndim + j]) * (x[j] - codebook[i * ndim + j]);
476 dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]); 204 if (dist < min_dist) {
477 if (dist<min_dist)
478 {
479 min_dist = dist; 205 min_dist = dist;
480 nearest = i; 206 nearest = i;
481 } 207 }
@@ -483,19 +209,18 @@ int find_nearest(const float *codebook, int nb_entries, float *x, int ndim)
483 return nearest; 209 return nearest;
484} 210}
485 211
486int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim) 212int find_nearest_weighted(const float *codebook, int nb_entries, float *x,
487{ 213 const float *w, int ndim) {
488 int i, j; 214 int i, j;
489 float min_dist = 1e15; 215 float min_dist = 1e15;
490 int nearest = 0; 216 int nearest = 0;
491 217
492 for (i=0;i<nb_entries;i++) 218 for (i = 0; i < nb_entries; i++) {
493 { 219 float dist = 0;
494 float dist=0; 220 for (j = 0; j < ndim; j++)
495 for (j=0;j<ndim;j++) 221 dist += w[j] * (x[j] - codebook[i * ndim + j]) *
496 dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]); 222 (x[j] - codebook[i * ndim + j]);
497 if (dist<min_dist) 223 if (dist < min_dist) {
498 {
499 min_dist = dist; 224 min_dist = dist;
500 nearest = i; 225 nearest = i;
501 } 226 }
@@ -503,212 +228,74 @@ int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const
503 return nearest; 228 return nearest;
504} 229}
505 230
506void lspjvm_quantise(float *x, float *xq, int order) 231void lspjmv_quantise(float *x, float *xq, int order) {
507{
508 int i, n1, n2, n3; 232 int i, n1, n2, n3;
509 float err[order], err2[order], err3[order]; 233 float err[order], err2[order], err3[order];
510 float w[order], w2[order], w3[order]; 234 float w[order], w2[order], w3[order];
511 const float *codebook1 = lsp_cbjvm[0].cb; 235 const float *codebook1 = lsp_cbjmv[0].cb;
512 const float *codebook2 = lsp_cbjvm[1].cb; 236 const float *codebook2 = lsp_cbjmv[1].cb;
513 const float *codebook3 = lsp_cbjvm[2].cb; 237 const float *codebook3 = lsp_cbjmv[2].cb;
514 238
515 w[0] = MIN(x[0], x[1]-x[0]); 239 w[0] = MIN(x[0], x[1] - x[0]);
516 for (i=1;i<order-1;i++) 240 for (i = 1; i < order - 1; i++) w[i] = MIN(x[i] - x[i - 1], x[i + 1] - x[i]);
517 w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]); 241 w[order - 1] = MIN(x[order - 1] - x[order - 2], PI - x[order - 1]);
518 w[order-1] = MIN(x[order-1]-x[order-2], PI-x[order-1]);
519 242
520 compute_weights(x, w, order); 243 compute_weights(x, w, order);
521 244
522 n1 = find_nearest(codebook1, lsp_cbjvm[0].m, x, order); 245 n1 = find_nearest(codebook1, lsp_cbjmv[0].m, x, order);
523 246
524 for (i=0;i<order;i++) 247 for (i = 0; i < order; i++) {
525 { 248 xq[i] = codebook1[order * n1 + i];
526 xq[i] = codebook1[order*n1+i];
527 err[i] = x[i] - xq[i]; 249 err[i] = x[i] - xq[i];
528 } 250 }
529 for (i=0;i<order/2;i++) 251 for (i = 0; i < order / 2; i++) {
530 { 252 err2[i] = err[2 * i];
531 err2[i] = err[2*i]; 253 err3[i] = err[2 * i + 1];
532 err3[i] = err[2*i+1]; 254 w2[i] = w[2 * i];
533 w2[i] = w[2*i]; 255 w3[i] = w[2 * i + 1];
534 w3[i] = w[2*i+1];
535 } 256 }
536 n2 = find_nearest_weighted(codebook2, lsp_cbjvm[1].m, err2, w2, order/2); 257 n2 = find_nearest_weighted(codebook2, lsp_cbjmv[1].m, err2, w2, order / 2);
537 n3 = find_nearest_weighted(codebook3, lsp_cbjvm[2].m, err3, w3, order/2); 258 n3 = find_nearest_weighted(codebook3, lsp_cbjmv[2].m, err3, w3, order / 2);
538 259
539 for (i=0;i<order/2;i++) 260 for (i = 0; i < order / 2; i++) {
540 { 261 xq[2 * i] += codebook2[order * n2 / 2 + i];
541 xq[2*i] += codebook2[order*n2/2+i]; 262 xq[2 * i + 1] += codebook3[order * n3 / 2 + i];
542 xq[2*i+1] += codebook3[order*n3/2+i];
543 } 263 }
544} 264}
545 265
266int check_lsp_order(float lsp[], int order) {
267 int i;
268 float tmp;
269 int swaps = 0;
270
271 for (i = 1; i < order; i++)
272 if (lsp[i] < lsp[i - 1]) {
273 // fprintf(stderr, "swap %d\n",i);
274 swaps++;
275 tmp = lsp[i - 1];
276 lsp[i - 1] = lsp[i] - 0.1;
277 lsp[i] = tmp + 0.1;
278 i = 1; /* start check again, as swap may have caused out of order */
279 }
546 280
547#ifndef CORTEX_M4 281 return swaps;
548/* simple (non mbest) 6th order LSP MEL VQ quantiser. Returns MSE of result */
549
550float lspmelvq_quantise(float *x, float *xq, int order)
551{
552 int i, n1, n2, n3;
553 float err[order];
554 const float *codebook1 = lspmelvq_cb[0].cb;
555 const float *codebook2 = lspmelvq_cb[1].cb;
556 const float *codebook3 = lspmelvq_cb[2].cb;
557 float tmp[order];
558 float mse;
559
560 assert(order == lspmelvq_cb[0].k);
561
562 n1 = find_nearest(codebook1, lspmelvq_cb[0].m, x, order);
563
564 for (i=0; i<order; i++) {
565 tmp[i] = codebook1[order*n1+i];
566 err[i] = x[i] - tmp[i];
567 }
568
569 n2 = find_nearest(codebook2, lspmelvq_cb[1].m, err, order);
570
571 for (i=0; i<order; i++) {
572 tmp[i] += codebook2[order*n2+i];
573 err[i] = x[i] - tmp[i];
574 }
575
576 n3 = find_nearest(codebook3, lspmelvq_cb[2].m, err, order);
577
578 mse = 0.0;
579 for (i=0; i<order; i++) {
580 tmp[i] += codebook3[order*n3+i];
581 err[i] = x[i] - tmp[i];
582 mse += err[i]*err[i];
583 }
584
585 for (i=0; i<order; i++) {
586 xq[i] = tmp[i];
587 }
588
589 return mse;
590}
591
592/* 3 stage VQ LSP quantiser useing mbest search. Design and guidance kindly submitted by Anssi, OH3GDD */
593
594float lspmelvq_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest_entries)
595{
596 int i, j, n1, n2, n3;
597 const float *codebook1 = lspmelvq_cb[0].cb;
598 const float *codebook2 = lspmelvq_cb[1].cb;
599 const float *codebook3 = lspmelvq_cb[2].cb;
600 struct MBEST *mbest_stage1, *mbest_stage2, *mbest_stage3;
601 float target[ndim];
602 float w[ndim];
603 int index[MBEST_STAGES];
604 float mse, tmp;
605
606 for(i=0; i<ndim; i++)
607 w[i] = 1.0;
608
609 mbest_stage1 = mbest_create(mbest_entries);
610 mbest_stage2 = mbest_create(mbest_entries);
611 mbest_stage3 = mbest_create(mbest_entries);
612 for(i=0; i<MBEST_STAGES; i++)
613 index[i] = 0;
614
615 /* Stage 1 */
616
617 mbest_search(codebook1, x, w, ndim, lspmelvq_cb[0].m, mbest_stage1, index);
618 MBEST_PRINT("Stage 1:", mbest_stage1);
619
620 /* Stage 2 */
621
622 for (j=0; j<mbest_entries; j++) {
623 index[1] = n1 = mbest_stage1->list[j].index[0];
624 for(i=0; i<ndim; i++)
625 target[i] = x[i] - codebook1[ndim*n1+i];
626 mbest_search(codebook2, target, w, ndim, lspmelvq_cb[1].m, mbest_stage2, index);
627 }
628 MBEST_PRINT("Stage 2:", mbest_stage2);
629
630 /* Stage 3 */
631
632 for (j=0; j<mbest_entries; j++) {
633 index[2] = n1 = mbest_stage2->list[j].index[1];
634 index[1] = n2 = mbest_stage2->list[j].index[0];
635 for(i=0; i<ndim; i++)
636 target[i] = x[i] - codebook1[ndim*n1+i] - codebook2[ndim*n2+i];
637 mbest_search(codebook3, target, w, ndim, lspmelvq_cb[2].m, mbest_stage3, index);
638 }
639 MBEST_PRINT("Stage 3:", mbest_stage3);
640
641 n1 = mbest_stage3->list[0].index[2];
642 n2 = mbest_stage3->list[0].index[1];
643 n3 = mbest_stage3->list[0].index[0];
644 mse = 0.0;
645 for (i=0;i<ndim;i++) {
646 tmp = codebook1[ndim*n1+i] + codebook2[ndim*n2+i] + codebook3[ndim*n3+i];
647 mse += (x[i]-tmp)*(x[i]-tmp);
648 xq[i] = tmp;
649 }
650
651 mbest_destroy(mbest_stage1);
652 mbest_destroy(mbest_stage2);
653 mbest_destroy(mbest_stage3);
654
655 indexes[0] = n1; indexes[1] = n2; indexes[2] = n3;
656
657 return mse;
658}
659
660
661void lspmelvq_decode(int *indexes, float *xq, int ndim)
662{
663 int i, n1, n2, n3;
664 const float *codebook1 = lspmelvq_cb[0].cb;
665 const float *codebook2 = lspmelvq_cb[1].cb;
666 const float *codebook3 = lspmelvq_cb[2].cb;
667
668 n1 = indexes[0]; n2 = indexes[1]; n3 = indexes[2];
669 for (i=0;i<ndim;i++) {
670 xq[i] = codebook1[ndim*n1+i] + codebook2[ndim*n2+i] + codebook3[ndim*n3+i];
671 }
672}
673#endif
674
675
676int check_lsp_order(float lsp[], int order)
677{
678 int i;
679 float tmp;
680 int swaps = 0;
681
682 for(i=1; i<order; i++)
683 if (lsp[i] < lsp[i-1]) {
684 //fprintf(stderr, "swap %d\n",i);
685 swaps++;
686 tmp = lsp[i-1];
687 lsp[i-1] = lsp[i]-0.1;
688 lsp[i] = tmp+0.1;
689 i = 1; /* start check again, as swap may have caused out of order */
690 }
691
692 return swaps;
693} 282}
694 283
695void force_min_lsp_dist(float lsp[], int order) 284void force_min_lsp_dist(float lsp[], int order) {
696{ 285 int i;
697 int i;
698 286
699 for(i=1; i<order; i++) 287 for (i = 1; i < order; i++)
700 if ((lsp[i]-lsp[i-1]) < 0.01) { 288 if ((lsp[i] - lsp[i - 1]) < 0.01) {
701 lsp[i] += 0.01; 289 lsp[i] += 0.01;
702 } 290 }
703} 291}
704 292
705
706/*---------------------------------------------------------------------------*\ 293/*---------------------------------------------------------------------------*\
707 294
708 lpc_post_filter() 295 lpc_post_filter()
709 296
710 Applies a post filter to the LPC synthesis filter power spectrum 297 Applies a post filter to the LPC synthesis filter power spectrum
711 Pw, which supresses the inter-formant energy. 298 Pw, which suppresses the inter-formant energy.
712 299
713 The algorithm is from p267 (Section 8.6) of "Digital Speech", 300 The algorithm is from p267 (Section 8.6) of "Digital Speech",
714 edited by A.M. Kondoz, 1994 published by Wiley and Sons. Chapter 8 301 edited by A.M. Kondoz, 1994 published by Wiley and Sons. Chapter 8
@@ -722,7 +309,7 @@ void force_min_lsp_dist(float lsp[], int order)
722 it should be possible to implement this more efficiently in the 309 it should be possible to implement this more efficiently in the
723 time domain. Just not sure how to handle relative time delays 310 time domain. Just not sure how to handle relative time delays
724 between the synthesis stage and updating these coeffs. A smaller 311 between the synthesis stage and updating these coeffs. A smaller
725 FFT size might also be accetable to save CPU. 312 FFT size might also be acceptable to save CPU.
726 313
727 TODO: 314 TODO:
728 [ ] sync var names between Octave and C version 315 [ ] sync var names between Octave and C version
@@ -733,104 +320,97 @@ void force_min_lsp_dist(float lsp[], int order)
733\*---------------------------------------------------------------------------*/ 320\*---------------------------------------------------------------------------*/
734 321
735void lpc_post_filter(codec2_fftr_cfg fftr_fwd_cfg, float Pw[], float ak[], 322void lpc_post_filter(codec2_fftr_cfg fftr_fwd_cfg, float Pw[], float ak[],
736 int order, int dump, float beta, float gamma, int bass_boost, float E) 323 int order, int dump, float beta, float gamma,
737{ 324 int bass_boost, float E) {
738 int i; 325 int i;
739 float x[FFT_ENC]; /* input to FFTs */ 326 float x[FFT_ENC]; /* input to FFTs */
740 COMP Ww[FFT_ENC/2+1]; /* weighting spectrum */ 327 COMP Ww[FFT_ENC / 2 + 1]; /* weighting spectrum */
741 float Rw[FFT_ENC/2+1]; /* R = WA */ 328 float Rw[FFT_ENC / 2 + 1]; /* R = WA */
742 float e_before, e_after, gain; 329 float e_before, e_after, gain;
743 float Pfw; 330 float Pfw;
744 float max_Rw, min_Rw; 331 float max_Rw, min_Rw;
745 float coeff; 332 float coeff;
746 PROFILE_VAR(tstart, tfft1, taw, tfft2, tww, tr); 333 PROFILE_VAR(tstart, tfft1, taw, tfft2, tww, tr);
747
748 PROFILE_SAMPLE(tstart);
749
750 /* Determine weighting filter spectrum W(exp(jw)) ---------------*/
751
752 for(i=0; i<FFT_ENC; i++) {
753 x[i] = 0.0;
754 }
755 334
756 x[0] = ak[0]; 335 PROFILE_SAMPLE(tstart);
757 coeff = gamma;
758 for(i=1; i<=order; i++) {
759 x[i] = ak[i] * coeff;
760 coeff *= gamma;
761 }
762 codec2_fftr(fftr_fwd_cfg, x, Ww);
763 336
764 PROFILE_SAMPLE_AND_LOG(tfft2, taw, " fft2"); 337 /* Determine weighting filter spectrum W(exp(jw)) ---------------*/
765 338
766 for(i=0; i<FFT_ENC/2; i++) { 339 for (i = 0; i < FFT_ENC; i++) {
767 Ww[i].real = Ww[i].real*Ww[i].real + Ww[i].imag*Ww[i].imag; 340 x[i] = 0.0;
768 } 341 }
769 342
770 PROFILE_SAMPLE_AND_LOG(tww, tfft2, " Ww"); 343 x[0] = ak[0];
344 coeff = gamma;
345 for (i = 1; i <= order; i++) {
346 x[i] = ak[i] * coeff;
347 coeff *= gamma;
348 }
349 codec2_fftr(fftr_fwd_cfg, x, Ww);
771 350
772 /* Determined combined filter R = WA ---------------------------*/ 351 PROFILE_SAMPLE_AND_LOG(tfft2, taw, " fft2");
773 352
774 max_Rw = 0.0; min_Rw = 1E32; 353 for (i = 0; i < FFT_ENC / 2; i++) {
775 for(i=0; i<FFT_ENC/2; i++) { 354 Ww[i].real = Ww[i].real * Ww[i].real + Ww[i].imag * Ww[i].imag;
776 Rw[i] = sqrtf(Ww[i].real * Pw[i]); 355 }
777 if (Rw[i] > max_Rw)
778 max_Rw = Rw[i];
779 if (Rw[i] < min_Rw)
780 min_Rw = Rw[i];
781 356
782 } 357 PROFILE_SAMPLE_AND_LOG(tww, tfft2, " Ww");
783 358
784 PROFILE_SAMPLE_AND_LOG(tr, tww, " R"); 359 /* Determined combined filter R = WA ---------------------------*/
785 360
786 #ifdef DUMP 361 max_Rw = 0.0;
787 if (dump) 362 min_Rw = 1E32;
788 dump_Rw(Rw); 363 for (i = 0; i < FFT_ENC / 2; i++) {
789 #endif 364 Rw[i] = sqrtf(Ww[i].real * Pw[i]);
365 if (Rw[i] > max_Rw) max_Rw = Rw[i];
366 if (Rw[i] < min_Rw) min_Rw = Rw[i];
367 }
790 368
791 /* create post filter mag spectrum and apply ------------------*/ 369 PROFILE_SAMPLE_AND_LOG(tr, tww, " R");
792 370
793 /* measure energy before post filtering */ 371#ifdef DUMP
372 if (dump) dump_Rw(Rw);
373#endif
794 374
795 e_before = 1E-4; 375 /* create post filter mag spectrum and apply ------------------*/
796 for(i=0; i<FFT_ENC/2; i++)
797 e_before += Pw[i];
798 376
799 /* apply post filter and measure energy */ 377 /* measure energy before post filtering */
800 378
801 #ifdef DUMP 379 e_before = 1E-4;
802 if (dump) 380 for (i = 0; i < FFT_ENC / 2; i++) e_before += Pw[i];
803 dump_Pwb(Pw);
804 #endif
805 381
382 /* apply post filter and measure energy */
806 383
807 e_after = 1E-4; 384#ifdef DUMP
808 for(i=0; i<FFT_ENC/2; i++) { 385 if (dump) dump_Pwb(Pw);
809 Pfw = powf(Rw[i], beta); 386#endif
810 Pw[i] *= Pfw * Pfw;
811 e_after += Pw[i];
812 }
813 gain = e_before/e_after;
814 387
815 /* apply gain factor to normalise energy, and LPC Energy */ 388 e_after = 1E-4;
389 for (i = 0; i < FFT_ENC / 2; i++) {
390 Pfw = powf(Rw[i], beta);
391 Pw[i] *= Pfw * Pfw;
392 e_after += Pw[i];
393 }
394 gain = e_before / e_after;
816 395
817 gain *= E; 396 /* apply gain factor to normalise energy, and LPC Energy */
818 for(i=0; i<FFT_ENC/2; i++) { 397
819 Pw[i] *= gain; 398 gain *= E;
820 } 399 for (i = 0; i < FFT_ENC / 2; i++) {
400 Pw[i] *= gain;
401 }
821 402
822 if (bass_boost) { 403 if (bass_boost) {
823 /* add 3dB to first 1 kHz to account for LP effect of PF */ 404 /* add 3dB to first 1 kHz to account for LP effect of PF */
824 405
825 for(i=0; i<FFT_ENC/8; i++) { 406 for (i = 0; i < FFT_ENC / 8; i++) {
826 Pw[i] *= 1.4*1.4; 407 Pw[i] *= 1.4 * 1.4;
827 }
828 } 408 }
409 }
829 410
830 PROFILE_SAMPLE_AND_LOG2(tr, " filt"); 411 PROFILE_SAMPLE_AND_LOG2(tr, " filt");
831} 412}
832 413
833
834/*---------------------------------------------------------------------------*\ 414/*---------------------------------------------------------------------------*\
835 415
836 aks_to_M2() 416 aks_to_M2()
@@ -841,134 +421,125 @@ void lpc_post_filter(codec2_fftr_cfg fftr_fwd_cfg, float Pw[], float ak[],
841 421
842\*---------------------------------------------------------------------------*/ 422\*---------------------------------------------------------------------------*/
843 423
844void aks_to_M2( 424void aks_to_M2(codec2_fftr_cfg fftr_fwd_cfg, float ak[], /* LPC's */
845 codec2_fftr_cfg fftr_fwd_cfg, 425 int order,
846 float ak[], /* LPC's */ 426 MODEL *model, /* sinusoidal model parameters for this frame */
847 int order, 427 float E, /* energy term */
848 MODEL *model, /* sinusoidal model parameters for this frame */ 428 float *snr, /* signal to noise ratio for this frame in dB */
849 float E, /* energy term */ 429 int dump, /* true to dump sample to dump file */
850 float *snr, /* signal to noise ratio for this frame in dB */ 430 int sim_pf, /* true to simulate a post filter */
851 int dump, /* true to dump sample to dump file */ 431 int pf, /* true to enable actual LPC post filter */
852 int sim_pf, /* true to simulate a post filter */ 432 int bass_boost, /* enable LPC filter 0-1kHz 3dB boost */
853 int pf, /* true to enable actual LPC post filter */ 433 float beta, float gamma, /* LPC post filter parameters */
854 int bass_boost, /* enable LPC filter 0-1kHz 3dB boost */ 434 COMP Aw[] /* output power spectrum */
855 float beta, 435) {
856 float gamma, /* LPC post filter parameters */ 436 int i, m; /* loop variables */
857 COMP Aw[] /* output power spectrum */ 437 int am, bm; /* limits of current band */
858) 438 float r; /* no. rads/bin */
859{ 439 float Em; /* energy in band */
860 int i,m; /* loop variables */ 440 float Am; /* spectral amplitude sample */
861 int am,bm; /* limits of current band */
862 float r; /* no. rads/bin */
863 float Em; /* energy in band */
864 float Am; /* spectral amplitude sample */
865 float signal, noise; 441 float signal, noise;
866 PROFILE_VAR(tstart, tfft, tpw, tpf); 442 PROFILE_VAR(tstart, tfft, tpw, tpf);
867 443
868 PROFILE_SAMPLE(tstart); 444 PROFILE_SAMPLE(tstart);
869 445
870 r = TWO_PI/(FFT_ENC); 446 r = TWO_PI / (FFT_ENC);
871 447
872 /* Determine DFT of A(exp(jw)) --------------------------------------------*/ 448 /* Determine DFT of A(exp(jw)) --------------------------------------------*/
873 { 449 {
874 float a[FFT_ENC]; /* input to FFT for power spectrum */ 450 float a[FFT_ENC]; /* input to FFT for power spectrum */
875 451
876 for(i=0; i<FFT_ENC; i++) { 452 for (i = 0; i < FFT_ENC; i++) {
877 a[i] = 0.0; 453 a[i] = 0.0;
878 } 454 }
879 455
880 for(i=0; i<=order; i++) 456 for (i = 0; i <= order; i++) a[i] = ak[i];
881 a[i] = ak[i]; 457 codec2_fftr(fftr_fwd_cfg, a, Aw);
882 codec2_fftr(fftr_fwd_cfg, a, Aw);
883 } 458 }
884 PROFILE_SAMPLE_AND_LOG(tfft, tstart, " fft"); 459 PROFILE_SAMPLE_AND_LOG(tfft, tstart, " fft");
885 460
886 /* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/ 461 /* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/
887 462
888 float Pw[FFT_ENC/2]; 463 float Pw[FFT_ENC / 2];
889 464
890#ifndef FDV_ARM_MATH 465#ifndef FDV_ARM_MATH
891 for(i=0; i<FFT_ENC/2; i++) { 466 for (i = 0; i < FFT_ENC / 2; i++) {
892 Pw[i] = 1.0/(Aw[i].real*Aw[i].real + Aw[i].imag*Aw[i].imag + 1E-6); 467 Pw[i] = 1.0 / (Aw[i].real * Aw[i].real + Aw[i].imag * Aw[i].imag + 1E-6);
893 } 468 }
894#else 469#else
895 // this difference may seem strange, but the gcc for STM32F4 generates almost 5 times 470 // this difference may seem strange, but the gcc for STM32F4 generates almost
896 // faster code with the two loops: 1120 ms -> 242 ms 471 // 5 times faster code with the two loops: 1120 ms -> 242 ms so please leave
897 // so please leave it as is or improve further 472 // it as is or improve further since this code is called 4 times it results in
898 // since this code is called 4 times it results in almost 4ms gain (21ms -> 17ms per audio frame decode @ 1300 ) 473 // almost 4ms gain (21ms -> 17ms per audio frame decode @ 1300 )
899 474
900 for(i=0; i<FFT_ENC/2; i++) 475 for (i = 0; i < FFT_ENC / 2; i++) {
901 { 476 Pw[i] = Aw[i].real * Aw[i].real + Aw[i].imag * Aw[i].imag + 1E-6;
902 Pw[i] = Aw[i].real * Aw[i].real + Aw[i].imag * Aw[i].imag + 1E-6;
903 } 477 }
904 for(i=0; i<FFT_ENC/2; i++) { 478 for (i = 0; i < FFT_ENC / 2; i++) {
905 Pw[i] = 1.0/(Pw[i]); 479 Pw[i] = 1.0 / (Pw[i]);
906 } 480 }
907#endif 481#endif
908 482
909 PROFILE_SAMPLE_AND_LOG(tpw, tfft, " Pw"); 483 PROFILE_SAMPLE_AND_LOG(tpw, tfft, " Pw");
910 484
911 if (pf) 485 if (pf)
912 lpc_post_filter(fftr_fwd_cfg, Pw, ak, order, dump, beta, gamma, bass_boost, E); 486 lpc_post_filter(fftr_fwd_cfg, Pw, ak, order, dump, beta, gamma, bass_boost,
487 E);
913 else { 488 else {
914 for(i=0; i<FFT_ENC/2; i++) { 489 for (i = 0; i < FFT_ENC / 2; i++) {
915 Pw[i] *= E; 490 Pw[i] *= E;
916 } 491 }
917 } 492 }
918 493
919 PROFILE_SAMPLE_AND_LOG(tpf, tpw, " LPC post filter"); 494 PROFILE_SAMPLE_AND_LOG(tpf, tpw, " LPC post filter");
920 495
921 #ifdef DUMP 496#ifdef DUMP
922 if (dump) 497 if (dump) dump_Pw(Pw);
923 dump_Pw(Pw); 498#endif
924 #endif
925 499
926 /* Determine magnitudes from P(w) ----------------------------------------*/ 500 /* Determine magnitudes from P(w) ----------------------------------------*/
927 501
928 /* when used just by decoder {A} might be all zeroes so init signal 502 /* when used just by decoder {A} might be all zeroes so init signal
929 and noise to prevent log(0) errors */ 503 and noise to prevent log(0) errors */
930 504
931 signal = 1E-30; noise = 1E-32; 505 signal = 1E-30;
932 506 noise = 1E-32;
933 for(m=1; m<=model->L; m++) { 507
934 am = (int)((m - 0.5)*model->Wo/r + 0.5); 508 for (m = 1; m <= model->L; m++) {
935 bm = (int)((m + 0.5)*model->Wo/r + 0.5); 509 am = (int)((m - 0.5) * model->Wo / r + 0.5);
936 510 bm = (int)((m + 0.5) * model->Wo / r + 0.5);
937 // FIXME: With arm_rfft_fast_f32 we have to use this 511
938 // otherwise sometimes a to high bm is calculated 512 // FIXME: With arm_rfft_fast_f32 we have to use this
939 // which causes trouble later in the calculation 513 // otherwise sometimes a to high bm is calculated
940 // chain 514 // which causes trouble later in the calculation
941 // it seems for some reason model->Wo is calculated somewhat too high 515 // chain
942 if (bm>FFT_ENC/2) 516 // it seems for some reason model->Wo is calculated somewhat too high
943 { 517 if (bm > FFT_ENC / 2) {
944 bm = FFT_ENC/2; 518 bm = FFT_ENC / 2;
945 } 519 }
946 Em = 0.0; 520 Em = 0.0;
947 521
948 for(i=am; i<bm; i++) 522 for (i = am; i < bm; i++) Em += Pw[i];
949 Em += Pw[i]; 523 Am = sqrtf(Em);
950 Am = sqrtf(Em); 524
951 525 signal += model->A[m] * model->A[m];
952 signal += model->A[m]*model->A[m]; 526 noise += (model->A[m] - Am) * (model->A[m] - Am);
953 noise += (model->A[m] - Am)*(model->A[m] - Am); 527
954 528 /* This code significantly improves perf of LPC model, in
955 /* This code significantly improves perf of LPC model, in 529 particular when combined with phase0. The LPC spectrum tends
956 particular when combined with phase0. The LPC spectrum tends 530 to track just under the peaks of the spectral envelope, and
957 to track just under the peaks of the spectral envelope, and 531 just above nulls. This algorithm does the reverse to
958 just above nulls. This algorithm does the reverse to 532 compensate - raising the amplitudes of spectral peaks, while
959 compensate - raising the amplitudes of spectral peaks, while 533 attenuating the null. This enhances the formants, and
960 attenuating the null. This enhances the formants, and 534 suppresses the energy between formants. */
961 supresses the energy between formants. */ 535
962 536 if (sim_pf) {
963 if (sim_pf) { 537 if (Am > model->A[m]) Am *= 0.7;
964 if (Am > model->A[m]) 538 if (Am < model->A[m]) Am *= 1.4;
965 Am *= 0.7; 539 }
966 if (Am < model->A[m]) 540 model->A[m] = Am;
967 Am *= 1.4;
968 }
969 model->A[m] = Am;
970 } 541 }
971 *snr = 10.0*log10f(signal/noise); 542 *snr = 10.0 * log10f(signal / noise);
972 543
973 PROFILE_SAMPLE_AND_LOG2(tpf, " rec"); 544 PROFILE_SAMPLE_AND_LOG2(tpf, " rec");
974} 545}
@@ -983,19 +554,18 @@ void aks_to_M2(
983 554
984\*---------------------------------------------------------------------------*/ 555\*---------------------------------------------------------------------------*/
985 556
986int encode_Wo(C2CONST *c2const, float Wo, int bits) 557int encode_Wo(C2CONST *c2const, float Wo, int bits) {
987{ 558 int index, Wo_levels = 1 << bits;
988 int index, Wo_levels = 1<<bits; 559 float Wo_min = c2const->Wo_min;
989 float Wo_min = c2const->Wo_min; 560 float Wo_max = c2const->Wo_max;
990 float Wo_max = c2const->Wo_max; 561 float norm;
991 float norm;
992 562
993 norm = (Wo - Wo_min)/(Wo_max - Wo_min); 563 norm = (Wo - Wo_min) / (Wo_max - Wo_min);
994 index = floorf(Wo_levels * norm + 0.5); 564 index = floorf(Wo_levels * norm + 0.5);
995 if (index < 0 ) index = 0; 565 if (index < 0) index = 0;
996 if (index > (Wo_levels-1)) index = Wo_levels-1; 566 if (index > (Wo_levels - 1)) index = Wo_levels - 1;
997 567
998 return index; 568 return index;
999} 569}
1000 570
1001/*---------------------------------------------------------------------------*\ 571/*---------------------------------------------------------------------------*\
@@ -1008,18 +578,17 @@ int encode_Wo(C2CONST *c2const, float Wo, int bits)
1008 578
1009\*---------------------------------------------------------------------------*/ 579\*---------------------------------------------------------------------------*/
1010 580
1011float decode_Wo(C2CONST *c2const, int index, int bits) 581float decode_Wo(C2CONST *c2const, int index, int bits) {
1012{ 582 float Wo_min = c2const->Wo_min;
1013 float Wo_min = c2const->Wo_min; 583 float Wo_max = c2const->Wo_max;
1014 float Wo_max = c2const->Wo_max; 584 float step;
1015 float step; 585 float Wo;
1016 float Wo; 586 int Wo_levels = 1 << bits;
1017 int Wo_levels = 1<<bits;
1018 587
1019 step = (Wo_max - Wo_min)/Wo_levels; 588 step = (Wo_max - Wo_min) / Wo_levels;
1020 Wo = Wo_min + step*(index); 589 Wo = Wo_min + step * (index);
1021 590
1022 return Wo; 591 return Wo;
1023} 592}
1024 593
1025/*---------------------------------------------------------------------------*\ 594/*---------------------------------------------------------------------------*\
@@ -1032,19 +601,18 @@ float decode_Wo(C2CONST *c2const, int index, int bits)
1032 601
1033\*---------------------------------------------------------------------------*/ 602\*---------------------------------------------------------------------------*/
1034 603
1035int encode_log_Wo(C2CONST *c2const, float Wo, int bits) 604int encode_log_Wo(C2CONST *c2const, float Wo, int bits) {
1036{ 605 int index, Wo_levels = 1 << bits;
1037 int index, Wo_levels = 1<<bits; 606 float Wo_min = c2const->Wo_min;
1038 float Wo_min = c2const->Wo_min; 607 float Wo_max = c2const->Wo_max;
1039 float Wo_max = c2const->Wo_max; 608 float norm;
1040 float norm;
1041 609
1042 norm = (log10f(Wo) - log10f(Wo_min))/(log10f(Wo_max) - log10f(Wo_min)); 610 norm = (log10f(Wo) - log10f(Wo_min)) / (log10f(Wo_max) - log10f(Wo_min));
1043 index = floorf(Wo_levels * norm + 0.5); 611 index = floorf(Wo_levels * norm + 0.5);
1044 if (index < 0 ) index = 0; 612 if (index < 0) index = 0;
1045 if (index > (Wo_levels-1)) index = Wo_levels-1; 613 if (index > (Wo_levels - 1)) index = Wo_levels - 1;
1046 614
1047 return index; 615 return index;
1048} 616}
1049 617
1050/*---------------------------------------------------------------------------*\ 618/*---------------------------------------------------------------------------*\
@@ -1057,99 +625,18 @@ int encode_log_Wo(C2CONST *c2const, float Wo, int bits)
1057 625
1058\*---------------------------------------------------------------------------*/ 626\*---------------------------------------------------------------------------*/
1059 627
1060float decode_log_Wo(C2CONST *c2const, int index, int bits) 628float decode_log_Wo(C2CONST *c2const, int index, int bits) {
1061{ 629 float Wo_min = c2const->Wo_min;
1062 float Wo_min = c2const->Wo_min; 630 float Wo_max = c2const->Wo_max;
1063 float Wo_max = c2const->Wo_max; 631 float step;
1064 float step; 632 float Wo;
1065 float Wo; 633 int Wo_levels = 1 << bits;
1066 int Wo_levels = 1<<bits;
1067
1068 step = (log10f(Wo_max) - log10f(Wo_min))/Wo_levels;
1069 Wo = log10f(Wo_min) + step*(index);
1070
1071 return POW10F(Wo);
1072}
1073
1074#if 0
1075/*---------------------------------------------------------------------------*\
1076
1077 FUNCTION....: encode_Wo_dt()
1078 AUTHOR......: David Rowe
1079 DATE CREATED: 6 Nov 2011
1080
1081 Encodes Wo difference from last frame.
1082
1083\*---------------------------------------------------------------------------*/
1084
1085int encode_Wo_dt(C2CONST *c2const, float Wo, float prev_Wo)
1086{
1087 int index, mask, max_index, min_index;
1088 float Wo_min = c2const->Wo_min;
1089 float Wo_max = c2const->Wo_max;
1090 float norm;
1091
1092 norm = (Wo - prev_Wo)/(Wo_max - Wo_min);
1093 index = floorf(WO_LEVELS * norm + 0.5);
1094 //printf("ENC index: %d ", index);
1095
1096 /* hard limit */
1097
1098 max_index = (1 << (WO_DT_BITS-1)) - 1;
1099 min_index = - (max_index+1);
1100 if (index > max_index) index = max_index;
1101 if (index < min_index) index = min_index;
1102 //printf("max_index: %d min_index: %d hard index: %d ",
1103 // max_index, min_index, index);
1104
1105 /* mask so that only LSB WO_DT_BITS remain, bit WO_DT_BITS is the sign bit */
1106
1107 mask = ((1 << WO_DT_BITS) - 1);
1108 index &= mask;
1109 //printf("mask: 0x%x index: 0x%x\n", mask, index);
1110
1111 return index;
1112}
1113
1114/*---------------------------------------------------------------------------*\
1115
1116 FUNCTION....: decode_Wo_dt()
1117 AUTHOR......: David Rowe
1118 DATE CREATED: 6 Nov 2011
1119
1120 Decodes Wo using WO_DT_BITS difference from last frame.
1121
1122\*---------------------------------------------------------------------------*/
1123
1124float decode_Wo_dt(C2CONST *c2const, int index, float prev_Wo)
1125{
1126 float Wo_min = c2const->Wo_min;
1127 float Wo_max = c2const->Wo_max;
1128 float step;
1129 float Wo;
1130 int mask;
1131
1132 /* sign extend index */
1133
1134 //printf("DEC index: %d ");
1135 if (index & (1 << (WO_DT_BITS-1))) {
1136 mask = ~((1 << WO_DT_BITS) - 1);
1137 index |= mask;
1138 }
1139 //printf("DEC mask: 0x%x index: %d \n", mask, index);
1140
1141 step = (Wo_max - Wo_min)/WO_LEVELS;
1142 Wo = prev_Wo + step*(index);
1143
1144 /* bit errors can make us go out of range leading to all sorts of
1145 probs like seg faults */
1146 634
1147 if (Wo > Wo_max) Wo = Wo_max; 635 step = (log10f(Wo_max) - log10f(Wo_min)) / Wo_levels;
1148 if (Wo < Wo_min) Wo = Wo_min; 636 Wo = log10f(Wo_min) + step * (index);
1149 637
1150 return Wo; 638 return POW10F(Wo);
1151} 639}
1152#endif
1153 640
1154/*---------------------------------------------------------------------------*\ 641/*---------------------------------------------------------------------------*\
1155 642
@@ -1163,56 +650,46 @@ float decode_Wo_dt(C2CONST *c2const, int index, float prev_Wo)
1163 650
1164\*---------------------------------------------------------------------------*/ 651\*---------------------------------------------------------------------------*/
1165 652
1166float speech_to_uq_lsps(float lsp[], 653float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[],
1167 float ak[], 654 int m_pitch, int order) {
1168 float Sn[], 655 int i, roots;
1169 float w[], 656 float Wn[m_pitch];
1170 int m_pitch, 657 float R[order + 1];
1171 int order 658 float e, E;
1172) 659
1173{ 660 e = 0.0;
1174 int i, roots; 661 for (i = 0; i < m_pitch; i++) {
1175 float Wn[m_pitch]; 662 Wn[i] = Sn[i] * w[i];
1176 float R[order+1]; 663 e += Wn[i] * Wn[i];
1177 float e, E; 664 }
1178
1179 e = 0.0;
1180 for(i=0; i<m_pitch; i++) {
1181 Wn[i] = Sn[i]*w[i];
1182 e += Wn[i]*Wn[i];
1183 }
1184 665
1185 /* trap 0 energy case as LPC analysis will fail */ 666 /* trap 0 energy case as LPC analysis will fail */
1186 667
1187 if (e == 0.0) { 668 if (e == 0.0) {
1188 for(i=0; i<order; i++) 669 for (i = 0; i < order; i++) lsp[i] = (PI / order) * (float)i;
1189 lsp[i] = (PI/order)*(float)i; 670 return 0.0;
1190 return 0.0; 671 }
1191 }
1192 672
1193 autocorrelate(Wn, R, m_pitch, order); 673 autocorrelate(Wn, R, m_pitch, order);
1194 levinson_durbin(R, ak, order); 674 levinson_durbin(R, ak, order);
1195 675
1196 E = 0.0; 676 E = 0.0;
1197 for(i=0; i<=order; i++) 677 for (i = 0; i <= order; i++) E += ak[i] * R[i];
1198 E += ak[i]*R[i];
1199 678
1200 /* 15 Hz BW expansion as I can't hear the difference and it may help 679 /* 15 Hz BW expansion as I can't hear the difference and it may help
1201 help occasional fails in the LSP root finding. Important to do this 680 help occasional fails in the LSP root finding. Important to do this
1202 after energy calculation to avoid -ve energy values. 681 after energy calculation to avoid -ve energy values.
1203 */ 682 */
1204 683
1205 for(i=0; i<=order; i++) 684 for (i = 0; i <= order; i++) ak[i] *= powf(0.994, (float)i);
1206 ak[i] *= powf(0.994,(float)i);
1207 685
1208 roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1); 686 roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1);
1209 if (roots != order) { 687 if (roots != order) {
1210 /* if root finding fails use some benign LSP values instead */ 688 /* if root finding fails use some benign LSP values instead */
1211 for(i=0; i<order; i++) 689 for (i = 0; i < order; i++) lsp[i] = (PI / order) * (float)i;
1212 lsp[i] = (PI/order)*(float)i; 690 }
1213 }
1214 691
1215 return E; 692 return E;
1216} 693}
1217 694
1218/*---------------------------------------------------------------------------*\ 695/*---------------------------------------------------------------------------*\
@@ -1221,317 +698,60 @@ float speech_to_uq_lsps(float lsp[],
1221 AUTHOR......: David Rowe 698 AUTHOR......: David Rowe
1222 DATE CREATED: 22/8/2010 699 DATE CREATED: 22/8/2010
1223 700
1224 Thirty-six bit sclar LSP quantiser. From a vector of unquantised 701 Scalar LSP quantiser. From a vector of unquantised (floating point)
1225 (floating point) LSPs finds the quantised LSP indexes. 702 LSPs finds the quantised LSP indexes.
1226
1227\*---------------------------------------------------------------------------*/
1228
1229void encode_lsps_scalar(int indexes[], float lsp[], int order)
1230{
1231 int i,k,m;
1232 float wt[1];
1233 float lsp_hz[order];
1234 const float * cb;
1235 float se;
1236
1237 /* convert from radians to Hz so we can use human readable
1238 frequencies */
1239
1240 for(i=0; i<order; i++)
1241 lsp_hz[i] = (4000.0/PI)*lsp[i];
1242
1243 /* scalar quantisers */
1244
1245 wt[0] = 1.0;
1246 for(i=0; i<order; i++) {
1247 k = lsp_cb[i].k;
1248 m = lsp_cb[i].m;
1249 cb = lsp_cb[i].cb;
1250 indexes[i] = quantise(cb, &lsp_hz[i], wt, k, m, &se);
1251 }
1252}
1253
1254/*---------------------------------------------------------------------------*\
1255
1256 FUNCTION....: decode_lsps_scalar()
1257 AUTHOR......: David Rowe
1258 DATE CREATED: 22/8/2010
1259
1260 From a vector of quantised LSP indexes, returns the quantised
1261 (floating point) LSPs.
1262
1263\*---------------------------------------------------------------------------*/
1264
1265void decode_lsps_scalar(float lsp[], int indexes[], int order)
1266{
1267 int i,k;
1268 float lsp_hz[order];
1269 const float * cb;
1270
1271 for(i=0; i<order; i++) {
1272 k = lsp_cb[i].k;
1273 cb = lsp_cb[i].cb;
1274 lsp_hz[i] = cb[indexes[i]*k];
1275 }
1276
1277 /* convert back to radians */
1278
1279 for(i=0; i<order; i++)
1280 lsp[i] = (PI/4000.0)*lsp_hz[i];
1281}
1282
1283
1284/*---------------------------------------------------------------------------*\
1285
1286 FUNCTION....: encode_mels_scalar()
1287 AUTHOR......: David Rowe
1288 DATE CREATED: April 2015
1289
1290 Low bit rate mel coeff encoder.
1291 703
1292\*---------------------------------------------------------------------------*/ 704\*---------------------------------------------------------------------------*/
1293 705
1294void encode_mels_scalar(int indexes[], float mels[], int order) 706void encode_lsps_scalar(int indexes[], float lsp[], int order) {
1295{ 707 int i, k, m;
1296 int i,m; 708 float wt[1];
1297 float wt[1]; 709 float lsp_hz[order];
1298 const float * cb; 710 const float *cb;
1299 float se, mel_, dmel; 711 float se;
1300
1301 /* scalar quantisers */
1302
1303 wt[0] = 1.0;
1304 for(i=0; i<order; i++) {
1305 m = mel_cb[i].m;
1306 cb = mel_cb[i].cb;
1307 if (i%2) {
1308 /* on odd mels quantise difference */
1309 mel_ = mel_cb[i-1].cb[indexes[i-1]];
1310 dmel = mels[i] - mel_;
1311 indexes[i] = quantise(cb, &dmel, wt, 1, m, &se);
1312 //printf("%d mel: %f mel_: %f dmel: %f index: %d\n", i, mels[i], mel_, dmel, indexes[i]);
1313 }
1314 else {
1315 indexes[i] = quantise(cb, &mels[i], wt, 1, m, &se);
1316 //printf("%d mel: %f dmel: %f index: %d\n", i, mels[i], 0.0, indexes[i]);
1317 }
1318 712
1319 } 713 /* convert from radians to Hz so we can use human readable
1320} 714 frequencies */
1321 715
716 for (i = 0; i < order; i++) lsp_hz[i] = (4000.0 / PI) * lsp[i];
1322 717
1323/*---------------------------------------------------------------------------*\ 718 /* scalar quantisers */
1324 719
1325 FUNCTION....: decode_mels_scalar() 720 wt[0] = 1.0;
1326 AUTHOR......: David Rowe 721 for (i = 0; i < order; i++) {
1327 DATE CREATED: April 2015 722 k = lsp_cb[i].k;
1328 723 m = lsp_cb[i].m;
1329 From a vector of quantised mel indexes, returns the quantised 724 cb = lsp_cb[i].cb;
1330 (floating point) mels. 725 indexes[i] = quantise(cb, &lsp_hz[i], wt, k, m, &se);
1331 726 }
1332\*---------------------------------------------------------------------------*/
1333
1334void decode_mels_scalar(float mels[], int indexes[], int order)
1335{
1336 int i;
1337 const float * cb;
1338
1339 for(i=0; i<order; i++) {
1340 cb = mel_cb[i].cb;
1341 if (i%2) {
1342 /* on odd mels quantise difference */
1343 mels[i] = mels[i-1] + cb[indexes[i]];
1344 }
1345 else
1346 mels[i] = cb[indexes[i]];
1347 }
1348
1349}
1350
1351
1352#ifdef __EXPERIMENTAL__
1353
1354/*---------------------------------------------------------------------------*\
1355
1356 FUNCTION....: encode_lsps_diff_freq_vq()
1357 AUTHOR......: David Rowe
1358 DATE CREATED: 15 November 2011
1359
1360 Twenty-five bit LSP quantiser. LSPs 1-4 are quantised with scalar
1361 LSP differences (in frequency, i.e difference from the previous
1362 LSP). LSPs 5-10 are quantised with a VQ trained generated using
1363 vqtrainjnd.c
1364
1365\*---------------------------------------------------------------------------*/
1366
1367void encode_lsps_diff_freq_vq(int indexes[], float lsp[], int order)
1368{
1369 int i,k,m;
1370 float lsp_hz[order];
1371 float lsp__hz[order];
1372 float dlsp[order];
1373 float dlsp_[order];
1374 float wt[order];
1375 const float * cb;
1376 float se;
1377
1378 for(i=0; i<order; i++) {
1379 wt[i] = 1.0;
1380 }
1381
1382 /* convert from radians to Hz so we can use human readable
1383 frequencies */
1384
1385 for(i=0; i<order; i++)
1386 lsp_hz[i] = (4000.0/PI)*lsp[i];
1387
1388 /* scalar quantisers for LSP differences 1..4 */
1389
1390 wt[0] = 1.0;
1391 for(i=0; i<4; i++) {
1392 if (i)
1393 dlsp[i] = lsp_hz[i] - lsp__hz[i-1];
1394 else
1395 dlsp[0] = lsp_hz[0];
1396
1397 k = lsp_cbd[i].k;
1398 m = lsp_cbd[i].m;
1399 cb = lsp_cbd[i].cb;
1400 indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se);
1401 dlsp_[i] = cb[indexes[i]*k];
1402
1403 if (i)
1404 lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
1405 else
1406 lsp__hz[0] = dlsp_[0];
1407 }
1408
1409 /* VQ LSPs 5,6,7,8,9,10 */
1410
1411 k = lsp_cbjnd[4].k;
1412 m = lsp_cbjnd[4].m;
1413 cb = lsp_cbjnd[4].cb;
1414 indexes[4] = quantise(cb, &lsp_hz[4], &wt[4], k, m, &se);
1415} 727}
1416 728
1417
1418/*---------------------------------------------------------------------------*\ 729/*---------------------------------------------------------------------------*\
1419 730
1420 FUNCTION....: decode_lsps_diff_freq_vq() 731 FUNCTION....: decode_lsps_scalar()
1421 AUTHOR......: David Rowe 732 AUTHOR......: David Rowe
1422 DATE CREATED: 15 Nov 2011 733 DATE CREATED: 22/8/2010
1423 734
1424 From a vector of quantised LSP indexes, returns the quantised 735 From a vector of quantised LSP indexes, returns the quantised
1425 (floating point) LSPs. 736 (floating point) LSPs.
1426 737
1427\*---------------------------------------------------------------------------*/ 738\*---------------------------------------------------------------------------*/
1428 739
1429void decode_lsps_diff_freq_vq(float lsp_[], int indexes[], int order) 740void decode_lsps_scalar(float lsp[], int indexes[], int order) {
1430{ 741 int i, k;
1431 int i,k,m; 742 float lsp_hz[order];
1432 float dlsp_[order]; 743 const float *cb;
1433 float lsp__hz[order];
1434 const float * cb;
1435
1436 /* scalar LSP differences */
1437
1438 for(i=0; i<4; i++) {
1439 cb = lsp_cbd[i].cb;
1440 dlsp_[i] = cb[indexes[i]];
1441 if (i)
1442 lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
1443 else
1444 lsp__hz[0] = dlsp_[0];
1445 }
1446
1447 /* VQ */
1448
1449 k = lsp_cbjnd[4].k;
1450 m = lsp_cbjnd[4].m;
1451 cb = lsp_cbjnd[4].cb;
1452 for(i=4; i<order; i++)
1453 lsp__hz[i] = cb[indexes[4]*k+i-4];
1454
1455 /* convert back to radians */
1456
1457 for(i=0; i<order; i++)
1458 lsp_[i] = (PI/4000.0)*lsp__hz[i];
1459}
1460
1461
1462/*---------------------------------------------------------------------------*\
1463
1464 FUNCTION....: encode_lsps_diff_time()
1465 AUTHOR......: David Rowe
1466 DATE CREATED: 12 Sep 2012
1467
1468 Encode difference from preious frames's LSPs using
1469 3,3,2,2,2,2,1,1,1,1 scalar quantisers (18 bits total).
1470
1471\*---------------------------------------------------------------------------*/
1472 744
1473void encode_lsps_diff_time(int indexes[], 745 for (i = 0; i < order; i++) {
1474 float lsps[], 746 k = lsp_cb[i].k;
1475 float lsps__prev[], 747 cb = lsp_cb[i].cb;
1476 int order) 748 lsp_hz[i] = cb[indexes[i] * k];
1477{ 749 }
1478 int i,k,m;
1479 float lsps_dt[order];
1480 float wt[LPC_MAX];
1481 const float * cb;
1482 float se;
1483
1484 /* Determine difference in time and convert from radians to Hz so
1485 we can use human readable frequencies */
1486
1487 for(i=0; i<order; i++) {
1488 lsps_dt[i] = (4000/PI)*(lsps[i] - lsps__prev[i]);
1489 }
1490
1491 /* scalar quantisers */
1492
1493 wt[0] = 1.0;
1494 for(i=0; i<order; i++) {
1495 k = lsp_cbdt[i].k;
1496 m = lsp_cbdt[i].m;
1497 cb = lsp_cbdt[i].cb;
1498 indexes[i] = quantise(cb, &lsps_dt[i], wt, k, m, &se);
1499 }
1500
1501}
1502
1503
1504/*---------------------------------------------------------------------------*\
1505
1506 FUNCTION....: decode_lsps_diff_time()
1507 AUTHOR......: David Rowe
1508 DATE CREATED: 15 Nov 2011
1509
1510 From a quantised LSP indexes, returns the quantised
1511 (floating point) LSPs.
1512
1513\*---------------------------------------------------------------------------*/
1514
1515void decode_lsps_diff_time(
1516 float lsps_[],
1517 int indexes[],
1518 float lsps__prev[],
1519 int order)
1520{
1521 int i,k,m;
1522 const float * cb;
1523
1524 for(i=0; i<order; i++)
1525 lsps_[i] = lsps__prev[i];
1526 750
1527 for(i=0; i<order; i++) { 751 /* convert back to radians */
1528 k = lsp_cbdt[i].k;
1529 cb = lsp_cbdt[i].cb;
1530 lsps_[i] += (PI/4000.0)*cb[indexes[i]*k];
1531 }
1532 752
753 for (i = 0; i < order; i++) lsp[i] = (PI / 4000.0) * lsp_hz[i];
1533} 754}
1534#endif
1535 755
1536/*---------------------------------------------------------------------------*\ 756/*---------------------------------------------------------------------------*\
1537 757
@@ -1543,45 +763,40 @@ void decode_lsps_diff_time(
1543 763
1544\*---------------------------------------------------------------------------*/ 764\*---------------------------------------------------------------------------*/
1545 765
1546void encode_lsps_vq(int *indexes, float *x, float *xq, int order) 766void encode_lsps_vq(int *indexes, float *x, float *xq, int order) {
1547{
1548 int i, n1, n2, n3; 767 int i, n1, n2, n3;
1549 float err[order], err2[order], err3[order]; 768 float err[order], err2[order], err3[order];
1550 float w[order], w2[order], w3[order]; 769 float w[order], w2[order], w3[order];
1551 const float *codebook1 = lsp_cbjvm[0].cb; 770 const float *codebook1 = lsp_cbjmv[0].cb;
1552 const float *codebook2 = lsp_cbjvm[1].cb; 771 const float *codebook2 = lsp_cbjmv[1].cb;
1553 const float *codebook3 = lsp_cbjvm[2].cb; 772 const float *codebook3 = lsp_cbjmv[2].cb;
1554 773
1555 w[0] = MIN(x[0], x[1]-x[0]); 774 w[0] = MIN(x[0], x[1] - x[0]);
1556 for (i=1;i<order-1;i++) 775 for (i = 1; i < order - 1; i++) w[i] = MIN(x[i] - x[i - 1], x[i + 1] - x[i]);
1557 w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]); 776 w[order - 1] = MIN(x[order - 1] - x[order - 2], PI - x[order - 1]);
1558 w[order-1] = MIN(x[order-1]-x[order-2], PI-x[order-1]);
1559 777
1560 compute_weights(x, w, order); 778 compute_weights(x, w, order);
1561 779
1562 n1 = find_nearest(codebook1, lsp_cbjvm[0].m, x, order); 780 n1 = find_nearest(codebook1, lsp_cbjmv[0].m, x, order);
1563 781
1564 for (i=0;i<order;i++) 782 for (i = 0; i < order; i++) {
1565 { 783 xq[i] = codebook1[order * n1 + i];
1566 xq[i] = codebook1[order*n1+i];
1567 err[i] = x[i] - xq[i]; 784 err[i] = x[i] - xq[i];
1568 } 785 }
1569 for (i=0;i<order/2;i++) 786 for (i = 0; i < order / 2; i++) {
1570 { 787 err2[i] = err[2 * i];
1571 err2[i] = err[2*i]; 788 err3[i] = err[2 * i + 1];
1572 err3[i] = err[2*i+1]; 789 w2[i] = w[2 * i];
1573 w2[i] = w[2*i]; 790 w3[i] = w[2 * i + 1];
1574 w3[i] = w[2*i+1];
1575 } 791 }
1576 n2 = find_nearest_weighted(codebook2, lsp_cbjvm[1].m, err2, w2, order/2); 792 n2 = find_nearest_weighted(codebook2, lsp_cbjmv[1].m, err2, w2, order / 2);
1577 n3 = find_nearest_weighted(codebook3, lsp_cbjvm[2].m, err3, w3, order/2); 793 n3 = find_nearest_weighted(codebook3, lsp_cbjmv[2].m, err3, w3, order / 2);
1578 794
1579 indexes[0] = n1; 795 indexes[0] = n1;
1580 indexes[1] = n2; 796 indexes[1] = n2;
1581 indexes[2] = n3; 797 indexes[2] = n3;
1582} 798}
1583 799
1584
1585/*---------------------------------------------------------------------------*\ 800/*---------------------------------------------------------------------------*\
1586 801
1587 FUNCTION....: decode_lsps_vq() 802 FUNCTION....: decode_lsps_vq()
@@ -1590,31 +805,28 @@ void encode_lsps_vq(int *indexes, float *x, float *xq, int order)
1590 805
1591\*---------------------------------------------------------------------------*/ 806\*---------------------------------------------------------------------------*/
1592 807
1593void decode_lsps_vq(int *indexes, float *xq, int order, int stages) 808void decode_lsps_vq(int *indexes, float *xq, int order, int stages) {
1594{
1595 int i, n1, n2, n3; 809 int i, n1, n2, n3;
1596 const float *codebook1 = lsp_cbjvm[0].cb; 810 const float *codebook1 = lsp_cbjmv[0].cb;
1597 const float *codebook2 = lsp_cbjvm[1].cb; 811 const float *codebook2 = lsp_cbjmv[1].cb;
1598 const float *codebook3 = lsp_cbjvm[2].cb; 812 const float *codebook3 = lsp_cbjmv[2].cb;
1599 813
1600 n1 = indexes[0]; 814 n1 = indexes[0];
1601 n2 = indexes[1]; 815 n2 = indexes[1];
1602 n3 = indexes[2]; 816 n3 = indexes[2];
1603 817
1604 for (i=0;i<order;i++) { 818 for (i = 0; i < order; i++) {
1605 xq[i] = codebook1[order*n1+i]; 819 xq[i] = codebook1[order * n1 + i];
1606 } 820 }
1607 821
1608 if (stages != 1) { 822 if (stages != 1) {
1609 for (i=0;i<order/2;i++) { 823 for (i = 0; i < order / 2; i++) {
1610 xq[2*i] += codebook2[order*n2/2+i]; 824 xq[2 * i] += codebook2[order * n2 / 2 + i];
1611 xq[2*i+1] += codebook3[order*n3/2+i]; 825 xq[2 * i + 1] += codebook3[order * n3 / 2 + i];
1612 } 826 }
1613 } 827 }
1614
1615} 828}
1616 829
1617
1618/*---------------------------------------------------------------------------*\ 830/*---------------------------------------------------------------------------*\
1619 831
1620 FUNCTION....: bw_expand_lsps() 832 FUNCTION....: bw_expand_lsps()
@@ -1628,124 +840,45 @@ void decode_lsps_vq(int *indexes, float *xq, int order, int stages)
1628 840
1629\*---------------------------------------------------------------------------*/ 841\*---------------------------------------------------------------------------*/
1630 842
1631void bw_expand_lsps(float lsp[], int order, float min_sep_low, float min_sep_high) 843void bw_expand_lsps(float lsp[], int order, float min_sep_low,
1632{ 844 float min_sep_high) {
1633 int i; 845 int i;
1634
1635 for(i=1; i<4; i++) {
1636
1637 if ((lsp[i] - lsp[i-1]) < min_sep_low*(PI/4000.0))
1638 lsp[i] = lsp[i-1] + min_sep_low*(PI/4000.0);
1639
1640 }
1641
1642 /* As quantiser gaps increased, larger BW expansion was required
1643 to prevent twinkly noises. This may need more experiment for
1644 different quanstisers.
1645 */
1646
1647 for(i=4; i<order; i++) {
1648 if (lsp[i] - lsp[i-1] < min_sep_high*(PI/4000.0))
1649 lsp[i] = lsp[i-1] + min_sep_high*(PI/4000.0);
1650 }
1651}
1652
1653void bw_expand_lsps2(float lsp[],
1654 int order
1655)
1656{
1657 int i;
1658
1659 for(i=1; i<4; i++) {
1660
1661 if ((lsp[i] - lsp[i-1]) < 100.0*(PI/4000.0))
1662 lsp[i] = lsp[i-1] + 100.0*(PI/4000.0);
1663 846
1664 } 847 for (i = 1; i < 4; i++) {
848 if ((lsp[i] - lsp[i - 1]) < min_sep_low * (PI / 4000.0))
849 lsp[i] = lsp[i - 1] + min_sep_low * (PI / 4000.0);
850 }
1665 851
1666 /* As quantiser gaps increased, larger BW expansion was required 852 /* As quantiser gaps increased, larger BW expansion was required
1667 to prevent twinkly noises. This may need more experiment for 853 to prevent twinkly noises. This may need more experiment for
1668 different quanstisers. 854 different quanstisers.
1669 */ 855 */
1670 856
1671 for(i=4; i<order; i++) { 857 for (i = 4; i < order; i++) {
1672 if (lsp[i] - lsp[i-1] < 200.0*(PI/4000.0)) 858 if (lsp[i] - lsp[i - 1] < min_sep_high * (PI / 4000.0))
1673 lsp[i] = lsp[i-1] + 200.0*(PI/4000.0); 859 lsp[i] = lsp[i - 1] + min_sep_high * (PI / 4000.0);
1674 } 860 }
1675} 861}
1676 862
1677/*---------------------------------------------------------------------------*\ 863void bw_expand_lsps2(float lsp[], int order) {
1678 864 int i;
1679 FUNCTION....: locate_lsps_jnd_steps()
1680 AUTHOR......: David Rowe
1681 DATE CREATED: 27/10/2011
1682
1683 Applies a form of Bandwidth Expansion (BW) to a vector of LSPs.
1684 Listening tests have determined that "quantising" the position of
1685 each LSP to the non-linear steps below introduces a "just noticable
1686 difference" in the synthesised speech.
1687
1688 This operation can be used before quantisation to limit the input
1689 data to the quantiser to a number of discrete steps.
1690
1691 This operation can also be used during quantisation as a form of
1692 hysteresis in the calculation of quantiser error. For example if
1693 the quantiser target of lsp1 is 500 Hz, candidate vectors with lsp1
1694 of 515 and 495 Hz sound effectively the same.
1695
1696\*---------------------------------------------------------------------------*/
1697
1698void locate_lsps_jnd_steps(float lsps[], int order)
1699{
1700 int i;
1701 float lsp_hz, step;
1702
1703 assert(order == 10);
1704
1705 /* quantise to 25Hz steps */
1706
1707 step = 25;
1708 for(i=0; i<2; i++) {
1709 lsp_hz = lsps[i]*4000.0/PI;
1710 lsp_hz = floorf(lsp_hz/step + 0.5)*step;
1711 lsps[i] = lsp_hz*PI/4000.0;
1712 if (i) {
1713 if (lsps[i] == lsps[i-1])
1714 lsps[i] += step*PI/4000.0;
1715
1716 }
1717 }
1718
1719 /* quantise to 50Hz steps */
1720
1721 step = 50;
1722 for(i=2; i<4; i++) {
1723 lsp_hz = lsps[i]*4000.0/PI;
1724 lsp_hz = floorf(lsp_hz/step + 0.5)*step;
1725 lsps[i] = lsp_hz*PI/4000.0;
1726 if (i) {
1727 if (lsps[i] == lsps[i-1])
1728 lsps[i] += step*PI/4000.0;
1729
1730 }
1731 }
1732 865
1733 /* quantise to 100Hz steps */ 866 for (i = 1; i < 4; i++) {
867 if ((lsp[i] - lsp[i - 1]) < 100.0 * (PI / 4000.0))
868 lsp[i] = lsp[i - 1] + 100.0 * (PI / 4000.0);
869 }
1734 870
1735 step = 100; 871 /* As quantiser gaps increased, larger BW expansion was required
1736 for(i=4; i<10; i++) { 872 to prevent twinkly noises. This may need more experiment for
1737 lsp_hz = lsps[i]*4000.0/PI; 873 different quanstisers.
1738 lsp_hz = floorf(lsp_hz/step + 0.5)*step; 874 */
1739 lsps[i] = lsp_hz*PI/4000.0;
1740 if (i) {
1741 if (lsps[i] == lsps[i-1])
1742 lsps[i] += step*PI/4000.0;
1743 875
1744 } 876 for (i = 4; i < order; i++) {
1745 } 877 if (lsp[i] - lsp[i - 1] < 200.0 * (PI / 4000.0))
878 lsp[i] = lsp[i - 1] + 200.0 * (PI / 4000.0);
879 }
1746} 880}
1747 881
1748
1749/*---------------------------------------------------------------------------*\ 882/*---------------------------------------------------------------------------*\
1750 883
1751 FUNCTION....: apply_lpc_correction() 884 FUNCTION....: apply_lpc_correction()
@@ -1757,11 +890,10 @@ void locate_lsps_jnd_steps(float lsps[], int order)
1757 890
1758\*---------------------------------------------------------------------------*/ 891\*---------------------------------------------------------------------------*/
1759 892
1760void apply_lpc_correction(MODEL *model) 893void apply_lpc_correction(MODEL *model) {
1761{ 894 if (model->Wo < (PI * 150.0 / 4000)) {
1762 if (model->Wo < (PI*150.0/4000)) { 895 model->A[1] *= 0.032;
1763 model->A[1] *= 0.032; 896 }
1764 }
1765} 897}
1766 898
1767/*---------------------------------------------------------------------------*\ 899/*---------------------------------------------------------------------------*\
@@ -1774,124 +906,81 @@ void apply_lpc_correction(MODEL *model)
1774 906
1775\*---------------------------------------------------------------------------*/ 907\*---------------------------------------------------------------------------*/
1776 908
1777int encode_energy(float e, int bits) 909int encode_energy(float e, int bits) {
1778{ 910 int index, e_levels = 1 << bits;
1779 int index, e_levels = 1<<bits; 911 float e_min = E_MIN_DB;
1780 float e_min = E_MIN_DB; 912 float e_max = E_MAX_DB;
1781 float e_max = E_MAX_DB; 913 float norm;
1782 float norm;
1783
1784 e = 10.0*log10f(e);
1785 norm = (e - e_min)/(e_max - e_min);
1786 index = floorf(e_levels * norm + 0.5);
1787 if (index < 0 ) index = 0;
1788 if (index > (e_levels-1)) index = e_levels-1;
1789
1790 return index;
1791}
1792
1793/*---------------------------------------------------------------------------*\
1794
1795 FUNCTION....: decode_energy()
1796 AUTHOR......: David Rowe
1797 DATE CREATED: 22/8/2010
1798
1799 Decodes energy using a E_LEVELS quantiser.
1800
1801\*---------------------------------------------------------------------------*/
1802
1803float decode_energy(int index, int bits)
1804{
1805 float e_min = E_MIN_DB;
1806 float e_max = E_MAX_DB;
1807 float step;
1808 float e;
1809 int e_levels = 1<<bits;
1810 914
1811 step = (e_max - e_min)/e_levels; 915 e = 10.0 * log10f(e);
1812 e = e_min + step*(index); 916 norm = (e - e_min) / (e_max - e_min);
1813 e = POW10F(e/10.0); 917 index = floorf(e_levels * norm + 0.5);
918 if (index < 0) index = 0;
919 if (index > (e_levels - 1)) index = e_levels - 1;
1814 920
1815 return e; 921 return index;
1816} 922}
1817 923
1818#ifdef NOT_USED
1819/*---------------------------------------------------------------------------*\ 924/*---------------------------------------------------------------------------*\
1820 925
1821 FUNCTION....: decode_amplitudes() 926 FUNCTION....: decode_energy()
1822 AUTHOR......: David Rowe 927 AUTHOR......: David Rowe
1823 DATE CREATED: 22/8/2010 928 DATE CREATED: 22/8/2010
1824 929
1825 Given the amplitude quantiser indexes recovers the harmonic 930 Decodes energy using a E_LEVELS quantiser.
1826 amplitudes.
1827 931
1828\*---------------------------------------------------------------------------*/ 932\*---------------------------------------------------------------------------*/
1829 933
1830float decode_amplitudes(codec2_fft_cfg fft_fwd_cfg, 934float decode_energy(int index, int bits) {
1831 MODEL *model, 935 float e_min = E_MIN_DB;
1832 float ak[], 936 float e_max = E_MAX_DB;
1833 int lsp_indexes[], 937 float step;
1834 int energy_index, 938 float e;
1835 float lsps[], 939 int e_levels = 1 << bits;
1836 float *e
1837)
1838{
1839 float snr;
1840 940
1841 decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD); 941 step = (e_max - e_min) / e_levels;
1842 bw_expand_lsps(lsps, LPC_ORD); 942 e = e_min + step * (index);
1843 lsp_to_lpc(lsps, ak, LPC_ORD); 943 e = POW10F(e / 10.0);
1844 *e = decode_energy(energy_index);
1845 aks_to_M2(ak, LPC_ORD, model, *e, &snr, 1, 0, 0, 1);
1846 apply_lpc_correction(model);
1847 944
1848 return snr; 945 return e;
1849} 946}
1850#endif
1851 947
1852static float ge_coeff[2] = {0.8, 0.9}; 948static float ge_coeff[2] = {0.8, 0.9};
1853 949
1854void compute_weights2(const float *x, const float *xp, float *w) 950void compute_weights2(const float *x, const float *xp, float *w) {
1855{
1856 w[0] = 30; 951 w[0] = 30;
1857 w[1] = 1; 952 w[1] = 1;
1858 if (x[1]<0) 953 if (x[1] < 0) {
1859 { 954 w[0] *= .6;
1860 w[0] *= .6; 955 w[1] *= .3;
1861 w[1] *= .3;
1862 } 956 }
1863 if (x[1]<-10) 957 if (x[1] < -10) {
1864 { 958 w[0] *= .3;
1865 w[0] *= .3; 959 w[1] *= .3;
1866 w[1] *= .3;
1867 } 960 }
1868 /* Higher weight if pitch is stable */ 961 /* Higher weight if pitch is stable */
1869 if (fabsf(x[0]-xp[0])<.2) 962 if (fabsf(x[0] - xp[0]) < .2) {
1870 { 963 w[0] *= 2;
1871 w[0] *= 2; 964 w[1] *= 1.5;
1872 w[1] *= 1.5; 965 } else if (fabsf(x[0] - xp[0]) > .5) /* Lower if not stable */
1873 } else if (fabsf(x[0]-xp[0])>.5) /* Lower if not stable */
1874 { 966 {
1875 w[0] *= .5; 967 w[0] *= .5;
1876 } 968 }
1877 969
1878 /* Lower weight for low energy */ 970 /* Lower weight for low energy */
1879 if (x[1] < xp[1]-10) 971 if (x[1] < xp[1] - 10) {
1880 { 972 w[1] *= .5;
1881 w[1] *= .5;
1882 } 973 }
1883 if (x[1] < xp[1]-20) 974 if (x[1] < xp[1] - 20) {
1884 { 975 w[1] *= .5;
1885 w[1] *= .5;
1886 } 976 }
1887 977
1888 //w[0] = 30; 978 // w[0] = 30;
1889 //w[1] = 1; 979 // w[1] = 1;
1890 980
1891 /* Square the weights because it's applied on the squared error */ 981 /* Square the weights because it's applied on the squared error */
1892 w[0] *= w[0]; 982 w[0] *= w[0];
1893 w[1] *= w[1]; 983 w[1] *= w[1];
1894
1895} 984}
1896 985
1897/*---------------------------------------------------------------------------*\ 986/*---------------------------------------------------------------------------*\
@@ -1906,7 +995,7 @@ void compute_weights2(const float *x, const float *xp, float *w)
1906 both the pitch and energy tend to only change by small amounts 995 both the pitch and energy tend to only change by small amounts
1907 during voiced speech, however it is important that these changes be 996 during voiced speech, however it is important that these changes be
1908 coded carefully. During unvoiced speech they both change a lot but 997 coded carefully. During unvoiced speech they both change a lot but
1909 the ear is less sensitve to errors so coarser quantisation is OK. 998 the ear is less sensitive to errors so coarser quantisation is OK.
1910 999
1911 The ear is sensitive to log energy and loq pitch so we quantise in 1000 The ear is sensitive to log energy and loq pitch so we quantise in
1912 these domains. That way the error measure used to quantise the 1001 these domains. That way the error measure used to quantise the
@@ -1916,15 +1005,14 @@ void compute_weights2(const float *x, const float *xp, float *w)
1916 1005
1917\*---------------------------------------------------------------------------*/ 1006\*---------------------------------------------------------------------------*/
1918 1007
1919void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) 1008void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) {
1920{ 1009 int i, n1;
1921 int i, n1; 1010 float x[2];
1922 float x[2]; 1011 float err[2];
1923 float err[2]; 1012 float w[2];
1924 float w[2];
1925 const float *codebook1 = ge_cb[0].cb; 1013 const float *codebook1 = ge_cb[0].cb;
1926 int nb_entries = ge_cb[0].m; 1014 int nb_entries = ge_cb[0].m;
1927 int ndim = ge_cb[0].k; 1015 int ndim = ge_cb[0].k;
1928 float Wo_min = c2const->Wo_min; 1016 float Wo_min = c2const->Wo_min;
1929 float Wo_max = c2const->Wo_max; 1017 float Wo_max = c2const->Wo_max;
1930 float Fs = c2const->Fs; 1018 float Fs = c2const->Fs;
@@ -1933,18 +1021,16 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[])
1933 1021
1934 assert(Fs == 8000); 1022 assert(Fs == 8000);
1935 1023
1936 x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2); 1024 x[0] = log10f((model->Wo / PI) * 4000.0 / 50.0) / log10f(2);
1937 x[1] = 10.0*log10f(1e-4 + *e); 1025 x[1] = 10.0 * log10f(1e-4 + *e);
1938 1026
1939 compute_weights2(x, xq, w); 1027 compute_weights2(x, xq, w);
1940 for (i=0;i<ndim;i++) 1028 for (i = 0; i < ndim; i++) err[i] = x[i] - ge_coeff[i] * xq[i];
1941 err[i] = x[i]-ge_coeff[i]*xq[i];
1942 n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim); 1029 n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim);
1943 1030
1944 for (i=0;i<ndim;i++) 1031 for (i = 0; i < ndim; i++) {
1945 { 1032 xq[i] = ge_coeff[i] * xq[i] + codebook1[ndim * n1 + i];
1946 xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i]; 1033 err[i] -= codebook1[ndim * n1 + i];
1947 err[i] -= codebook1[ndim*n1+i];
1948 } 1034 }
1949 1035
1950 /* 1036 /*
@@ -1953,7 +1039,7 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[])
1953 Wo = (2^x)*(PI*50)/4000; 1039 Wo = (2^x)*(PI*50)/4000;
1954 */ 1040 */
1955 1041
1956 model->Wo = powf(2.0, xq[0])*(PI*50.0)/4000.0; 1042 model->Wo = powf(2.0, xq[0]) * (PI * 50.0) / 4000.0;
1957 1043
1958 /* bit errors can make us go out of range leading to all sorts of 1044 /* bit errors can make us go out of range leading to all sorts of
1959 probs like seg faults */ 1045 probs like seg faults */
@@ -1961,9 +1047,9 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[])
1961 if (model->Wo > Wo_max) model->Wo = Wo_max; 1047 if (model->Wo > Wo_max) model->Wo = Wo_max;
1962 if (model->Wo < Wo_min) model->Wo = Wo_min; 1048 if (model->Wo < Wo_min) model->Wo = Wo_min;
1963 1049
1964 model->L = PI/model->Wo; /* if we quantise Wo re-compute L */ 1050 model->L = PI / model->Wo; /* if we quantise Wo re-compute L */
1965 1051
1966 *e = POW10F(xq[1]/10.0); 1052 *e = POW10F(xq[1] / 10.0);
1967} 1053}
1968 1054
1969/*---------------------------------------------------------------------------*\ 1055/*---------------------------------------------------------------------------*\
@@ -1977,39 +1063,36 @@ void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[])
1977 1063
1978\*---------------------------------------------------------------------------*/ 1064\*---------------------------------------------------------------------------*/
1979 1065
1980int encode_WoE(MODEL *model, float e, float xq[]) 1066int encode_WoE(MODEL *model, float e, float xq[]) {
1981{ 1067 int i, n1;
1982 int i, n1; 1068 float x[2];
1983 float x[2]; 1069 float err[2];
1984 float err[2]; 1070 float w[2];
1985 float w[2];
1986 const float *codebook1 = ge_cb[0].cb; 1071 const float *codebook1 = ge_cb[0].cb;
1987 int nb_entries = ge_cb[0].m; 1072 int nb_entries = ge_cb[0].m;
1988 int ndim = ge_cb[0].k; 1073 int ndim = ge_cb[0].k;
1989 1074
1990 assert((1<<WO_E_BITS) == nb_entries); 1075 assert((1 << WO_E_BITS) == nb_entries);
1991 1076
1992 if (e < 0.0) e = 0; /* occasional small negative energies due LPC round off I guess */ 1077 if (e < 0.0)
1078 e = 0; /* occasional small negative energies due LPC round off I guess */
1993 1079
1994 x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2); 1080 x[0] = log10f((model->Wo / PI) * 4000.0 / 50.0) / log10f(2);
1995 x[1] = 10.0*log10f(1e-4 + e); 1081 x[1] = 10.0 * log10f(1e-4 + e);
1996 1082
1997 compute_weights2(x, xq, w); 1083 compute_weights2(x, xq, w);
1998 for (i=0;i<ndim;i++) 1084 for (i = 0; i < ndim; i++) err[i] = x[i] - ge_coeff[i] * xq[i];
1999 err[i] = x[i]-ge_coeff[i]*xq[i];
2000 n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim); 1085 n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim);
2001 1086
2002 for (i=0;i<ndim;i++) 1087 for (i = 0; i < ndim; i++) {
2003 { 1088 xq[i] = ge_coeff[i] * xq[i] + codebook1[ndim * n1 + i];
2004 xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i]; 1089 err[i] -= codebook1[ndim * n1 + i];
2005 err[i] -= codebook1[ndim*n1+i];
2006 } 1090 }
2007 1091
2008 //printf("enc: %f %f (%f)(%f) \n", xq[0], xq[1], e, 10.0*log10(1e-4 + e)); 1092 // printf("enc: %f %f (%f)(%f) \n", xq[0], xq[1], e, 10.0*log10(1e-4 + e));
2009 return n1; 1093 return n1;
2010} 1094}
2011 1095
2012
2013/*---------------------------------------------------------------------------*\ 1096/*---------------------------------------------------------------------------*\
2014 1097
2015 FUNCTION....: decode_WoE() 1098 FUNCTION....: decode_WoE()
@@ -2022,21 +1105,19 @@ int encode_WoE(MODEL *model, float e, float xq[])
2022 1105
2023\*---------------------------------------------------------------------------*/ 1106\*---------------------------------------------------------------------------*/
2024 1107
2025void decode_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[], int n1) 1108void decode_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[], int n1) {
2026{ 1109 int i;
2027 int i;
2028 const float *codebook1 = ge_cb[0].cb; 1110 const float *codebook1 = ge_cb[0].cb;
2029 int ndim = ge_cb[0].k; 1111 int ndim = ge_cb[0].k;
2030 float Wo_min = c2const->Wo_min; 1112 float Wo_min = c2const->Wo_min;
2031 float Wo_max = c2const->Wo_max; 1113 float Wo_max = c2const->Wo_max;
2032 1114
2033 for (i=0;i<ndim;i++) 1115 for (i = 0; i < ndim; i++) {
2034 { 1116 xq[i] = ge_coeff[i] * xq[i] + codebook1[ndim * n1 + i];
2035 xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i];
2036 } 1117 }
2037 1118
2038 //printf("dec: %f %f\n", xq[0], xq[1]); 1119 // printf("dec: %f %f\n", xq[0], xq[1]);
2039 model->Wo = powf(2.0, xq[0])*(PI*50.0)/4000.0; 1120 model->Wo = powf(2.0, xq[0]) * (PI * 50.0) / 4000.0;
2040 1121
2041 /* bit errors can make us go out of range leading to all sorts of 1122 /* bit errors can make us go out of range leading to all sorts of
2042 probs like seg faults */ 1123 probs like seg faults */
@@ -2044,8 +1125,7 @@ void decode_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[], int n1)
2044 if (model->Wo > Wo_max) model->Wo = Wo_max; 1125 if (model->Wo > Wo_max) model->Wo = Wo_max;
2045 if (model->Wo < Wo_min) model->Wo = Wo_min; 1126 if (model->Wo < Wo_min) model->Wo = Wo_min;
2046 1127
2047 model->L = PI/model->Wo; /* if we quantise Wo re-compute L */ 1128 model->L = PI / model->Wo; /* if we quantise Wo re-compute L */
2048 1129
2049 *e = POW10F(xq[1]/10.0); 1130 *e = POW10F(xq[1] / 10.0);
2050} 1131}
2051
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-12 02:17:13 +0000