« back to all changes in this revision
Viewing changes to gr-vocoder/lib/codec2/quantise.c
- browse files at revision 20
- compare with another revision
- download diff
- download tarball
- view history from revision 20
- Committer: Package Import Robot
- Author(s): A. Maitland Bottoms
- Date: 2014-01-05 12:14:43 UTC
- mfrom: (2.1.23 sid)
- Revision ID: package-import@ubuntu.com-20140105121443-je18dkof6uhox808
Tags: 3.7.2.1-5
use correct compiler for unstable build
- files added:
- gr-vocoder/lib/codec2/codebook/lsp8910.txt
- files removed:
- .pc/Feature606
- .pc/Feature606/grc
- .pc/update-codec2
- .pc/update-codec2/gr-vocoder
- .pc/update-codec2/gr-vocoder/examples
- .pc/update-codec2/gr-vocoder/grc
- .pc/update-codec2/gr-vocoder/include
- .pc/update-codec2/gr-vocoder/include/gnuradio
- .pc/update-codec2/gr-vocoder/include/gnuradio/vocoder
- .pc/update-codec2/gr-vocoder/lib
- .pc/update-codec2/gr-vocoder/lib/codec2
- .pc/update-codec2/gr-vocoder/lib/codec2/codebook
- .pc/update-codec2/gr-vocoder/python
- .pc/update-codec2/gr-vocoder/python/vocoder
- files modified:
- .pc/applied-patches
added
removed
gr-vocoder/lib/codec2/quantise.c
36
36
#include "quantise.h"
37
37
#include "lpc.h"
38
38
#include "lsp.h"
39
#include "kiss_fft.h"
40
#undef TIMER
41
#include "machdep.h"
39
#include "fft.h"
42
40
43
41
#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */
44
42
61
59
return lsp_cb[i].log2m;
62
60
}
63
61
64
int lspd_bits(int i) {
65
return lsp_cbd[i].log2m;
66
}
67
68
#ifdef __EXPERIMENTAL__
69
int lspdt_bits(int i) {
70
return lsp_cbdt[i].log2m;
71
}
62
#if VECTOR_QUANTISATION
63
/*---------------------------------------------------------------------------*\
64
65
quantise_uniform
66
67
Simulates uniform quantising of a float.
68
69
\*---------------------------------------------------------------------------*/
70
71
void quantise_uniform(float *val, float min, float max, int bits)
72
{
73
int levels = 1 << (bits-1);
74
float norm;
75
int index;
76
77
/* hard limit to quantiser range */
78
79
printf("min: %f max: %f val: %f ", min, max, val[0]);
80
if (val[0] < min) val[0] = min;
81
if (val[0] > max) val[0] = max;
82
83
norm = (*val - min)/(max-min);
84
printf("%f norm: %f ", val[0], norm);
85
index = fabs(levels*norm + 0.5);
86
87
*val = min + index*(max-min)/levels;
88
89
printf("index %d val_: %f\n", index, val[0]);
90
}
91
72
92
#endif
73
93
74
int lsp_pred_vq_bits(int i) {
75
return lsp_cbjvm[i].log2m;
76
}
77
78
94
/*---------------------------------------------------------------------------*\
79
95
80
96
quantise_init
111
127
float beste; /* best error so far */
112
128
long j;
113
129
int i;
114
float diff;
115
130
116
131
besti = 0;
117
132
beste = 1E32;
118
133
for(j=0; j<m; j++) {
119
134
e = 0.0;
120
for(i=0; i<k; i++) {
121
diff = cb[j*k+i]-vec[i];
122
e += powf(diff*w[i],2.0);
123
}
135
for(i=0; i<k; i++)
136
e += pow((cb[j*k+i]-vec[i])*w[i],2.0);
124
137
if (e < beste) {
125
138
beste = e;
126
139
besti = j;
134
147
135
148
/*---------------------------------------------------------------------------*\
136
149
137
encode_lspds_scalar()
150
lspd_quantise
138
151
139
Scalar/VQ LSP difference quantiser.
152
Scalar lsp difference quantiser.
140
153
141
154
\*---------------------------------------------------------------------------*/
142
155
143
void encode_lspds_scalar(
144
int indexes[],
145
float lsp[],
146
int order
156
void lspd_quantise(
157
float lsp[],
158
float lsp_[],
159
int order
147
160
)
148
161
{
149
162
int i,k,m;
151
164
float lsp__hz[LPC_MAX];
152
165
float dlsp[LPC_MAX];
153
166
float dlsp_[LPC_MAX];
154
float wt[LPC_MAX];
167
float wt[1];
155
168
const float *cb;
156
float se;
157
158
assert(order == LPC_ORD);
159
160
for(i=0; i<order; i++) {
161
wt[i] = 1.0;
162
}
169
float se = 0.0;
170
int indexes[LPC_MAX];
163
171
164
172
/* convert from radians to Hz so we can use human readable
165
173
frequencies */
167
175
for(i=0; i<order; i++)
168
176
lsp_hz[i] = (4000.0/PI)*lsp[i];
169
177
170
//printf("\n");
178
dlsp[0] = lsp_hz[0];
179
for(i=1; i<order; i++)
180
dlsp[i] = lsp_hz[i] - lsp_hz[i-1];
181
182
/* simple uniform scalar quantisers */
171
183
172
184
wt[0] = 1.0;
173
185
for(i=0; i<order; i++) {
174
175
/* find difference from previous qunatised lsp */
176
177
186
if (i)
178
187
dlsp[i] = lsp_hz[i] - lsp__hz[i-1];
179
188
else
185
194
indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se);
186
195
dlsp_[i] = cb[indexes[i]*k];
187
196
188
189
197
if (i)
190
198
lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
191
199
else
192
200
lsp__hz[0] = dlsp_[0];
193
194
//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]);
195
201
}
196
197
}
198
199
void decode_lspds_scalar(
200
float lsp_[],
201
int indexes[],
202
int order
203
)
204
{
205
int i,k;
206
float lsp__hz[LPC_MAX];
207
float dlsp_[LPC_MAX];
208
const float *cb;
209
210
assert(order == LPC_ORD);
211
212
for(i=0; i<order; i++) {
213
214
k = lsp_cbd[i].k;
215
cb = lsp_cbd[i].cb;
216
dlsp_[i] = cb[indexes[i]*k];
217
218
if (i)
219
lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
220
else
221
lsp__hz[0] = dlsp_[0];
222
202
for(; i<order; i++)
203
lsp__hz[i] = lsp__hz[i-1] + dlsp[i];
204
205
/* convert back to radians */
206
207
for(i=0; i<order; i++)
223
208
lsp_[i] = (PI/4000.0)*lsp__hz[i];
224
225
//printf("%d dlsp_ %3.2f lsp_ %3.2f\n", i, dlsp_[i], lsp__hz[i]);
226
}
227
228
209
}
229
210
230
#ifdef __EXPERIMENTAL__
231
211
/*---------------------------------------------------------------------------*\
232
212
233
lspvq_quantise
213
lspd_vq_quantise
234
214
235
Vector LSP quantiser.
215
Vector lsp difference quantiser.
236
216
237
217
\*---------------------------------------------------------------------------*/
238
218
239
void lspvq_quantise(
219
void lspdvq_quantise(
240
220
float lsp[],
241
221
float lsp_[],
242
222
int order
243
223
)
244
224
{
245
225
int i,k,m,ncb, nlsp;
246
float wt[LPC_ORD], lsp_hz[LPC_ORD];
247
const float *cb;
248
float se;
249
int index;
250
251
for(i=0; i<LPC_ORD; i++) {
252
wt[i] = 1.0;
253
lsp_hz[i] = 4000.0*lsp[i]/PI;
254
}
255
256
/* scalar quantise LSPs 1,2,3,4 */
257
258
/* simple uniform scalar quantisers */
259
260
for(i=0; i<4; i++) {
261
k = lsp_cb[i].k;
262
m = lsp_cb[i].m;
263
cb = lsp_cb[i].cb;
264
index = quantise(cb, &lsp_hz[i], wt, k, m, &se);
265
lsp_[i] = cb[index*k]*PI/4000.0;
266
}
267
268
//#define WGHT
269
#ifdef WGHT
270
for(i=4; i<9; i++) {
271
wt[i] = 1.0/(lsp[i]-lsp[i-1]) + 1.0/(lsp[i+1]-lsp[i]);
272
//printf("wt[%d] = %f\n", i, wt[i]);
273
}
274
wt[9] = 1.0/(lsp[i]-lsp[i-1]);
275
#endif
276
277
/* VQ LSPs 5,6,7,8,9,10 */
278
279
ncb = 4;
280
nlsp = 4;
281
k = lsp_cbjnd[ncb].k;
282
m = lsp_cbjnd[ncb].m;
283
cb = lsp_cbjnd[ncb].cb;
284
index = quantise(cb, &lsp_hz[nlsp], &wt[nlsp], k, m, &se);
285
for(i=4; i<LPC_ORD; i++) {
286
lsp_[i] = cb[index*k+i-4]*(PI/4000.0);
287
//printf("%4.f (%4.f) ", lsp_hz[i], cb[index*k+i-4]);
288
}
289
}
290
291
/*---------------------------------------------------------------------------*\
292
293
lspjnd_quantise
294
295
Experimental JND LSP quantiser.
296
297
\*---------------------------------------------------------------------------*/
298
299
void lspjnd_quantise(float lsps[], float lsps_[], int order)
300
{
301
int i,k,m;
302
float wt[LPC_ORD], lsps_hz[LPC_ORD];
303
const float *cb;
304
float se = 0.0;
305
int index;
306
307
for(i=0; i<LPC_ORD; i++) {
308
wt[i] = 1.0;
309
}
310
311
/* convert to Hz */
312
313
for(i=0; i<LPC_ORD; i++) {
314
lsps_hz[i] = lsps[i]*(4000.0/PI);
315
lsps_[i] = lsps[i];
316
}
317
318
/* simple uniform scalar quantisers */
319
320
for(i=0; i<4; i++) {
321
k = lsp_cbjnd[i].k;
322
m = lsp_cbjnd[i].m;
323
cb = lsp_cbjnd[i].cb;
324
index = quantise(cb, &lsps_hz[i], wt, k, m, &se);
325
lsps_[i] = cb[index*k]*(PI/4000.0);
326
}
327
328
/* VQ LSPs 5,6,7,8,9,10 */
329
330
k = lsp_cbjnd[4].k;
331
m = lsp_cbjnd[4].m;
332
cb = lsp_cbjnd[4].cb;
333
index = quantise(cb, &lsps_hz[4], &wt[4], k, m, &se);
334
//printf("k = %d m = %d c[0] %f cb[k] %f\n", k,m,cb[0],cb[k]);
335
//printf("index = %4d: ", index);
336
for(i=4; i<LPC_ORD; i++) {
337
lsps_[i] = cb[index*k+i-4]*(PI/4000.0);
338
//printf("%4.f (%4.f) ", lsps_hz[i], cb[index*k+i-4]);
339
}
340
//printf("\n");
341
}
342
343
void compute_weights(const float *x, float *w, int ndim);
344
345
/*---------------------------------------------------------------------------*\
346
347
lspdt_quantise
348
349
LSP difference in time quantiser. Split VQ, encoding LSPs 1-4 with
350
one VQ, and LSPs 5-10 with a second. Update of previous lsp memory
351
is done outside of this function to handle dT between 10 or 20ms
352
frames.
353
354
mode action
355
------------------
356
357
LSPDT_ALL VQ LSPs 1-4 and 5-10
358
LSPDT_LOW Just VQ LSPs 1-4, for LSPs 5-10 just copy previous
359
LSPDT_HIGH Just VQ LSPs 5-10, for LSPs 1-4 just copy previous
360
361
\*---------------------------------------------------------------------------*/
362
363
void lspdt_quantise(float lsps[], float lsps_[], float lsps__prev[], int mode)
364
{
365
int i;
366
float wt[LPC_ORD];
367
float lsps_dt[LPC_ORD];
368
#ifdef TRY_LSPDT_VQ
369
int k,m;
370
int index;
371
const float *cb;
372
float se = 0.0;
373
#endif // TRY_LSPDT_VQ
374
375
//compute_weights(lsps, wt, LPC_ORD);
376
for(i=0; i<LPC_ORD; i++) {
377
wt[i] = 1.0;
378
}
379
380
//compute_weights(lsps, wt, LPC_ORD );
381
382
for(i=0; i<LPC_ORD; i++) {
383
lsps_dt[i] = lsps[i] - lsps__prev[i];
384
lsps_[i] = lsps__prev[i];
385
}
386
387
//#define TRY_LSPDT_VQ
388
#ifdef TRY_LSPDT_VQ
389
/* this actually improves speech a bit, but 40ms updates works surprsingly well.... */
390
k = lsp_cbdt[0].k;
391
m = lsp_cbdt[0].m;
392
cb = lsp_cbdt[0].cb;
393
index = quantise(cb, lsps_dt, wt, k, m, &se);
394
for(i=0; i<LPC_ORD; i++) {
395
lsps_[i] += cb[index*k + i];
396
}
397
#endif
398
399
}
400
#endif
401
402
#define MIN(a,b) ((a)<(b)?(a):(b))
403
#define MAX_ENTRIES 16384
404
405
void compute_weights(const float *x, float *w, int ndim)
406
{
407
int i;
408
w[0] = MIN(x[0], x[1]-x[0]);
409
for (i=1;i<ndim-1;i++)
410
w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
411
w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], PI-x[ndim-1]);
412
413
for (i=0;i<ndim;i++)
414
w[i] = 1./(.01+w[i]);
415
//w[0]*=3;
416
//w[1]*=2;
417
}
418
419
/* LSP weight calculation ported from m-file function kindly submitted
420
by Anssi, OH3GDD */
421
422
void compute_weights_anssi_mode2(const float *x, float *w, int ndim)
423
{
424
int i;
425
float d[LPC_ORD];
426
427
assert(ndim == LPC_ORD);
428
429
for(i=0; i<LPC_ORD; i++)
430
d[i] = 1.0;
431
432
d[0] = x[1];
433
for (i=1; i<LPC_ORD-1; i++)
434
d[i] = x[i+1] - x[i-1];
435
d[LPC_ORD-1] = PI - x[8];
436
for (i=0; i<LPC_ORD; i++) {
437
if (x[i]<((400.0/4000.0)*PI))
438
w[i]=5.0/(0.01+d[i]);
439
else if (x[i]<((700.0/4000.0)*PI))
440
w[i]=4.0/(0.01+d[i]);
441
else if (x[i]<((1200.0/4000.0)*PI))
442
w[i]=3.0/(0.01+d[i]);
443
else if (x[i]<((2000.0/4000.0)*PI))
444
w[i]=2.0/(0.01+d[i]);
445
else
446
w[i]=1.0/(0.01+d[i]);
447
448
w[i]=pow(w[i]+0.3, 0.66);
449
}
450
}
451
452
int find_nearest(const float *codebook, int nb_entries, float *x, int ndim)
453
{
454
int i, j;
455
float min_dist = 1e15;
456
int nearest = 0;
457
458
for (i=0;i<nb_entries;i++)
459
{
460
float dist=0;
461
for (j=0;j<ndim;j++)
462
dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
463
if (dist<min_dist)
464
{
465
min_dist = dist;
466
nearest = i;
467
}
468
}
469
return nearest;
470
}
471
472
int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim)
473
{
474
int i, j;
475
float min_dist = 1e15;
476
int nearest = 0;
477
478
for (i=0;i<nb_entries;i++)
479
{
480
float dist=0;
481
for (j=0;j<ndim;j++)
482
dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
483
if (dist<min_dist)
484
{
485
min_dist = dist;
486
nearest = i;
487
}
488
}
489
return nearest;
490
}
491
492
void lspjvm_quantise(float *x, float *xq, int ndim)
493
{
494
int i, n1, n2, n3;
495
float err[LPC_ORD], err2[LPC_ORD], err3[LPC_ORD];
496
float w[LPC_ORD], w2[LPC_ORD], w3[LPC_ORD];
497
const float *codebook1 = lsp_cbjvm[0].cb;
498
const float *codebook2 = lsp_cbjvm[1].cb;
499
const float *codebook3 = lsp_cbjvm[2].cb;
500
501
w[0] = MIN(x[0], x[1]-x[0]);
502
for (i=1;i<ndim-1;i++)
503
w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
504
w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], PI-x[ndim-1]);
505
506
compute_weights(x, w, ndim);
507
508
n1 = find_nearest(codebook1, lsp_cbjvm[0].m, x, ndim);
509
510
for (i=0;i<ndim;i++)
511
{
512
xq[i] = codebook1[ndim*n1+i];
513
err[i] = x[i] - xq[i];
514
}
515
for (i=0;i<ndim/2;i++)
516
{
517
err2[i] = err[2*i];
518
err3[i] = err[2*i+1];
519
w2[i] = w[2*i];
520
w3[i] = w[2*i+1];
521
}
522
n2 = find_nearest_weighted(codebook2, lsp_cbjvm[1].m, err2, w2, ndim/2);
523
n3 = find_nearest_weighted(codebook3, lsp_cbjvm[2].m, err3, w3, ndim/2);
524
525
for (i=0;i<ndim/2;i++)
526
{
527
xq[2*i] += codebook2[ndim*n2/2+i];
528
xq[2*i+1] += codebook3[ndim*n3/2+i];
529
}
530
}
531
532
#ifdef __EXPERIMENTAL__
533
534
#define MBEST_STAGES 4
535
536
struct MBEST_LIST {
537
int index[MBEST_STAGES]; /* index of each stage that lead us to this error */
538
float error;
539
};
540
541
struct MBEST {
542
int entries; /* number of entries in mbest list */
543
struct MBEST_LIST *list;
544
};
545
546
547
static struct MBEST *mbest_create(int entries) {
548
int i,j;
549
struct MBEST *mbest;
550
551
assert(entries > 0);
552
mbest = (struct MBEST *)malloc(sizeof(struct MBEST));
553
assert(mbest != NULL);
554
555
mbest->entries = entries;
556
mbest->list = (struct MBEST_LIST *)malloc(entries*sizeof(struct MBEST_LIST));
557
assert(mbest->list != NULL);
558
559
for(i=0; i<mbest->entries; i++) {
560
for(j=0; j<MBEST_STAGES; j++)
561
mbest->list[i].index[j] = 0;
562
mbest->list[i].error = 1E32;
563
}
564
565
return mbest;
566
}
567
568
569
static void mbest_destroy(struct MBEST *mbest) {
570
assert(mbest != NULL);
571
free(mbest->list);
572
free(mbest);
573
}
574
575
576
/*---------------------------------------------------------------------------*\
577
578
mbest_insert
579
580
Insert the results of a vector to codebook entry comparison. The
581
list is ordered in order or error, so those entries with the
582
smallest error will be first on the list.
583
584
\*---------------------------------------------------------------------------*/
585
586
static void mbest_insert(struct MBEST *mbest, int index[], float error) {
587
int i, j, found;
588
struct MBEST_LIST *list = mbest->list;
589
int entries = mbest->entries;
590
591
found = 0;
592
for(i=0; i<entries && !found; i++)
593
if (error < list[i].error) {
594
found = 1;
595
for(j=entries-1; j>i; j--)
596
list[j] = list[j-1];
597
for(j=0; j<MBEST_STAGES; j++)
598
list[i].index[j] = index[j];
599
list[i].error = error;
600
}
601
}
602
603
604
static void mbest_print(char title[], struct MBEST *mbest) {
605
int i,j;
606
607
printf("%s\n", title);
608
for(i=0; i<mbest->entries; i++) {
609
for(j=0; j<MBEST_STAGES; j++)
610
printf(" %4d ", mbest->list[i].index[j]);
611
printf(" %f\n", mbest->list[i].error);
612
}
613
}
614
615
616
/*---------------------------------------------------------------------------*\
617
618
mbest_search
619
620
Searches vec[] to a codebbook of vectors, and maintains a list of the mbest
621
closest matches.
622
623
\*---------------------------------------------------------------------------*/
624
625
static void mbest_search(
626
const float *cb, /* VQ codebook to search */
627
float vec[], /* target vector */
628
float w[], /* weighting vector */
629
int k, /* dimension of vector */
630
int m, /* number on entries in codebook */
631
struct MBEST *mbest, /* list of closest matches */
632
int index[] /* indexes that lead us here */
633
)
634
{
635
float e;
636
int i,j;
637
float diff;
638
639
for(j=0; j<m; j++) {
640
e = 0.0;
641
for(i=0; i<k; i++) {
642
diff = cb[j*k+i]-vec[i];
643
e += pow(diff*w[i],2.0);
644
}
645
index[0] = j;
646
mbest_insert(mbest, index, e);
647
}
648
}
649
650
651
/* 3 stage VQ LSP quantiser. Design and guidance kindly submitted by Anssi, OH3GDD */
652
653
void lspanssi_quantise(float *x, float *xq, int ndim, int mbest_entries)
654
{
655
int i, j, n1, n2, n3, n4;
656
float w[LPC_ORD];
657
const float *codebook1 = lsp_cbvqanssi[0].cb;
658
const float *codebook2 = lsp_cbvqanssi[1].cb;
659
const float *codebook3 = lsp_cbvqanssi[2].cb;
660
const float *codebook4 = lsp_cbvqanssi[3].cb;
661
struct MBEST *mbest_stage1, *mbest_stage2, *mbest_stage3, *mbest_stage4;
662
float target[LPC_ORD];
663
int index[MBEST_STAGES];
664
665
mbest_stage1 = mbest_create(mbest_entries);
666
mbest_stage2 = mbest_create(mbest_entries);
667
mbest_stage3 = mbest_create(mbest_entries);
668
mbest_stage4 = mbest_create(mbest_entries);
669
for(i=0; i<MBEST_STAGES; i++)
670
index[i] = 0;
671
672
compute_weights_anssi_mode2(x, w, ndim);
673
674
#ifdef DUMP
675
dump_weights(w, ndim);
676
#endif
677
678
/* Stage 1 */
679
680
mbest_search(codebook1, x, w, ndim, lsp_cbvqanssi[0].m, mbest_stage1, index);
681
mbest_print("Stage 1:", mbest_stage1);
682
683
/* Stage 2 */
684
685
for (j=0; j<mbest_entries; j++) {
686
index[1] = n1 = mbest_stage1->list[j].index[0];
687
for(i=0; i<ndim; i++)
688
target[i] = x[i] - codebook1[ndim*n1+i];
689
mbest_search(codebook2, target, w, ndim, lsp_cbvqanssi[1].m, mbest_stage2, index);
690
}
691
mbest_print("Stage 2:", mbest_stage2);
692
693
/* Stage 3 */
694
695
for (j=0; j<mbest_entries; j++) {
696
index[2] = n1 = mbest_stage2->list[j].index[1];
697
index[1] = n2 = mbest_stage2->list[j].index[0];
698
for(i=0; i<ndim; i++)
699
target[i] = x[i] - codebook1[ndim*n1+i] - codebook2[ndim*n2+i];
700
mbest_search(codebook3, target, w, ndim, lsp_cbvqanssi[2].m, mbest_stage3, index);
701
}
702
mbest_print("Stage 3:", mbest_stage3);
703
704
/* Stage 4 */
705
706
for (j=0; j<mbest_entries; j++) {
707
index[3] = n1 = mbest_stage3->list[j].index[2];
708
index[2] = n2 = mbest_stage3->list[j].index[1];
709
index[1] = n3 = mbest_stage3->list[j].index[0];
710
for(i=0; i<ndim; i++)
711
target[i] = x[i] - codebook1[ndim*n1+i] - codebook2[ndim*n2+i] - codebook3[ndim*n3+i];
712
mbest_search(codebook4, target, w, ndim, lsp_cbvqanssi[3].m, mbest_stage4, index);
713
}
714
mbest_print("Stage 4:", mbest_stage4);
715
716
n1 = mbest_stage4->list[0].index[3];
717
n2 = mbest_stage4->list[0].index[2];
718
n3 = mbest_stage4->list[0].index[1];
719
n4 = mbest_stage4->list[0].index[0];
720
for (i=0;i<ndim;i++)
721
xq[i] = codebook1[ndim*n1+i] + codebook2[ndim*n2+i] + codebook3[ndim*n3+i] + codebook4[ndim*n4+i];
722
723
mbest_destroy(mbest_stage1);
724
mbest_destroy(mbest_stage2);
725
mbest_destroy(mbest_stage3);
726
mbest_destroy(mbest_stage4);
727
}
728
#endif
729
730
int check_lsp_order(float lsp[], int lpc_order)
226
float dlsp[LPC_MAX];
227
float dlsp_[LPC_MAX];
228
float wt[LPC_ORD];
229
const float *cb;
230
float se = 0.0;
231
int index;
232
233
dlsp[0] = lsp[0];
234
for(i=1; i<order; i++)
235
dlsp[i] = lsp[i] - lsp[i-1];
236
237
for(i=0; i<order; i++)
238
dlsp_[i] = dlsp[i];
239
240
for(i=0; i<order; i++)
241
wt[i] = 1.0;
242
243
/* scalar quantise dLSPs 1,2,3,4,5 */
244
245
for(i=0; i<5; i++) {
246
if (i)
247
dlsp[i] = (lsp[i] - lsp_[i-1])*4000.0/PI;
248
else
249
dlsp[0] = lsp[0]*4000.0/PI;
250
251
k = lsp_cbdvq[i].k;
252
m = lsp_cbdvq[i].m;
253
cb = lsp_cbdvq[i].cb;
254
index = quantise(cb, &dlsp[i], wt, k, m, &se);
255
dlsp_[i] = cb[index*k]*PI/4000.0;
256
257
if (i)
258
lsp_[i] = lsp_[i-1] + dlsp_[i];
259
else
260
lsp_[0] = dlsp_[0];
261
}
262
dlsp[i] = lsp[i] - lsp_[i-1];
263
dlsp_[i] = dlsp[i];
264
265
//printf("lsp[0] %f lsp_[0] %f\n", lsp[0], lsp_[0]);
266
//printf("lsp[1] %f lsp_[1] %f\n", lsp[1], lsp_[1]);
267
268
#ifdef TT
269
/* VQ dLSPs 3,4,5 */
270
271
ncb = 2;
272
nlsp = 2;
273
k = lsp_cbdvq[ncb].k;
274
m = lsp_cbdvq[ncb].m;
275
cb = lsp_cbdvq[ncb].cb;
276
index = quantise(cb, &dlsp[nlsp], wt, k, m, &se);
277
dlsp_[nlsp] = cb[index*k];
278
dlsp_[nlsp+1] = cb[index*k+1];
279
dlsp_[nlsp+2] = cb[index*k+2];
280
281
lsp_[0] = dlsp_[0];
282
for(i=1; i<5; i++)
283
lsp_[i] = lsp_[i-1] + dlsp_[i];
284
dlsp[i] = lsp[i] - lsp_[i-1];
285
dlsp_[i] = dlsp[i];
286
#endif
287
/* VQ dLSPs 6,7,8,9,10 */
288
289
ncb = 5;
290
nlsp = 5;
291
k = lsp_cbdvq[ncb].k;
292
m = lsp_cbdvq[ncb].m;
293
cb = lsp_cbdvq[ncb].cb;
294
index = quantise(cb, &dlsp[nlsp], wt, k, m, &se);
295
dlsp_[nlsp] = cb[index*k];
296
dlsp_[nlsp+1] = cb[index*k+1];
297
dlsp_[nlsp+2] = cb[index*k+2];
298
dlsp_[nlsp+3] = cb[index*k+3];
299
dlsp_[nlsp+4] = cb[index*k+4];
300
301
/* rebuild LSPs for dLSPs */
302
303
lsp_[0] = dlsp_[0];
304
for(i=1; i<order; i++)
305
lsp_[i] = lsp_[i-1] + dlsp_[i];
306
}
307
308
void check_lsp_order(float lsp[], int lpc_order)
731
309
{
732
310
int i;
733
311
float tmp;
734
int swaps = 0;
735
312
736
313
for(i=1; i<lpc_order; i++)
737
314
if (lsp[i] < lsp[i-1]) {
738
//fprintf(stderr, "swap %d\n",i);
739
swaps++;
315
printf("swap %d\n",i);
740
316
tmp = lsp[i-1];
741
lsp[i-1] = lsp[i]-0.1;
742
lsp[i] = tmp+0.1;
743
i = 1; /* start check again, as swap may have caused out of order */
317
lsp[i-1] = lsp[i]-0.05;
318
lsp[i] = tmp+0.05;
744
319
}
745
746
return swaps;
747
320
}
748
321
749
322
void force_min_lsp_dist(float lsp[], int lpc_order)
756
329
}
757
330
}
758
331
759
760
332
/*---------------------------------------------------------------------------*\
761
333
762
lpc_post_filter()
763
764
Applies a post filter to the LPC synthesis filter power spectrum
765
Pw, which supresses the inter-formant energy.
766
767
The algorithm is from p267 (Section 8.6) of "Digital Speech",
768
edited by A.M. Kondoz, 1994 published by Wiley and Sons. Chapter 8
769
of this text is on the MBE vocoder, and this is a freq domain
770
adaptation of post filtering commonly used in CELP.
771
772
I used the Octave simulation lpcpf.m to get an understaing of the
773
algorithm.
774
775
Requires two more FFTs which is significantly more MIPs. However
776
it should be possible to implement this more efficiently in the
777
time domain. Just not sure how to handle relative time delays
778
between the synthesis stage and updating these coeffs. A smaller
779
FFT size might also be accetable to save CPU.
780
781
TODO:
782
[ ] sync var names between Octave and C version
783
[ ] doc gain normalisation
784
[ ] I think the first FFT is not rqd as we do the same
785
thing in aks_to_M2().
334
lpc_model_amplitudes
335
336
Derive a LPC model for amplitude samples then estimate amplitude samples
337
from this model with optional LSP quantisation.
338
339
Returns the spectral distortion for this frame.
786
340
787
341
\*---------------------------------------------------------------------------*/
788
342
789
void lpc_post_filter(kiss_fft_cfg fft_fwd_cfg, MODEL *model, COMP Pw[], float ak[],
790
int order, int dump, float beta, float gamma, int bass_boost)
343
float lpc_model_amplitudes(
344
float Sn[], /* Input frame of speech samples */
345
float w[],
346
MODEL *model, /* sinusoidal model parameters */
347
int order, /* LPC model order */
348
int lsp_quant, /* optional LSP quantisation if non-zero */
349
float ak[] /* output aks */
350
)
791
351
{
792
int i;
793
COMP x[FFT_ENC]; /* input to FFTs */
794
COMP Aw[FFT_ENC]; /* LPC analysis filter spectrum */
795
COMP Ww[FFT_ENC]; /* weighting spectrum */
796
float Rw[FFT_ENC]; /* R = WA */
797
float e_before, e_after, gain;
798
float Pfw[FFT_ENC]; /* Post filter mag spectrum */
799
float max_Rw, min_Rw;
800
float coeff;
801
TIMER_VAR(tstart, tfft1, taw, tfft2, tww, tr);
802
803
TIMER_SAMPLE(tstart);
804
805
/* Determine LPC inverse filter spectrum 1/A(exp(jw)) -----------*/
806
807
/* we actually want the synthesis filter A(exp(jw)) but the
808
inverse (analysis) filter is easier to find as it's FIR, we
809
just use the inverse of 1/A to get the synthesis filter
810
A(exp(jw)) */
811
812
for(i=0; i<FFT_ENC; i++) {
813
x[i].real = 0.0;
814
x[i].imag = 0.0;
815
}
816
817
for(i=0; i<=order; i++)
818
x[i].real = ak[i];
819
kiss_fft(fft_fwd_cfg, (kiss_fft_cpx *)x, (kiss_fft_cpx *)Aw);
820
821
TIMER_SAMPLE_AND_LOG(tfft1, tstart, " fft1");
822
823
for(i=0; i<FFT_ENC/2; i++) {
824
Aw[i].real = 1.0/(Aw[i].real*Aw[i].real + Aw[i].imag*Aw[i].imag);
825
}
826
827
TIMER_SAMPLE_AND_LOG(taw, tfft1, " Aw");
828
829
/* Determine weighting filter spectrum W(exp(jw)) ---------------*/
830
831
for(i=0; i<FFT_ENC; i++) {
832
x[i].real = 0.0;
833
x[i].imag = 0.0;
834
}
835
836
x[0].real = ak[0];
837
coeff = gamma;
838
for(i=1; i<=order; i++) {
839
x[i].real = ak[i] * coeff;
840
coeff *= gamma;
841
}
842
kiss_fft(fft_fwd_cfg, (kiss_fft_cpx *)x, (kiss_fft_cpx *)Ww);
843
844
TIMER_SAMPLE_AND_LOG(tfft2, taw, " fft2");
845
846
for(i=0; i<FFT_ENC/2; i++) {
847
Ww[i].real = Ww[i].real*Ww[i].real + Ww[i].imag*Ww[i].imag;
848
}
849
850
TIMER_SAMPLE_AND_LOG(tww, tfft2, " Ww");
851
852
/* Determined combined filter R = WA ---------------------------*/
853
854
max_Rw = 0.0; min_Rw = 1E32;
855
for(i=0; i<FFT_ENC/2; i++) {
856
Rw[i] = sqrtf(Ww[i].real * Aw[i].real);
857
if (Rw[i] > max_Rw)
858
max_Rw = Rw[i];
859
if (Rw[i] < min_Rw)
860
min_Rw = Rw[i];
861
862
}
863
864
TIMER_SAMPLE_AND_LOG(tr, tww, " R");
865
866
#ifdef DUMP
867
if (dump)
868
dump_Rw(Rw);
869
#endif
870
871
/* create post filter mag spectrum and apply ------------------*/
872
873
/* measure energy before post filtering */
874
875
e_before = 1E-4;
876
for(i=0; i<FFT_ENC/2; i++)
877
e_before += Pw[i].real;
878
879
/* apply post filter and measure energy */
880
881
#ifdef DUMP
882
if (dump)
883
dump_Pwb(Pw);
884
#endif
885
886
e_after = 1E-4;
887
for(i=0; i<FFT_ENC/2; i++) {
888
Pfw[i] = powf(Rw[i], beta);
889
Pw[i].real *= Pfw[i] * Pfw[i];
890
e_after += Pw[i].real;
891
}
892
gain = e_before/e_after;
893
894
/* apply gain factor to normalise energy */
895
896
for(i=0; i<FFT_ENC/2; i++) {
897
Pw[i].real *= gain;
898
}
899
900
if (bass_boost) {
901
/* add 3dB to first 1 kHz to account for LP effect of PF */
902
903
for(i=0; i<FFT_ENC/8; i++) {
904
Pw[i].real *= 1.4*1.4;
905
}
906
}
907
908
TIMER_SAMPLE_AND_LOG2(tr, " filt");
352
float Wn[M];
353
float R[LPC_MAX+1];
354
float E;
355
int i,j;
356
float snr;
357
float lsp[LPC_MAX];
358
float lsp_hz[LPC_MAX];
359
float lsp_[LPC_MAX];
360
int roots; /* number of LSP roots found */
361
int index;
362
float se = 0.0;
363
int k,m;
364
const float * cb;
365
float wt[LPC_MAX];
366
367
for(i=0; i<M; i++)
368
Wn[i] = Sn[i]*w[i];
369
autocorrelate(Wn,R,M,order);
370
levinson_durbin(R,ak,order);
371
372
E = 0.0;
373
for(i=0; i<=order; i++)
374
E += ak[i]*R[i];
375
376
for(i=0; i<order; i++)
377
wt[i] = 1.0;
378
379
if (lsp_quant) {
380
roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1);
381
if (roots != order)
382
printf("LSP roots not found\n");
383
384
/* convert from radians to Hz to make quantisers more
385
human readable */
386
387
for(i=0; i<order; i++)
388
lsp_hz[i] = (4000.0/PI)*lsp[i];
389
390
/* simple uniform scalar quantisers */
391
392
for(i=0; i<10; i++) {
393
k = lsp_cb[i].k;
394
m = lsp_cb[i].m;
395
cb = lsp_cb[i].cb;
396
index = quantise(cb, &lsp_hz[i], wt, k, m, &se);
397
lsp_hz[i] = cb[index*k];
398
}
399
400
/* experiment: simulating uniform quantisation error
401
for(i=0; i<order; i++)
402
lsp[i] += PI*(12.5/4000.0)*(1.0 - 2.0*(float)rand()/RAND_MAX);
403
*/
404
405
for(i=0; i<order; i++)
406
lsp[i] = (PI/4000.0)*lsp_hz[i];
407
408
/* Bandwidth Expansion (BW). Prevents any two LSPs getting too
409
close together after quantisation. We know from experiment
410
that LSP quantisation errors < 12.5Hz (25Hz setp size) are
411
inaudible so we use that as the minimum LSP separation.
412
*/
413
414
for(i=1; i<5; i++) {
415
if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0))
416
lsp[i] = lsp[i-1] + PI*(12.5/4000.0);
417
}
418
419
/* as quantiser gaps increased, larger BW expansion was required
420
to prevent twinkly noises */
421
422
for(i=5; i<8; i++) {
423
if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0))
424
lsp[i] = lsp[i-1] + PI*(25.0/4000.0);
425
}
426
for(i=8; i<order; i++) {
427
if (lsp[i] - lsp[i-1] < PI*(75.0/4000.0))
428
lsp[i] = lsp[i-1] + PI*(75.0/4000.0);
429
}
430
431
for(j=0; j<order; j++)
432
lsp_[j] = lsp[j];
433
434
lsp_to_lpc(lsp_, ak, order);
435
#ifdef DUMP
436
dump_lsp(lsp);
437
#endif
438
}
439
440
#ifdef DUMP
441
dump_E(E);
442
#endif
443
#ifdef SIM_QUANT
444
/* simulated LPC energy quantisation */
445
{
446
float e = 10.0*log10(E);
447
e += 2.0*(1.0 - 2.0*(float)rand()/RAND_MAX);
448
E = pow(10.0,e/10.0);
449
}
450
#endif
451
452
aks_to_M2(ak,order,model,E,&snr, 1); /* {ak} -> {Am} LPC decode */
453
454
return snr;
909
455
}
910
456
911
912
457
/*---------------------------------------------------------------------------*\
913
458
914
459
aks_to_M2()
920
465
\*---------------------------------------------------------------------------*/
921
466
922
467
void aks_to_M2(
923
kiss_fft_cfg fft_fwd_cfg,
924
float ak[], /* LPC's */
925
int order,
926
MODEL *model, /* sinusoidal model parameters for this frame */
927
float E, /* energy term */
928
float *snr, /* signal to noise ratio for this frame in dB */
929
int dump, /* true to dump sample to dump file */
930
int sim_pf, /* true to simulate a post filter */
931
int pf, /* true to LPC post filter */
932
int bass_boost, /* enable LPC filter 0-1khz 3dB boost */
933
float beta,
934
float gamma /* LPC post filter parameters */
468
float ak[], /* LPC's */
469
int order,
470
MODEL *model, /* sinusoidal model parameters for this frame */
471
float E, /* energy term */
472
float *snr, /* signal to noise ratio for this frame in dB */
473
int dump /* true to dump sample to dump file */
935
474
)
936
475
{
937
COMP pw[FFT_ENC]; /* input to FFT for power spectrum */
938
COMP Pw[FFT_ENC]; /* output power spectrum */
476
COMP Pw[FFT_DEC]; /* power spectrum */
939
477
int i,m; /* loop variables */
940
478
int am,bm; /* limits of current band */
941
479
float r; /* no. rads/bin */
942
480
float Em; /* energy in band */
943
481
float Am; /* spectral amplitude sample */
944
482
float signal, noise;
945
TIMER_VAR(tstart, tfft, tpw, tpf);
946
947
TIMER_SAMPLE(tstart);
948
949
r = TWO_PI/(FFT_ENC);
483
484
r = TWO_PI/(FFT_DEC);
950
485
951
486
/* Determine DFT of A(exp(jw)) --------------------------------------------*/
952
487
953
for(i=0; i<FFT_ENC; i++) {
954
pw[i].real = 0.0;
955
pw[i].imag = 0.0;
488
for(i=0; i<FFT_DEC; i++) {
489
Pw[i].real = 0.0;
490
Pw[i].imag = 0.0;
956
491
}
957
492
958
493
for(i=0; i<=order; i++)
959
pw[i].real = ak[i];
960
kiss_fft(fft_fwd_cfg, (kiss_fft_cpx *)pw, (kiss_fft_cpx *)Pw);
961
962
TIMER_SAMPLE_AND_LOG(tfft, tstart, " fft");
494
Pw[i].real = ak[i];
495
fft(&Pw[0].real,FFT_DEC,1);
963
496
964
497
/* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/
965
498
966
for(i=0; i<FFT_ENC/2; i++)
499
for(i=0; i<FFT_DEC/2; i++)
967
500
Pw[i].real = E/(Pw[i].real*Pw[i].real + Pw[i].imag*Pw[i].imag);
968
969
TIMER_SAMPLE_AND_LOG(tpw, tfft, " Pw");
970
971
if (pf)
972
lpc_post_filter(fft_fwd_cfg, model, Pw, ak, order, dump, beta, gamma, bass_boost);
973
974
TIMER_SAMPLE_AND_LOG(tpf, tpw, " LPC post filter");
975
976
#ifdef DUMP
501
#ifdef DUMP
977
502
if (dump)
978
503
dump_Pw(Pw);
979
#endif
980
981
/* Determine magnitudes from P(w) ----------------------------------------*/
982
983
/* when used just by decoder {A} might be all zeroes so init signal
984
and noise to prevent log(0) errors */
985
986
signal = 1E-30; noise = 1E-32;
987
504
#endif
505
506
/* Determine magnitudes by linear interpolation of P(w) -------------------*/
507
508
signal = noise = 0.0;
988
509
for(m=1; m<=model->L; m++) {
989
am = (int)((m - 0.5)*model->Wo/r + 0.5);
990
bm = (int)((m + 0.5)*model->Wo/r + 0.5);
991
Em = 0.0;
992
993
for(i=am; i<bm; i++)
994
Em += Pw[i].real;
995
Am = sqrtf(Em);
996
997
signal += model->A[m]*model->A[m];
998
noise += (model->A[m] - Am)*(model->A[m] - Am);
999
1000
/* This code significantly improves perf of LPC model, in
1001
particular when combined with phase0. The LPC spectrum tends
1002
to track just under the peaks of the spectral envelope, and
1003
just above nulls. This algorithm does the reverse to
1004
compensate - raising the amplitudes of spectral peaks, while
1005
attenuating the null. This enhances the formants, and
1006
supresses the energy between formants. */
1007
1008
if (sim_pf) {
1009
if (Am > model->A[m])
1010
Am *= 0.7;
1011
if (Am < model->A[m])
1012
Am *= 1.4;
1013
}
1014
1015
model->A[m] = Am;
510
am = floor((m - 0.5)*model->Wo/r + 0.5);
511
bm = floor((m + 0.5)*model->Wo/r + 0.5);
512
Em = 0.0;
513
514
for(i=am; i<bm; i++)
515
Em += Pw[i].real;
516
Am = sqrt(Em);
517
518
signal += pow(model->A[m],2.0);
519
noise += pow(model->A[m] - Am,2.0);
520
model->A[m] = Am;
1016
521
}
1017
*snr = 10.0*log10f(signal/noise);
1018
1019
TIMER_SAMPLE_AND_LOG2(tpf, " rec");
522
*snr = 10.0*log10(signal/noise);
1020
523
}
1021
524
1022
525
/*---------------------------------------------------------------------------*\
1037
540
float norm;
1038
541
1039
542
norm = (Wo - Wo_min)/(Wo_max - Wo_min);
1040
index = floorf(WO_LEVELS * norm + 0.5);
543
index = floor(WO_LEVELS * norm + 0.5);
1041
544
if (index < 0 ) index = 0;
1042
545
if (index > (WO_LEVELS-1)) index = WO_LEVELS-1;
1043
546
1069
572
1070
573
/*---------------------------------------------------------------------------*\
1071
574
1072
FUNCTION....: encode_Wo_dt()
1073
AUTHOR......: David Rowe
1074
DATE CREATED: 6 Nov 2011
1075
1076
Encodes Wo difference from last frame.
1077
1078
\*---------------------------------------------------------------------------*/
1079
1080
int encode_Wo_dt(float Wo, float prev_Wo)
1081
{
1082
int index, mask, max_index, min_index;
1083
float Wo_min = TWO_PI/P_MAX;
1084
float Wo_max = TWO_PI/P_MIN;
1085
float norm;
1086
1087
norm = (Wo - prev_Wo)/(Wo_max - Wo_min);
1088
index = floor(WO_LEVELS * norm + 0.5);
1089
//printf("ENC index: %d ", index);
1090
1091
/* hard limit */
1092
1093
max_index = (1 << (WO_DT_BITS-1)) - 1;
1094
min_index = - (max_index+1);
1095
if (index > max_index) index = max_index;
1096
if (index < min_index) index = min_index;
1097
//printf("max_index: %d min_index: %d hard index: %d ",
1098
// max_index, min_index, index);
1099
1100
/* mask so that only LSB WO_DT_BITS remain, bit WO_DT_BITS is the sign bit */
1101
1102
mask = ((1 << WO_DT_BITS) - 1);
1103
index &= mask;
1104
//printf("mask: 0x%x index: 0x%x\n", mask, index);
1105
1106
return index;
1107
}
1108
1109
/*---------------------------------------------------------------------------*\
1110
1111
FUNCTION....: decode_Wo_dt()
1112
AUTHOR......: David Rowe
1113
DATE CREATED: 6 Nov 2011
1114
1115
Decodes Wo using WO_DT_BITS difference from last frame.
1116
1117
\*---------------------------------------------------------------------------*/
1118
1119
float decode_Wo_dt(int index, float prev_Wo)
1120
{
1121
float Wo_min = TWO_PI/P_MAX;
1122
float Wo_max = TWO_PI/P_MIN;
1123
float step;
1124
float Wo;
1125
int mask;
1126
1127
/* sign extend index */
1128
1129
//printf("DEC index: %d ");
1130
if (index & (1 << (WO_DT_BITS-1))) {
1131
mask = ~((1 << WO_DT_BITS) - 1);
1132
index |= mask;
1133
}
1134
//printf("DEC mask: 0x%x index: %d \n", mask, index);
1135
1136
step = (Wo_max - Wo_min)/WO_LEVELS;
1137
Wo = prev_Wo + step*(index);
1138
1139
/* bit errors can make us go out of range leading to all sorts of
1140
probs like seg faults */
1141
1142
if (Wo > Wo_max) Wo = Wo_max;
1143
if (Wo < Wo_min) Wo = Wo_min;
1144
1145
return Wo;
1146
}
1147
1148
/*---------------------------------------------------------------------------*\
1149
1150
575
FUNCTION....: speech_to_uq_lsps()
1151
576
AUTHOR......: David Rowe
1152
577
DATE CREATED: 22/8/2010
1167
592
int i, roots;
1168
593
float Wn[M];
1169
594
float R[LPC_MAX+1];
1170
float e, E;
595
float E;
1171
596
1172
e = 0.0;
1173
for(i=0; i<M; i++) {
597
for(i=0; i<M; i++)
1174
598
Wn[i] = Sn[i]*w[i];
1175
e += Wn[i]*Wn[i];
1176
}
1177
1178
/* trap 0 energy case as LPC analysis will fail */
1179
1180
if (e == 0.0) {
1181
for(i=0; i<order; i++)
1182
lsp[i] = (PI/order)*(float)i;
1183
return 0.0;
1184
}
1185
1186
599
autocorrelate(Wn, R, M, order);
1187
600
levinson_durbin(R, ak, order);
1188
601
1190
603
for(i=0; i<=order; i++)
1191
604
E += ak[i]*R[i];
1192
605
1193
/* 15 Hz BW expansion as I can't hear the difference and it may help
1194
help occasional fails in the LSP root finding. Important to do this
1195
after energy calculation to avoid -ve energy values.
1196
*/
1197
1198
for(i=0; i<=order; i++)
1199
ak[i] *= powf(0.994,(float)i);
1200
1201
606
roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1);
1202
607
if (roots != order) {
1203
/* if root finding fails use some benign LSP values instead */
608
/* for some reason LSP roots could not be found */
609
/* some alpha testers are reporting this condition */
610
fprintf(stderr, "LSP roots not found!\nroots = %d\n", roots);
611
for(i=0; i<=order; i++)
612
fprintf(stderr, "a[%d] = %f\n", i, ak[i]);
613
614
/* some benign LSP values we can use instead */
1204
615
for(i=0; i<order; i++)
1205
616
lsp[i] = (PI/order)*(float)i;
1206
617
}
1210
621
1211
622
/*---------------------------------------------------------------------------*\
1212
623
1213
FUNCTION....: encode_lsps_scalar()
624
FUNCTION....: encode_lsps()
1214
625
AUTHOR......: David Rowe
1215
626
DATE CREATED: 22/8/2010
1216
627
1217
Thirty-six bit sclar LSP quantiser. From a vector of unquantised
1218
(floating point) LSPs finds the quantised LSP indexes.
628
From a vector of unquantised (floating point) LSPs finds the quantised
629
LSP indexes.
1219
630
1220
631
\*---------------------------------------------------------------------------*/
1221
632
1222
void encode_lsps_scalar(int indexes[], float lsp[], int order)
633
void encode_lsps(int indexes[], float lsp[], int order)
1223
634
{
1224
635
int i,k,m;
1225
636
float wt[1];
1226
637
float lsp_hz[LPC_MAX];
1227
638
const float * cb;
1228
float se;
639
float se = 0.0;
1229
640
1230
641
/* convert from radians to Hz so we can use human readable
1231
642
frequencies */
1233
644
for(i=0; i<order; i++)
1234
645
lsp_hz[i] = (4000.0/PI)*lsp[i];
1235
646
1236
/* scalar quantisers */
647
/* simple uniform scalar quantisers */
1237
648
1238
649
wt[0] = 1.0;
1239
650
for(i=0; i<order; i++) {
1246
657
1247
658
/*---------------------------------------------------------------------------*\
1248
659
1249
FUNCTION....: decode_lsps_scalar()
660
FUNCTION....: decode_lsps()
1250
661
AUTHOR......: David Rowe
1251
662
DATE CREATED: 22/8/2010
1252
663
1255
666
1256
667
\*---------------------------------------------------------------------------*/
1257
668
1258
void decode_lsps_scalar(float lsp[], int indexes[], int order)
669
void decode_lsps(float lsp[], int indexes[], int order)
1259
670
{
1260
671
int i,k;
1261
672
float lsp_hz[LPC_MAX];
1273
684
lsp[i] = (PI/4000.0)*lsp_hz[i];
1274
685
}
1275
686
1276
1277
#ifdef __EXPERIMENTAL__
1278
1279
/*---------------------------------------------------------------------------*\
1280
1281
FUNCTION....: encode_lsps_diff_freq_vq()
1282
AUTHOR......: David Rowe
1283
DATE CREATED: 15 November 2011
1284
1285
Twenty-five bit LSP quantiser. LSPs 1-4 are quantised with scalar
1286
LSP differences (in frequency, i.e difference from the previous
1287
LSP). LSPs 5-10 are quantised with a VQ trained generated using
1288
vqtrainjnd.c
1289
1290
\*---------------------------------------------------------------------------*/
1291
1292
void encode_lsps_diff_freq_vq(int indexes[], float lsp[], int order)
1293
{
1294
int i,k,m;
1295
float lsp_hz[LPC_MAX];
1296
float lsp__hz[LPC_MAX];
1297
float dlsp[LPC_MAX];
1298
float dlsp_[LPC_MAX];
1299
float wt[LPC_MAX];
1300
const float * cb;
1301
float se;
1302
1303
for(i=0; i<LPC_ORD; i++) {
1304
wt[i] = 1.0;
1305
}
1306
1307
/* convert from radians to Hz so we can use human readable
1308
frequencies */
1309
1310
for(i=0; i<order; i++)
1311
lsp_hz[i] = (4000.0/PI)*lsp[i];
1312
1313
/* scalar quantisers for LSP differences 1..4 */
1314
1315
wt[0] = 1.0;
1316
for(i=0; i<4; i++) {
1317
if (i)
1318
dlsp[i] = lsp_hz[i] - lsp__hz[i-1];
1319
else
1320
dlsp[0] = lsp_hz[0];
1321
1322
k = lsp_cbd[i].k;
1323
m = lsp_cbd[i].m;
1324
cb = lsp_cbd[i].cb;
1325
indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se);
1326
dlsp_[i] = cb[indexes[i]*k];
1327
1328
if (i)
1329
lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
1330
else
1331
lsp__hz[0] = dlsp_[0];
1332
}
1333
1334
/* VQ LSPs 5,6,7,8,9,10 */
1335
1336
k = lsp_cbjnd[4].k;
1337
m = lsp_cbjnd[4].m;
1338
cb = lsp_cbjnd[4].cb;
1339
indexes[4] = quantise(cb, &lsp_hz[4], &wt[4], k, m, &se);
1340
}
1341
1342
1343
/*---------------------------------------------------------------------------*\
1344
1345
FUNCTION....: decode_lsps_diff_freq_vq()
1346
AUTHOR......: David Rowe
1347
DATE CREATED: 15 Nov 2011
1348
1349
From a vector of quantised LSP indexes, returns the quantised
1350
(floating point) LSPs.
1351
1352
\*---------------------------------------------------------------------------*/
1353
1354
void decode_lsps_diff_freq_vq(float lsp_[], int indexes[], int order)
1355
{
1356
int i,k,m;
1357
float dlsp_[LPC_MAX];
1358
float lsp__hz[LPC_MAX];
1359
const float * cb;
1360
1361
/* scalar LSP differences */
1362
1363
for(i=0; i<4; i++) {
1364
cb = lsp_cbd[i].cb;
1365
dlsp_[i] = cb[indexes[i]];
1366
if (i)
1367
lsp__hz[i] = lsp__hz[i-1] + dlsp_[i];
1368
else
1369
lsp__hz[0] = dlsp_[0];
1370
}
1371
1372
/* VQ */
1373
1374
k = lsp_cbjnd[4].k;
1375
m = lsp_cbjnd[4].m;
1376
cb = lsp_cbjnd[4].cb;
1377
for(i=4; i<order; i++)
1378
lsp__hz[i] = cb[indexes[4]*k+i-4];
1379
1380
/* convert back to radians */
1381
1382
for(i=0; i<order; i++)
1383
lsp_[i] = (PI/4000.0)*lsp__hz[i];
1384
}
1385
1386
1387
/*---------------------------------------------------------------------------*\
1388
1389
FUNCTION....: encode_lsps_diff_time()
1390
AUTHOR......: David Rowe
1391
DATE CREATED: 12 Sep 2012
1392
1393
Encode difference from preious frames's LSPs using
1394
3,3,2,2,2,2,1,1,1,1 scalar quantisers (18 bits total).
1395
1396
\*---------------------------------------------------------------------------*/
1397
1398
void encode_lsps_diff_time(int indexes[],
1399
float lsps[],
1400
float lsps__prev[],
1401
int order)
1402
{
1403
int i,k,m;
1404
float lsps_dt[LPC_ORD];
1405
float wt[LPC_MAX];
1406
const float * cb;
1407
float se;
1408
1409
/* Determine difference in time and convert from radians to Hz so
1410
we can use human readable frequencies */
1411
1412
for(i=0; i<LPC_ORD; i++) {
1413
lsps_dt[i] = (4000/PI)*(lsps[i] - lsps__prev[i]);
1414
}
1415
1416
/* scalar quantisers */
1417
1418
wt[0] = 1.0;
1419
for(i=0; i<order; i++) {
1420
k = lsp_cbdt[i].k;
1421
m = lsp_cbdt[i].m;
1422
cb = lsp_cbdt[i].cb;
1423
indexes[i] = quantise(cb, &lsps_dt[i], wt, k, m, &se);
1424
}
1425
1426
}
1427
1428
1429
/*---------------------------------------------------------------------------*\
1430
1431
FUNCTION....: decode_lsps_diff_time()
1432
AUTHOR......: David Rowe
1433
DATE CREATED: 15 Nov 2011
1434
1435
From a quantised LSP indexes, returns the quantised
1436
(floating point) LSPs.
1437
1438
\*---------------------------------------------------------------------------*/
1439
1440
void decode_lsps_diff_time(
1441
float lsps_[],
1442
int indexes[],
1443
float lsps__prev[],
1444
int order)
1445
{
1446
int i,k,m;
1447
const float * cb;
1448
1449
for(i=0; i<order; i++)
1450
lsps_[i] = lsps__prev[i];
1451
1452
for(i=0; i<order; i++) {
1453
k = lsp_cbdt[i].k;
1454
cb = lsp_cbdt[i].cb;
1455
lsps_[i] += (PI/4000.0)*cb[indexes[i]*k];
1456
}
1457
1458
}
1459
#endif
1460
1461
/*---------------------------------------------------------------------------*\
1462
1463
FUNCTION....: encode_lsps_vq()
1464
AUTHOR......: David Rowe
1465
DATE CREATED: 15 Feb 2012
1466
1467
Multi-stage VQ LSP quantiser developed by Jean-Marc Valin.
1468
1469
\*---------------------------------------------------------------------------*/
1470
1471
void encode_lsps_vq(int *indexes, float *x, float *xq, int ndim)
1472
{
1473
int i, n1, n2, n3;
1474
float err[LPC_ORD], err2[LPC_ORD], err3[LPC_ORD];
1475
float w[LPC_ORD], w2[LPC_ORD], w3[LPC_ORD];
1476
const float *codebook1 = lsp_cbjvm[0].cb;
1477
const float *codebook2 = lsp_cbjvm[1].cb;
1478
const float *codebook3 = lsp_cbjvm[2].cb;
1479
1480
assert(ndim <= LPC_ORD);
1481
1482
w[0] = MIN(x[0], x[1]-x[0]);
1483
for (i=1;i<ndim-1;i++)
1484
w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
1485
w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], PI-x[ndim-1]);
1486
1487
compute_weights(x, w, ndim);
1488
1489
n1 = find_nearest(codebook1, lsp_cbjvm[0].m, x, ndim);
1490
1491
for (i=0;i<ndim;i++)
1492
{
1493
xq[i] = codebook1[ndim*n1+i];
1494
err[i] = x[i] - xq[i];
1495
}
1496
for (i=0;i<ndim/2;i++)
1497
{
1498
err2[i] = err[2*i];
1499
err3[i] = err[2*i+1];
1500
w2[i] = w[2*i];
1501
w3[i] = w[2*i+1];
1502
}
1503
n2 = find_nearest_weighted(codebook2, lsp_cbjvm[1].m, err2, w2, ndim/2);
1504
n3 = find_nearest_weighted(codebook3, lsp_cbjvm[2].m, err3, w3, ndim/2);
1505
1506
indexes[0] = n1;
1507
indexes[1] = n2;
1508
indexes[2] = n3;
1509
}
1510
1511
1512
/*---------------------------------------------------------------------------*\
1513
1514
FUNCTION....: decode_lsps_vq()
1515
AUTHOR......: David Rowe
1516
DATE CREATED: 15 Feb 2012
1517
1518
\*---------------------------------------------------------------------------*/
1519
1520
void decode_lsps_vq(int *indexes, float *xq, int ndim)
1521
{
1522
int i, n1, n2, n3;
1523
const float *codebook1 = lsp_cbjvm[0].cb;
1524
const float *codebook2 = lsp_cbjvm[1].cb;
1525
const float *codebook3 = lsp_cbjvm[2].cb;
1526
1527
n1 = indexes[0];
1528
n2 = indexes[1];
1529
n3 = indexes[2];
1530
1531
for (i=0;i<ndim;i++)
1532
{
1533
xq[i] = codebook1[ndim*n1+i];
1534
}
1535
for (i=0;i<ndim/2;i++)
1536
{
1537
xq[2*i] += codebook2[ndim*n2/2+i];
1538
xq[2*i+1] += codebook3[ndim*n3/2+i];
1539
}
1540
}
1541
1542
1543
687
/*---------------------------------------------------------------------------*\
1544
688
1545
689
FUNCTION....: bw_expand_lsps()
1548
692
1549
693
Applies Bandwidth Expansion (BW) to a vector of LSPs. Prevents any
1550
694
two LSPs getting too close together after quantisation. We know
1551
from experiment that LSP quantisation errors < 12.5Hz (25Hz step
695
from experiment that LSP quantisation errors < 12.5Hz (25Hz setp
1552
696
size) are inaudible so we use that as the minimum LSP separation.
1553
697
1554
698
\*---------------------------------------------------------------------------*/
1559
703
{
1560
704
int i;
1561
705
1562
for(i=1; i<4; i++) {
1563
1564
if ((lsp[i] - lsp[i-1]) < 50.0*(PI/4000.0))
1565
lsp[i] = lsp[i-1] + 50.0*(PI/4000.0);
1566
1567
}
1568
1569
/* As quantiser gaps increased, larger BW expansion was required
1570
to prevent twinkly noises. This may need more experiment for
1571
different quanstisers.
1572
*/
1573
1574
for(i=4; i<order; i++) {
1575
if (lsp[i] - lsp[i-1] < 100.0*(PI/4000.0))
1576
lsp[i] = lsp[i-1] + 100.0*(PI/4000.0);
1577
}
1578
}
1579
1580
void bw_expand_lsps2(float lsp[],
1581
int order
1582
)
1583
{
1584
int i;
1585
1586
for(i=1; i<4; i++) {
1587
1588
if ((lsp[i] - lsp[i-1]) < 100.0*(PI/4000.0))
1589
lsp[i] = lsp[i-1] + 100.0*(PI/4000.0);
1590
1591
}
1592
1593
/* As quantiser gaps increased, larger BW expansion was required
1594
to prevent twinkly noises. This may need more experiment for
1595
different quanstisers.
1596
*/
1597
1598
for(i=4; i<order; i++) {
1599
if (lsp[i] - lsp[i-1] < 200.0*(PI/4000.0))
1600
lsp[i] = lsp[i-1] + 200.0*(PI/4000.0);
1601
}
1602
}
1603
1604
/*---------------------------------------------------------------------------*\
1605
1606
FUNCTION....: locate_lsps_jnd_steps()
1607
AUTHOR......: David Rowe
1608
DATE CREATED: 27/10/2011
1609
1610
Applies a form of Bandwidth Expansion (BW) to a vector of LSPs.
1611
Listening tests have determined that "quantising" the position of
1612
each LSP to the non-linear steps below introduces a "just noticable
1613
difference" in the synthesised speech.
1614
1615
This operation can be used before quantisation to limit the input
1616
data to the quantiser to a number of discrete steps.
1617
1618
This operation can also be used during quantisation as a form of
1619
hysteresis in the calculation of quantiser error. For example if
1620
the quantiser target of lsp1 is 500 Hz, candidate vectors with lsp1
1621
of 515 and 495 Hz sound effectively the same.
1622
1623
\*---------------------------------------------------------------------------*/
1624
1625
void locate_lsps_jnd_steps(float lsps[], int order)
1626
{
1627
int i;
1628
float lsp_hz, step;
1629
1630
assert(order == 10);
1631
1632
/* quantise to 25Hz steps */
1633
1634
step = 25;
1635
for(i=0; i<2; i++) {
1636
lsp_hz = lsps[i]*4000.0/PI;
1637
lsp_hz = floorf(lsp_hz/step + 0.5)*step;
1638
lsps[i] = lsp_hz*PI/4000.0;
1639
if (i) {
1640
if (lsps[i] == lsps[i-1])
1641
lsps[i] += step*PI/4000.0;
1642
1643
}
1644
}
1645
1646
/* quantise to 50Hz steps */
1647
1648
step = 50;
1649
for(i=2; i<4; i++) {
1650
lsp_hz = lsps[i]*4000.0/PI;
1651
lsp_hz = floorf(lsp_hz/step + 0.5)*step;
1652
lsps[i] = lsp_hz*PI/4000.0;
1653
if (i) {
1654
if (lsps[i] == lsps[i-1])
1655
lsps[i] += step*PI/4000.0;
1656
1657
}
1658
}
1659
1660
/* quantise to 100Hz steps */
1661
1662
step = 100;
1663
for(i=4; i<10; i++) {
1664
lsp_hz = lsps[i]*4000.0/PI;
1665
lsp_hz = floorf(lsp_hz/step + 0.5)*step;
1666
lsps[i] = lsp_hz*PI/4000.0;
1667
if (i) {
1668
if (lsps[i] == lsps[i-1])
1669
lsps[i] += step*PI/4000.0;
1670
1671
}
1672
}
1673
}
1674
706
for(i=1; i<5; i++) {
707
if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0))
708
lsp[i] = lsp[i-1] + PI*(12.5/4000.0);
709
}
710
711
/* As quantiser gaps increased, larger BW expansion was required
712
to prevent twinkly noises. This may need more experiment for
713
different quanstisers.
714
*/
715
716
for(i=5; i<8; i++) {
717
if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0))
718
lsp[i] = lsp[i-1] + PI*(25.0/4000.0);
719
}
720
for(i=8; i<order; i++) {
721
if (lsp[i] - lsp[i-1] < PI*(75.0/4000.0))
722
lsp[i] = lsp[i-1] + PI*(75.0/4000.0);
723
}
724
}
1675
725
1676
726
/*---------------------------------------------------------------------------*\
1677
727
1708
758
float e_max = E_MAX_DB;
1709
759
float norm;
1710
760
1711
e = 10.0*log10f(e);
761
e = 10.0*log10(e);
1712
762
norm = (e - e_min)/(e_max - e_min);
1713
index = floorf(E_LEVELS * norm + 0.5);
763
index = floor(E_LEVELS * norm + 0.5);
1714
764
if (index < 0 ) index = 0;
1715
765
if (index > (E_LEVELS-1)) index = E_LEVELS-1;
1716
766
1723
773
AUTHOR......: David Rowe
1724
774
DATE CREATED: 22/8/2010
1725
775
1726
Decodes energy using a E_LEVELS quantiser.
776
Decodes energy using a WO_BITS quantiser.
1727
777
1728
778
\*---------------------------------------------------------------------------*/
1729
779
1736
786
1737
787
step = (e_max - e_min)/E_LEVELS;
1738
788
e = e_min + step*(index);
1739
e = powf(10.0,e/10.0);
789
e = pow(10.0,e/10.0);
1740
790
1741
791
return e;
1742
792
}
1743
793
1744
#ifdef NOT_USED
794
/*---------------------------------------------------------------------------*\
795
796
FUNCTION....: encode_amplitudes()
797
AUTHOR......: David Rowe
798
DATE CREATED: 22/8/2010
799
800
Time domain LPC is used model the amplitudes which are then
801
converted to LSPs and quantised. So we don't actually encode the
802
amplitudes directly, rather we derive an equivalent representation
803
from the time domain speech.
804
805
\*---------------------------------------------------------------------------*/
806
807
void encode_amplitudes(int lsp_indexes[],
808
int *energy_index,
809
MODEL *model,
810
float Sn[],
811
float w[])
812
{
813
float lsps[LPC_ORD];
814
float ak[LPC_ORD+1];
815
float e;
816
817
e = speech_to_uq_lsps(lsps, ak, Sn, w, LPC_ORD);
818
encode_lsps(lsp_indexes, lsps, LPC_ORD);
819
*energy_index = encode_energy(e);
820
}
821
1745
822
/*---------------------------------------------------------------------------*\
1746
823
1747
824
FUNCTION....: decode_amplitudes()
1753
830
1754
831
\*---------------------------------------------------------------------------*/
1755
832
1756
float decode_amplitudes(kiss_fft_cfg fft_fwd_cfg,
1757
MODEL *model,
833
float decode_amplitudes(MODEL *model,
1758
834
float ak[],
1759
835
int lsp_indexes[],
1760
836
int energy_index,
1764
840
{
1765
841
float snr;
1766
842
1767
decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD);
843
decode_lsps(lsps, lsp_indexes, LPC_ORD);
1768
844
bw_expand_lsps(lsps, LPC_ORD);
1769
845
lsp_to_lpc(lsps, ak, LPC_ORD);
1770
846
*e = decode_energy(energy_index);
1771
aks_to_M2(ak, LPC_ORD, model, *e, &snr, 1, 0, 0, 1);
847
aks_to_M2(ak, LPC_ORD, model, *e, &snr, 1);
1772
848
apply_lpc_correction(model);
1773
849
1774
850
return snr;
1775
851
}
1776
#endif
1777
1778
static float ge_coeff[2] = {0.8, 0.9};
1779
1780
void compute_weights2(const float *x, const float *xp, float *w, int ndim)
1781
{
1782
w[0] = 30;
1783
w[1] = 1;
1784
if (x[1]<0)
1785
{
1786
w[0] *= .6;
1787
w[1] *= .3;
1788
}
1789
if (x[1]<-10)
1790
{
1791
w[0] *= .3;
1792
w[1] *= .3;
1793
}
1794
/* Higher weight if pitch is stable */
1795
if (fabsf(x[0]-xp[0])<.2)
1796
{
1797
w[0] *= 2;
1798
w[1] *= 1.5;
1799
} else if (fabsf(x[0]-xp[0])>.5) /* Lower if not stable */
1800
{
1801
w[0] *= .5;
1802
}
1803
1804
/* Lower weight for low energy */
1805
if (x[1] < xp[1]-10)
1806
{
1807
w[1] *= .5;
1808
}
1809
if (x[1] < xp[1]-20)
1810
{
1811
w[1] *= .5;
1812
}
1813
1814
//w[0] = 30;
1815
//w[1] = 1;
1816
1817
/* Square the weights because it's applied on the squared error */
1818
w[0] *= w[0];
1819
w[1] *= w[1];
1820
1821
}
1822
1823
/*---------------------------------------------------------------------------*\
1824
1825
FUNCTION....: quantise_WoE()
1826
AUTHOR......: Jean-Marc Valin & David Rowe
1827
DATE CREATED: 29 Feb 2012
1828
1829
Experimental joint Wo and LPC energy vector quantiser developed by
1830
Jean-Marc Valin. Exploits correlations between the difference in
1831
the log pitch and log energy from frame to frame. For example
1832
both the pitch and energy tend to only change by small amounts
1833
during voiced speech, however it is important that these changes be
1834
coded carefully. During unvoiced speech they both change a lot but
1835
the ear is less sensitve to errors so coarser quantisation is OK.
1836
1837
The ear is sensitive to log energy and loq pitch so we quantise in
1838
these domains. That way the error measure used to quantise the
1839
values is close to way the ear senses errors.
1840
1841
See http://jmspeex.livejournal.com/10446.html
1842
1843
\*---------------------------------------------------------------------------*/
1844
1845
void quantise_WoE(MODEL *model, float *e, float xq[])
1846
{
1847
int i, n1;
1848
float x[2];
1849
float err[2];
1850
float w[2];
1851
const float *codebook1 = ge_cb[0].cb;
1852
int nb_entries = ge_cb[0].m;
1853
int ndim = ge_cb[0].k;
1854
float Wo_min = TWO_PI/P_MAX;
1855
float Wo_max = TWO_PI/P_MIN;
1856
1857
x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2);
1858
x[1] = 10.0*log10f(1e-4 + *e);
1859
1860
compute_weights2(x, xq, w, ndim);
1861
for (i=0;i<ndim;i++)
1862
err[i] = x[i]-ge_coeff[i]*xq[i];
1863
n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim);
1864
1865
for (i=0;i<ndim;i++)
1866
{
1867
xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i];
1868
err[i] -= codebook1[ndim*n1+i];
1869
}
1870
1871
/*
1872
x = log2(4000*Wo/(PI*50));
1873
2^x = 4000*Wo/(PI*50)
1874
Wo = (2^x)*(PI*50)/4000;
1875
*/
1876
1877
model->Wo = powf(2.0, xq[0])*(PI*50.0)/4000.0;
1878
1879
/* bit errors can make us go out of range leading to all sorts of
1880
probs like seg faults */
1881
1882
if (model->Wo > Wo_max) model->Wo = Wo_max;
1883
if (model->Wo < Wo_min) model->Wo = Wo_min;
1884
1885
model->L = PI/model->Wo; /* if we quantise Wo re-compute L */
1886
1887
*e = powf(10.0, xq[1]/10.0);
1888
}
1889
1890
/*---------------------------------------------------------------------------*\
1891
1892
FUNCTION....: encode_WoE()
1893
AUTHOR......: Jean-Marc Valin & David Rowe
1894
DATE CREATED: 11 May 2012
1895
1896
Joint Wo and LPC energy vector quantiser developed my Jean-Marc
1897
Valin. Returns index, and updated states xq[].
1898
1899
\*---------------------------------------------------------------------------*/
1900
1901
int encode_WoE(MODEL *model, float e, float xq[])
1902
{
1903
int i, n1;
1904
float x[2];
1905
float err[2];
1906
float w[2];
1907
const float *codebook1 = ge_cb[0].cb;
1908
int nb_entries = ge_cb[0].m;
1909
int ndim = ge_cb[0].k;
1910
1911
assert((1<<WO_E_BITS) == nb_entries);
1912
1913
if (e < 0.0) e = 0; /* occasional small negative energies due LPC round off I guess */
1914
1915
x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2);
1916
x[1] = 10.0*log10f(1e-4 + e);
1917
1918
compute_weights2(x, xq, w, ndim);
1919
for (i=0;i<ndim;i++)
1920
err[i] = x[i]-ge_coeff[i]*xq[i];
1921
n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim);
1922
1923
for (i=0;i<ndim;i++)
1924
{
1925
xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i];
1926
err[i] -= codebook1[ndim*n1+i];
1927
}
1928
1929
//printf("enc: %f %f (%f)(%f) \n", xq[0], xq[1], e, 10.0*log10(1e-4 + e));
1930
return n1;
1931
}
1932
1933
1934
/*---------------------------------------------------------------------------*\
1935
1936
FUNCTION....: decode_WoE()
1937
AUTHOR......: Jean-Marc Valin & David Rowe
1938
DATE CREATED: 11 May 2012
1939
1940
Joint Wo and LPC energy vector quantiser developed my Jean-Marc
1941
Valin. Given index and states xq[], returns Wo & E, and updates
1942
states xq[].
1943
1944
\*---------------------------------------------------------------------------*/
1945
1946
void decode_WoE(MODEL *model, float *e, float xq[], int n1)
1947
{
1948
int i;
1949
float err[2];
1950
const float *codebook1 = ge_cb[0].cb;
1951
int ndim = ge_cb[0].k;
1952
float Wo_min = TWO_PI/P_MAX;
1953
float Wo_max = TWO_PI/P_MIN;
1954
1955
for (i=0;i<ndim;i++)
1956
{
1957
xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i];
1958
err[i] -= codebook1[ndim*n1+i];
1959
}
1960
1961
//printf("dec: %f %f\n", xq[0], xq[1]);
1962
model->Wo = powf(2.0, xq[0])*(PI*50.0)/4000.0;
1963
1964
/* bit errors can make us go out of range leading to all sorts of
1965
probs like seg faults */
1966
1967
if (model->Wo > Wo_max) model->Wo = Wo_max;
1968
if (model->Wo < Wo_min) model->Wo = Wo_min;
1969
1970
model->L = PI/model->Wo; /* if we quantise Wo re-compute L */
1971
1972
*e = powf(10.0, xq[1]/10.0);
1973
}
1974
Loggerhead is a web-based interface for Breezy
Version: 2.0.1