diff options
author | erdgeist@erdgeist.org <erdgeist@bauklotz.fritz.box> | 2019-07-04 23:26:09 +0200 |
---|---|---|
committer | erdgeist@erdgeist.org <erdgeist@bauklotz.fritz.box> | 2019-07-04 23:26:09 +0200 |
commit | f02dfce6e6c34b3d8a7b8a0e784b506178e331fa (patch) | |
tree | 45556e6104242d4702689760433d7321ae74ec17 /quantise.c |
stripdown of version 0.9
Diffstat (limited to 'quantise.c')
-rw-r--r-- | quantise.c | 2051 |
1 files changed, 2051 insertions, 0 deletions
diff --git a/quantise.c b/quantise.c new file mode 100644 index 0000000..37bf8be --- /dev/null +++ b/quantise.c | |||
@@ -0,0 +1,2051 @@ | |||
1 | /*---------------------------------------------------------------------------*\ | ||
2 | |||
3 | FILE........: quantise.c | ||
4 | AUTHOR......: David Rowe | ||
5 | DATE CREATED: 31/5/92 | ||
6 | |||
7 | Quantisation functions for the sinusoidal coder. | ||
8 | |||
9 | \*---------------------------------------------------------------------------*/ | ||
10 | |||
11 | /* | ||
12 | All rights reserved. | ||
13 | |||
14 | This program is free software; you can redistribute it and/or modify | ||
15 | it under the terms of the GNU Lesser General Public License version 2.1, as | ||
16 | published by the Free Software Foundation. This program is | ||
17 | distributed in the hope that it will be useful, but WITHOUT ANY | ||
18 | WARRANTY; without even the implied warranty of MERCHANTABILITY or | ||
19 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public | ||
20 | License for more details. | ||
21 | |||
22 | You should have received a copy of the GNU Lesser General Public License | ||
23 | along with this program; if not, see <http://www.gnu.org/licenses/>. | ||
24 | |||
25 | */ | ||
26 | |||
27 | #include <assert.h> | ||
28 | #include <ctype.h> | ||
29 | #include <stdio.h> | ||
30 | #include <stdlib.h> | ||
31 | #include <string.h> | ||
32 | #include <math.h> | ||
33 | |||
34 | #include "defines.h" | ||
35 | #include "dump.h" | ||
36 | #include "quantise.h" | ||
37 | #include "lpc.h" | ||
38 | #include "lsp.h" | ||
39 | #include "codec2_fft.h" | ||
40 | #include "phase.h" | ||
41 | #include "mbest.h" | ||
42 | |||
43 | #undef PROFILE | ||
44 | #include "machdep.h" | ||
45 | |||
46 | #define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */ | ||
47 | |||
48 | /*---------------------------------------------------------------------------*\ | ||
49 | |||
50 | FUNCTION HEADERS | ||
51 | |||
52 | \*---------------------------------------------------------------------------*/ | ||
53 | |||
54 | float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[], | ||
55 | int m_pitch, int order); | ||
56 | |||
57 | /*---------------------------------------------------------------------------*\ | ||
58 | |||
59 | FUNCTIONS | ||
60 | |||
61 | \*---------------------------------------------------------------------------*/ | ||
62 | |||
63 | int lsp_bits(int i) { | ||
64 | return lsp_cb[i].log2m; | ||
65 | } | ||
66 | |||
67 | int lspd_bits(int i) { | ||
68 | return lsp_cbd[i].log2m; | ||
69 | } | ||
70 | |||
71 | #ifndef CORTEX_M4 | ||
72 | int mel_bits(int i) { | ||
73 | return mel_cb[i].log2m; | ||
74 | } | ||
75 | |||
76 | int lspmelvq_cb_bits(int i) { | ||
77 | return lspmelvq_cb[i].log2m; | ||
78 | } | ||
79 | #endif | ||
80 | |||
81 | #ifdef __EXPERIMENTAL__ | ||
82 | int lspdt_bits(int i) { | ||
83 | return lsp_cbdt[i].log2m; | ||
84 | } | ||
85 | #endif | ||
86 | |||
87 | int 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 | |||
100 | void quantise_init() | ||
101 | { | ||
102 | } | ||
103 | |||
104 | /*---------------------------------------------------------------------------*\ | ||
105 | |||
106 | quantise | ||
107 | |||
108 | Quantises vec by choosing the nearest vector in codebook cb, and | ||
109 | returns the vector index. The squared error of the quantised vector | ||
110 | is added to se. | ||
111 | |||
112 | \*---------------------------------------------------------------------------*/ | ||
113 | |||
114 | long quantise(const float * cb, float vec[], float w[], int k, int m, float *se) | ||
115 | /* float cb[][K]; current VQ codebook */ | ||
116 | /* float vec[]; vector to quantise */ | ||
117 | /* float w[]; weighting vector */ | ||
118 | /* int k; dimension of vectors */ | ||
119 | /* int m; size of codebook */ | ||
120 | /* float *se; accumulated squared error */ | ||
121 | { | ||
122 | float e; /* current error */ | ||
123 | long besti; /* best index so far */ | ||
124 | float beste; /* best error so far */ | ||
125 | long j; | ||
126 | int i; | ||
127 | 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 | |||
158 | void encode_lspds_scalar( | ||
159 | int indexes[], | ||
160 | float lsp[], | ||
161 | int order | ||
162 | ) | ||
163 | { | ||
164 | int i,k,m; | ||
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 | } | ||
176 | |||
177 | /* convert from radians to Hz so we can use human readable | ||
178 | frequencies */ | ||
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 | } | ||
209 | |||
210 | } | ||
211 | |||
212 | |||
213 | void 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 | |||
240 | } | ||
241 | |||
242 | #ifdef __EXPERIMENTAL__ | ||
243 | /*---------------------------------------------------------------------------*\ | ||
244 | |||
245 | lspvq_quantise | ||
246 | |||
247 | Vector LSP quantiser. | ||
248 | |||
249 | \*---------------------------------------------------------------------------*/ | ||
250 | |||
251 | void lspvq_quantise( | ||
252 | float lsp[], | ||
253 | float lsp_[], | ||
254 | int order | ||
255 | ) | ||
256 | { | ||
257 | int i,k,m,ncb, nlsp; | ||
258 | float wt[order], lsp_hz[order]; | ||
259 | const float *cb; | ||
260 | float se; | ||
261 | int index; | ||
262 | |||
263 | for(i=0; i<order; i++) { | ||
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 | |||
307 | Experimental JND LSP quantiser. | ||
308 | |||
309 | \*---------------------------------------------------------------------------*/ | ||
310 | |||
311 | void lspjnd_quantise(float lsps[], float lsps_[], int order) | ||
312 | { | ||
313 | int i,k,m; | ||
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 | |||
325 | for(i=0; i<order; i++) { | ||
326 | lsps_hz[i] = lsps[i]*(4000.0/PI); | ||
327 | lsps_[i] = lsps[i]; | ||
328 | } | ||
329 | |||
330 | /* simple uniform scalar quantisers */ | ||
331 | |||
332 | for(i=0; i<4; i++) { | ||
333 | k = lsp_cbjnd[i].k; | ||
334 | m = lsp_cbjnd[i].m; | ||
335 | cb = lsp_cbjnd[i].cb; | ||
336 | index = quantise(cb, &lsps_hz[i], wt, k, m, &se); | ||
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 | } | ||
354 | |||
355 | void compute_weights(const float *x, float *w, int ndim); | ||
356 | |||
357 | /*---------------------------------------------------------------------------*\ | ||
358 | |||
359 | lspdt_quantise | ||
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 | |||
375 | void 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 | |||
394 | for(i=0; i<LPC_ORD; i++) { | ||
395 | lsps_dt[i] = lsps[i] - lsps__prev[i]; | ||
396 | lsps_[i] = lsps__prev[i]; | ||
397 | } | ||
398 | |||
399 | //#define TRY_LSPDT_VQ | ||
400 | #ifdef TRY_LSPDT_VQ | ||
401 | /* this actually improves speech a bit, but 40ms updates works surprsingly well.... */ | ||
402 | k = lsp_cbdt[0].k; | ||
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 | |||
411 | } | ||
412 | #endif | ||
413 | |||
414 | #define MIN(a,b) ((a)<(b)?(a):(b)) | ||
415 | #define MAX_ENTRIES 16384 | ||
416 | |||
417 | void compute_weights(const float *x, float *w, int ndim) | ||
418 | { | ||
419 | int i; | ||
420 | w[0] = MIN(x[0], x[1]-x[0]); | ||
421 | for (i=1;i<ndim-1;i++) | ||
422 | w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]); | ||
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 | |||
435 | void compute_weights_anssi_mode2(const float *x, float *w, int ndim) | ||
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 | } | ||
464 | #endif | ||
465 | |||
466 | int find_nearest(const float *codebook, int nb_entries, float *x, int ndim) | ||
467 | { | ||
468 | int i, j; | ||
469 | float min_dist = 1e15; | ||
470 | int nearest = 0; | ||
471 | |||
472 | for (i=0;i<nb_entries;i++) | ||
473 | { | ||
474 | float dist=0; | ||
475 | for (j=0;j<ndim;j++) | ||
476 | dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]); | ||
477 | if (dist<min_dist) | ||
478 | { | ||
479 | min_dist = dist; | ||
480 | nearest = i; | ||
481 | } | ||
482 | } | ||
483 | return nearest; | ||
484 | } | ||
485 | |||
486 | int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim) | ||
487 | { | ||
488 | int i, j; | ||
489 | float min_dist = 1e15; | ||
490 | int nearest = 0; | ||
491 | |||
492 | for (i=0;i<nb_entries;i++) | ||
493 | { | ||
494 | float dist=0; | ||
495 | for (j=0;j<ndim;j++) | ||
496 | dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]); | ||
497 | if (dist<min_dist) | ||
498 | { | ||
499 | min_dist = dist; | ||
500 | nearest = i; | ||
501 | } | ||
502 | } | ||
503 | return nearest; | ||
504 | } | ||
505 | |||
506 | void lspjvm_quantise(float *x, float *xq, int order) | ||
507 | { | ||
508 | int i, n1, n2, n3; | ||
509 | float err[order], err2[order], err3[order]; | ||
510 | float w[order], w2[order], w3[order]; | ||
511 | const float *codebook1 = lsp_cbjvm[0].cb; | ||
512 | const float *codebook2 = lsp_cbjvm[1].cb; | ||
513 | const float *codebook3 = lsp_cbjvm[2].cb; | ||
514 | |||
515 | w[0] = MIN(x[0], x[1]-x[0]); | ||
516 | for (i=1;i<order-1;i++) | ||
517 | w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]); | ||
518 | w[order-1] = MIN(x[order-1]-x[order-2], PI-x[order-1]); | ||
519 | |||
520 | compute_weights(x, w, order); | ||
521 | |||
522 | n1 = find_nearest(codebook1, lsp_cbjvm[0].m, x, order); | ||
523 | |||
524 | for (i=0;i<order;i++) | ||
525 | { | ||
526 | xq[i] = codebook1[order*n1+i]; | ||
527 | err[i] = x[i] - xq[i]; | ||
528 | } | ||
529 | for (i=0;i<order/2;i++) | ||
530 | { | ||
531 | err2[i] = err[2*i]; | ||
532 | err3[i] = err[2*i+1]; | ||
533 | w2[i] = w[2*i]; | ||
534 | w3[i] = w[2*i+1]; | ||
535 | } | ||
536 | n2 = find_nearest_weighted(codebook2, lsp_cbjvm[1].m, err2, w2, order/2); | ||
537 | n3 = find_nearest_weighted(codebook3, lsp_cbjvm[2].m, err3, w3, order/2); | ||
538 | |||
539 | for (i=0;i<order/2;i++) | ||
540 | { | ||
541 | xq[2*i] += codebook2[order*n2/2+i]; | ||
542 | xq[2*i+1] += codebook3[order*n3/2+i]; | ||
543 | } | ||
544 | } | ||
545 | |||
546 | |||
547 | #ifndef CORTEX_M4 | ||
548 | /* simple (non mbest) 6th order LSP MEL VQ quantiser. Returns MSE of result */ | ||
549 | |||
550 | float 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 | |||
594 | float 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 | |||
661 | void 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 | |||
676 | int 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 | } | ||
694 | |||
695 | void force_min_lsp_dist(float lsp[], int order) | ||
696 | { | ||
697 | int i; | ||
698 | |||
699 | for(i=1; i<order; i++) | ||
700 | if ((lsp[i]-lsp[i-1]) < 0.01) { | ||
701 | lsp[i] += 0.01; | ||
702 | } | ||
703 | } | ||
704 | |||
705 | |||
706 | /*---------------------------------------------------------------------------*\ | ||
707 | |||
708 | lpc_post_filter() | ||
709 | |||
710 | Applies a post filter to the LPC synthesis filter power spectrum | ||
711 | Pw, which supresses the inter-formant energy. | ||
712 | |||
713 | 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 | ||
715 | of this text is on the MBE vocoder, and this is a freq domain | ||
716 | adaptation of post filtering commonly used in CELP. | ||
717 | |||
718 | I used the Octave simulation lpcpf.m to get an understanding of the | ||
719 | algorithm. | ||
720 | |||
721 | Requires two more FFTs which is significantly more MIPs. However | ||
722 | it should be possible to implement this more efficiently in the | ||
723 | time domain. Just not sure how to handle relative time delays | ||
724 | between the synthesis stage and updating these coeffs. A smaller | ||
725 | FFT size might also be accetable to save CPU. | ||
726 | |||
727 | TODO: | ||
728 | [ ] sync var names between Octave and C version | ||
729 | [ ] doc gain normalisation | ||
730 | [ ] I think the first FFT is not rqd as we do the same | ||
731 | thing in aks_to_M2(). | ||
732 | |||
733 | \*---------------------------------------------------------------------------*/ | ||
734 | |||
735 | void 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) | ||
737 | { | ||
738 | int i; | ||
739 | float x[FFT_ENC]; /* input to FFTs */ | ||
740 | COMP Ww[FFT_ENC/2+1]; /* weighting spectrum */ | ||
741 | float Rw[FFT_ENC/2+1]; /* R = WA */ | ||
742 | float e_before, e_after, gain; | ||
743 | float Pfw; | ||
744 | float max_Rw, min_Rw; | ||
745 | float coeff; | ||
746 | 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 | |||
756 | x[0] = ak[0]; | ||
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 | |||
764 | PROFILE_SAMPLE_AND_LOG(tfft2, taw, " fft2"); | ||
765 | |||
766 | for(i=0; i<FFT_ENC/2; i++) { | ||
767 | Ww[i].real = Ww[i].real*Ww[i].real + Ww[i].imag*Ww[i].imag; | ||
768 | } | ||
769 | |||
770 | PROFILE_SAMPLE_AND_LOG(tww, tfft2, " Ww"); | ||
771 | |||
772 | /* Determined combined filter R = WA ---------------------------*/ | ||
773 | |||
774 | max_Rw = 0.0; min_Rw = 1E32; | ||
775 | for(i=0; i<FFT_ENC/2; i++) { | ||
776 | Rw[i] = sqrtf(Ww[i].real * Pw[i]); | ||
777 | if (Rw[i] > max_Rw) | ||
778 | max_Rw = Rw[i]; | ||
779 | if (Rw[i] < min_Rw) | ||
780 | min_Rw = Rw[i]; | ||
781 | |||
782 | } | ||
783 | |||
784 | PROFILE_SAMPLE_AND_LOG(tr, tww, " R"); | ||
785 | |||
786 | #ifdef DUMP | ||
787 | if (dump) | ||
788 | dump_Rw(Rw); | ||
789 | #endif | ||
790 | |||
791 | /* create post filter mag spectrum and apply ------------------*/ | ||
792 | |||
793 | /* measure energy before post filtering */ | ||
794 | |||
795 | e_before = 1E-4; | ||
796 | for(i=0; i<FFT_ENC/2; i++) | ||
797 | e_before += Pw[i]; | ||
798 | |||
799 | /* apply post filter and measure energy */ | ||
800 | |||
801 | #ifdef DUMP | ||
802 | if (dump) | ||
803 | dump_Pwb(Pw); | ||
804 | #endif | ||
805 | |||
806 | |||
807 | e_after = 1E-4; | ||
808 | for(i=0; i<FFT_ENC/2; i++) { | ||
809 | Pfw = powf(Rw[i], beta); | ||
810 | Pw[i] *= Pfw * Pfw; | ||
811 | e_after += Pw[i]; | ||
812 | } | ||
813 | gain = e_before/e_after; | ||
814 | |||
815 | /* apply gain factor to normalise energy, and LPC Energy */ | ||
816 | |||
817 | gain *= E; | ||
818 | for(i=0; i<FFT_ENC/2; i++) { | ||
819 | Pw[i] *= gain; | ||
820 | } | ||
821 | |||
822 | if (bass_boost) { | ||
823 | /* add 3dB to first 1 kHz to account for LP effect of PF */ | ||
824 | |||
825 | for(i=0; i<FFT_ENC/8; i++) { | ||
826 | Pw[i] *= 1.4*1.4; | ||
827 | } | ||
828 | } | ||
829 | |||
830 | PROFILE_SAMPLE_AND_LOG2(tr, " filt"); | ||
831 | } | ||
832 | |||
833 | |||
834 | /*---------------------------------------------------------------------------*\ | ||
835 | |||
836 | aks_to_M2() | ||
837 | |||
838 | Transforms the linear prediction coefficients to spectral amplitude | ||
839 | samples. This function determines A(m) from the average energy per | ||
840 | band using an FFT. | ||
841 | |||
842 | \*---------------------------------------------------------------------------*/ | ||
843 | |||
844 | void aks_to_M2( | ||
845 | codec2_fftr_cfg fftr_fwd_cfg, | ||
846 | float ak[], /* LPC's */ | ||
847 | int order, | ||
848 | MODEL *model, /* sinusoidal model parameters for this frame */ | ||
849 | float E, /* energy term */ | ||
850 | float *snr, /* signal to noise ratio for this frame in dB */ | ||
851 | int dump, /* true to dump sample to dump file */ | ||
852 | int sim_pf, /* true to simulate a post filter */ | ||
853 | int pf, /* true to enable actual LPC post filter */ | ||
854 | int bass_boost, /* enable LPC filter 0-1kHz 3dB boost */ | ||
855 | float beta, | ||
856 | float gamma, /* LPC post filter parameters */ | ||
857 | COMP Aw[] /* output power spectrum */ | ||
858 | ) | ||
859 | { | ||
860 | int i,m; /* loop variables */ | ||
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; | ||
866 | PROFILE_VAR(tstart, tfft, tpw, tpf); | ||
867 | |||
868 | PROFILE_SAMPLE(tstart); | ||
869 | |||
870 | r = TWO_PI/(FFT_ENC); | ||
871 | |||
872 | /* Determine DFT of A(exp(jw)) --------------------------------------------*/ | ||
873 | { | ||
874 | float a[FFT_ENC]; /* input to FFT for power spectrum */ | ||
875 | |||
876 | for(i=0; i<FFT_ENC; i++) { | ||
877 | a[i] = 0.0; | ||
878 | } | ||
879 | |||
880 | for(i=0; i<=order; i++) | ||
881 | a[i] = ak[i]; | ||
882 | codec2_fftr(fftr_fwd_cfg, a, Aw); | ||
883 | } | ||
884 | PROFILE_SAMPLE_AND_LOG(tfft, tstart, " fft"); | ||
885 | |||
886 | /* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/ | ||
887 | |||
888 | float Pw[FFT_ENC/2]; | ||
889 | |||
890 | #ifndef FDV_ARM_MATH | ||
891 | 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); | ||
893 | } | ||
894 | #else | ||
895 | // this difference may seem strange, but the gcc for STM32F4 generates almost 5 times | ||
896 | // faster code with the two loops: 1120 ms -> 242 ms | ||
897 | // so please leave it as is or improve further | ||
898 | // since this code is called 4 times it results in almost 4ms gain (21ms -> 17ms per audio frame decode @ 1300 ) | ||
899 | |||
900 | for(i=0; i<FFT_ENC/2; i++) | ||
901 | { | ||
902 | Pw[i] = Aw[i].real * Aw[i].real + Aw[i].imag * Aw[i].imag + 1E-6; | ||
903 | } | ||
904 | for(i=0; i<FFT_ENC/2; i++) { | ||
905 | Pw[i] = 1.0/(Pw[i]); | ||
906 | } | ||
907 | #endif | ||
908 | |||
909 | PROFILE_SAMPLE_AND_LOG(tpw, tfft, " Pw"); | ||
910 | |||
911 | if (pf) | ||
912 | lpc_post_filter(fftr_fwd_cfg, Pw, ak, order, dump, beta, gamma, bass_boost, E); | ||
913 | else { | ||
914 | for(i=0; i<FFT_ENC/2; i++) { | ||
915 | Pw[i] *= E; | ||
916 | } | ||
917 | } | ||
918 | |||
919 | PROFILE_SAMPLE_AND_LOG(tpf, tpw, " LPC post filter"); | ||
920 | |||
921 | #ifdef DUMP | ||
922 | if (dump) | ||
923 | dump_Pw(Pw); | ||
924 | #endif | ||
925 | |||
926 | /* Determine magnitudes from P(w) ----------------------------------------*/ | ||
927 | |||
928 | /* when used just by decoder {A} might be all zeroes so init signal | ||
929 | and noise to prevent log(0) errors */ | ||
930 | |||
931 | signal = 1E-30; noise = 1E-32; | ||
932 | |||
933 | for(m=1; m<=model->L; m++) { | ||
934 | am = (int)((m - 0.5)*model->Wo/r + 0.5); | ||
935 | bm = (int)((m + 0.5)*model->Wo/r + 0.5); | ||
936 | |||
937 | // FIXME: With arm_rfft_fast_f32 we have to use this | ||
938 | // otherwise sometimes a to high bm is calculated | ||
939 | // which causes trouble later in the calculation | ||
940 | // chain | ||
941 | // it seems for some reason model->Wo is calculated somewhat too high | ||
942 | if (bm>FFT_ENC/2) | ||
943 | { | ||
944 | bm = FFT_ENC/2; | ||
945 | } | ||
946 | Em = 0.0; | ||
947 | |||
948 | for(i=am; i<bm; i++) | ||
949 | Em += Pw[i]; | ||
950 | Am = sqrtf(Em); | ||
951 | |||
952 | signal += model->A[m]*model->A[m]; | ||
953 | noise += (model->A[m] - Am)*(model->A[m] - Am); | ||
954 | |||
955 | /* This code significantly improves perf of LPC model, in | ||
956 | particular when combined with phase0. The LPC spectrum tends | ||
957 | to track just under the peaks of the spectral envelope, and | ||
958 | just above nulls. This algorithm does the reverse to | ||
959 | compensate - raising the amplitudes of spectral peaks, while | ||
960 | attenuating the null. This enhances the formants, and | ||
961 | supresses the energy between formants. */ | ||
962 | |||
963 | if (sim_pf) { | ||
964 | if (Am > model->A[m]) | ||
965 | Am *= 0.7; | ||
966 | if (Am < model->A[m]) | ||
967 | Am *= 1.4; | ||
968 | } | ||
969 | model->A[m] = Am; | ||
970 | } | ||
971 | *snr = 10.0*log10f(signal/noise); | ||
972 | |||
973 | PROFILE_SAMPLE_AND_LOG2(tpf, " rec"); | ||
974 | } | ||
975 | |||
976 | /*---------------------------------------------------------------------------*\ | ||
977 | |||
978 | FUNCTION....: encode_Wo() | ||
979 | AUTHOR......: David Rowe | ||
980 | DATE CREATED: 22/8/2010 | ||
981 | |||
982 | Encodes Wo using a WO_LEVELS quantiser. | ||
983 | |||
984 | \*---------------------------------------------------------------------------*/ | ||
985 | |||
986 | int encode_Wo(C2CONST *c2const, float Wo, int bits) | ||
987 | { | ||
988 | int index, Wo_levels = 1<<bits; | ||
989 | float Wo_min = c2const->Wo_min; | ||
990 | float Wo_max = c2const->Wo_max; | ||
991 | float norm; | ||
992 | |||
993 | norm = (Wo - Wo_min)/(Wo_max - Wo_min); | ||
994 | index = floorf(Wo_levels * norm + 0.5); | ||
995 | if (index < 0 ) index = 0; | ||
996 | if (index > (Wo_levels-1)) index = Wo_levels-1; | ||
997 | |||
998 | return index; | ||
999 | } | ||
1000 | |||
1001 | /*---------------------------------------------------------------------------*\ | ||
1002 | |||
1003 | FUNCTION....: decode_Wo() | ||
1004 | AUTHOR......: David Rowe | ||
1005 | DATE CREATED: 22/8/2010 | ||
1006 | |||
1007 | Decodes Wo using a WO_LEVELS quantiser. | ||
1008 | |||
1009 | \*---------------------------------------------------------------------------*/ | ||
1010 | |||
1011 | float decode_Wo(C2CONST *c2const, int index, int bits) | ||
1012 | { | ||
1013 | float Wo_min = c2const->Wo_min; | ||
1014 | float Wo_max = c2const->Wo_max; | ||
1015 | float step; | ||
1016 | float Wo; | ||
1017 | int Wo_levels = 1<<bits; | ||
1018 | |||
1019 | step = (Wo_max - Wo_min)/Wo_levels; | ||
1020 | Wo = Wo_min + step*(index); | ||
1021 | |||
1022 | return Wo; | ||
1023 | } | ||
1024 | |||
1025 | /*---------------------------------------------------------------------------*\ | ||
1026 | |||
1027 | FUNCTION....: encode_log_Wo() | ||
1028 | AUTHOR......: David Rowe | ||
1029 | DATE CREATED: 22/8/2010 | ||
1030 | |||
1031 | Encodes Wo in the log domain using a WO_LEVELS quantiser. | ||
1032 | |||
1033 | \*---------------------------------------------------------------------------*/ | ||
1034 | |||
1035 | int encode_log_Wo(C2CONST *c2const, float Wo, int bits) | ||
1036 | { | ||
1037 | int index, Wo_levels = 1<<bits; | ||
1038 | float Wo_min = c2const->Wo_min; | ||
1039 | float Wo_max = c2const->Wo_max; | ||
1040 | float norm; | ||
1041 | |||
1042 | norm = (log10f(Wo) - log10f(Wo_min))/(log10f(Wo_max) - log10f(Wo_min)); | ||
1043 | index = floorf(Wo_levels * norm + 0.5); | ||
1044 | if (index < 0 ) index = 0; | ||
1045 | if (index > (Wo_levels-1)) index = Wo_levels-1; | ||
1046 | |||
1047 | return index; | ||
1048 | } | ||
1049 | |||
1050 | /*---------------------------------------------------------------------------*\ | ||
1051 | |||
1052 | FUNCTION....: decode_log_Wo() | ||
1053 | AUTHOR......: David Rowe | ||
1054 | DATE CREATED: 22/8/2010 | ||
1055 | |||
1056 | Decodes Wo using a WO_LEVELS quantiser in the log domain. | ||
1057 | |||
1058 | \*---------------------------------------------------------------------------*/ | ||
1059 | |||
1060 | float decode_log_Wo(C2CONST *c2const, int index, int bits) | ||
1061 | { | ||
1062 | float Wo_min = c2const->Wo_min; | ||
1063 | float Wo_max = c2const->Wo_max; | ||
1064 | float step; | ||
1065 | float Wo; | ||
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 | |||
1085 | int 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 | |||
1124 | float 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 | |||
1147 | if (Wo > Wo_max) Wo = Wo_max; | ||
1148 | if (Wo < Wo_min) Wo = Wo_min; | ||
1149 | |||
1150 | return Wo; | ||
1151 | } | ||
1152 | #endif | ||
1153 | |||
1154 | /*---------------------------------------------------------------------------*\ | ||
1155 | |||
1156 | FUNCTION....: speech_to_uq_lsps() | ||
1157 | AUTHOR......: David Rowe | ||
1158 | DATE CREATED: 22/8/2010 | ||
1159 | |||
1160 | Analyse a windowed frame of time domain speech to determine LPCs | ||
1161 | which are the converted to LSPs for quantisation and transmission | ||
1162 | over the channel. | ||
1163 | |||
1164 | \*---------------------------------------------------------------------------*/ | ||
1165 | |||
1166 | float speech_to_uq_lsps(float lsp[], | ||
1167 | float ak[], | ||
1168 | float Sn[], | ||
1169 | float w[], | ||
1170 | int m_pitch, | ||
1171 | int order | ||
1172 | ) | ||
1173 | { | ||
1174 | int i, roots; | ||
1175 | float Wn[m_pitch]; | ||
1176 | float R[order+1]; | ||
1177 | float e, E; | ||
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 | |||
1185 | /* trap 0 energy case as LPC analysis will fail */ | ||
1186 | |||
1187 | if (e == 0.0) { | ||
1188 | for(i=0; i<order; i++) | ||
1189 | lsp[i] = (PI/order)*(float)i; | ||
1190 | return 0.0; | ||
1191 | } | ||
1192 | |||
1193 | autocorrelate(Wn, R, m_pitch, order); | ||
1194 | levinson_durbin(R, ak, order); | ||
1195 | |||
1196 | E = 0.0; | ||
1197 | for(i=0; i<=order; i++) | ||
1198 | E += ak[i]*R[i]; | ||
1199 | |||
1200 | /* 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 | ||
1202 | after energy calculation to avoid -ve energy values. | ||
1203 | */ | ||
1204 | |||
1205 | for(i=0; i<=order; i++) | ||
1206 | ak[i] *= powf(0.994,(float)i); | ||
1207 | |||
1208 | roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1); | ||
1209 | if (roots != order) { | ||
1210 | /* if root finding fails use some benign LSP values instead */ | ||
1211 | for(i=0; i<order; i++) | ||
1212 | lsp[i] = (PI/order)*(float)i; | ||
1213 | } | ||
1214 | |||
1215 | return E; | ||
1216 | } | ||
1217 | |||
1218 | /*---------------------------------------------------------------------------*\ | ||
1219 | |||
1220 | FUNCTION....: encode_lsps_scalar() | ||
1221 | AUTHOR......: David Rowe | ||
1222 | DATE CREATED: 22/8/2010 | ||
1223 | |||
1224 | Thirty-six bit sclar LSP quantiser. From a vector of unquantised | ||
1225 | (floating point) LSPs finds the quantised LSP indexes. | ||
1226 | |||
1227 | \*---------------------------------------------------------------------------*/ | ||
1228 | |||
1229 | void 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 | |||
1265 | void 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 | |||
1292 | \*---------------------------------------------------------------------------*/ | ||
1293 | |||
1294 | void encode_mels_scalar(int indexes[], float mels[], int order) | ||
1295 | { | ||
1296 | int i,m; | ||
1297 | float wt[1]; | ||
1298 | const float * cb; | ||
1299 | float se, mel_, dmel; | ||
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 | |||
1319 | } | ||
1320 | } | ||
1321 | |||
1322 | |||
1323 | /*---------------------------------------------------------------------------*\ | ||
1324 | |||
1325 | FUNCTION....: decode_mels_scalar() | ||
1326 | AUTHOR......: David Rowe | ||
1327 | DATE CREATED: April 2015 | ||
1328 | |||
1329 | From a vector of quantised mel indexes, returns the quantised | ||
1330 | (floating point) mels. | ||
1331 | |||
1332 | \*---------------------------------------------------------------------------*/ | ||
1333 | |||
1334 | void 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 | |||
1367 | void 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 | } | ||
1416 | |||
1417 | |||
1418 | /*---------------------------------------------------------------------------*\ | ||
1419 | |||
1420 | FUNCTION....: decode_lsps_diff_freq_vq() | ||
1421 | AUTHOR......: David Rowe | ||
1422 | DATE CREATED: 15 Nov 2011 | ||
1423 | |||
1424 | From a vector of quantised LSP indexes, returns the quantised | ||
1425 | (floating point) LSPs. | ||
1426 | |||
1427 | \*---------------------------------------------------------------------------*/ | ||
1428 | |||
1429 | void decode_lsps_diff_freq_vq(float lsp_[], int indexes[], int order) | ||
1430 | { | ||
1431 | int i,k,m; | ||
1432 | float dlsp_[order]; | ||
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 | |||
1473 | void encode_lsps_diff_time(int indexes[], | ||
1474 | float lsps[], | ||
1475 | float lsps__prev[], | ||
1476 | int order) | ||
1477 | { | ||
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 | |||
1515 | void 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 | |||
1527 | for(i=0; i<order; i++) { | ||
1528 | k = lsp_cbdt[i].k; | ||
1529 | cb = lsp_cbdt[i].cb; | ||
1530 | lsps_[i] += (PI/4000.0)*cb[indexes[i]*k]; | ||
1531 | } | ||
1532 | |||
1533 | } | ||
1534 | #endif | ||
1535 | |||
1536 | /*---------------------------------------------------------------------------*\ | ||
1537 | |||
1538 | FUNCTION....: encode_lsps_vq() | ||
1539 | AUTHOR......: David Rowe | ||
1540 | DATE CREATED: 15 Feb 2012 | ||
1541 | |||
1542 | Multi-stage VQ LSP quantiser developed by Jean-Marc Valin. | ||
1543 | |||
1544 | \*---------------------------------------------------------------------------*/ | ||
1545 | |||
1546 | void encode_lsps_vq(int *indexes, float *x, float *xq, int order) | ||
1547 | { | ||
1548 | int i, n1, n2, n3; | ||
1549 | float err[order], err2[order], err3[order]; | ||
1550 | float w[order], w2[order], w3[order]; | ||
1551 | const float *codebook1 = lsp_cbjvm[0].cb; | ||
1552 | const float *codebook2 = lsp_cbjvm[1].cb; | ||
1553 | const float *codebook3 = lsp_cbjvm[2].cb; | ||
1554 | |||
1555 | w[0] = MIN(x[0], x[1]-x[0]); | ||
1556 | for (i=1;i<order-1;i++) | ||
1557 | w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]); | ||
1558 | w[order-1] = MIN(x[order-1]-x[order-2], PI-x[order-1]); | ||
1559 | |||
1560 | compute_weights(x, w, order); | ||
1561 | |||
1562 | n1 = find_nearest(codebook1, lsp_cbjvm[0].m, x, order); | ||
1563 | |||
1564 | for (i=0;i<order;i++) | ||
1565 | { | ||
1566 | xq[i] = codebook1[order*n1+i]; | ||
1567 | err[i] = x[i] - xq[i]; | ||
1568 | } | ||
1569 | for (i=0;i<order/2;i++) | ||
1570 | { | ||
1571 | err2[i] = err[2*i]; | ||
1572 | err3[i] = err[2*i+1]; | ||
1573 | w2[i] = w[2*i]; | ||
1574 | w3[i] = w[2*i+1]; | ||
1575 | } | ||
1576 | n2 = find_nearest_weighted(codebook2, lsp_cbjvm[1].m, err2, w2, order/2); | ||
1577 | n3 = find_nearest_weighted(codebook3, lsp_cbjvm[2].m, err3, w3, order/2); | ||
1578 | |||
1579 | indexes[0] = n1; | ||
1580 | indexes[1] = n2; | ||
1581 | indexes[2] = n3; | ||
1582 | } | ||
1583 | |||
1584 | |||
1585 | /*---------------------------------------------------------------------------*\ | ||
1586 | |||
1587 | FUNCTION....: decode_lsps_vq() | ||
1588 | AUTHOR......: David Rowe | ||
1589 | DATE CREATED: 15 Feb 2012 | ||
1590 | |||
1591 | \*---------------------------------------------------------------------------*/ | ||
1592 | |||
1593 | void decode_lsps_vq(int *indexes, float *xq, int order, int stages) | ||
1594 | { | ||
1595 | int i, n1, n2, n3; | ||
1596 | const float *codebook1 = lsp_cbjvm[0].cb; | ||
1597 | const float *codebook2 = lsp_cbjvm[1].cb; | ||
1598 | const float *codebook3 = lsp_cbjvm[2].cb; | ||
1599 | |||
1600 | n1 = indexes[0]; | ||
1601 | n2 = indexes[1]; | ||
1602 | n3 = indexes[2]; | ||
1603 | |||
1604 | for (i=0;i<order;i++) { | ||
1605 | xq[i] = codebook1[order*n1+i]; | ||
1606 | } | ||
1607 | |||
1608 | if (stages != 1) { | ||
1609 | for (i=0;i<order/2;i++) { | ||
1610 | xq[2*i] += codebook2[order*n2/2+i]; | ||
1611 | xq[2*i+1] += codebook3[order*n3/2+i]; | ||
1612 | } | ||
1613 | } | ||
1614 | |||
1615 | } | ||
1616 | |||
1617 | |||
1618 | /*---------------------------------------------------------------------------*\ | ||
1619 | |||
1620 | FUNCTION....: bw_expand_lsps() | ||
1621 | AUTHOR......: David Rowe | ||
1622 | DATE CREATED: 22/8/2010 | ||
1623 | |||
1624 | Applies Bandwidth Expansion (BW) to a vector of LSPs. Prevents any | ||
1625 | two LSPs getting too close together after quantisation. We know | ||
1626 | from experiment that LSP quantisation errors < 12.5Hz (25Hz step | ||
1627 | size) are inaudible so we use that as the minimum LSP separation. | ||
1628 | |||
1629 | \*---------------------------------------------------------------------------*/ | ||
1630 | |||
1631 | void bw_expand_lsps(float lsp[], int order, float min_sep_low, float min_sep_high) | ||
1632 | { | ||
1633 | 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 | |||
1653 | void 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 | |||
1664 | } | ||
1665 | |||
1666 | /* As quantiser gaps increased, larger BW expansion was required | ||
1667 | to prevent twinkly noises. This may need more experiment for | ||
1668 | different quanstisers. | ||
1669 | */ | ||
1670 | |||
1671 | for(i=4; i<order; i++) { | ||
1672 | if (lsp[i] - lsp[i-1] < 200.0*(PI/4000.0)) | ||
1673 | lsp[i] = lsp[i-1] + 200.0*(PI/4000.0); | ||
1674 | } | ||
1675 | } | ||
1676 | |||
1677 | /*---------------------------------------------------------------------------*\ | ||
1678 | |||
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 | |||
1698 | void 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 | |||
1733 | /* quantise to 100Hz steps */ | ||
1734 | |||
1735 | step = 100; | ||
1736 | for(i=4; i<10; i++) { | ||
1737 | lsp_hz = lsps[i]*4000.0/PI; | ||
1738 | lsp_hz = floorf(lsp_hz/step + 0.5)*step; | ||
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 | |||
1744 | } | ||
1745 | } | ||
1746 | } | ||
1747 | |||
1748 | |||
1749 | /*---------------------------------------------------------------------------*\ | ||
1750 | |||
1751 | FUNCTION....: apply_lpc_correction() | ||
1752 | AUTHOR......: David Rowe | ||
1753 | DATE CREATED: 22/8/2010 | ||
1754 | |||
1755 | Apply first harmonic LPC correction at decoder. This helps improve | ||
1756 | low pitch males after LPC modelling, like hts1a and morig. | ||
1757 | |||
1758 | \*---------------------------------------------------------------------------*/ | ||
1759 | |||
1760 | void apply_lpc_correction(MODEL *model) | ||
1761 | { | ||
1762 | if (model->Wo < (PI*150.0/4000)) { | ||
1763 | model->A[1] *= 0.032; | ||
1764 | } | ||
1765 | } | ||
1766 | |||
1767 | /*---------------------------------------------------------------------------*\ | ||
1768 | |||
1769 | FUNCTION....: encode_energy() | ||
1770 | AUTHOR......: David Rowe | ||
1771 | DATE CREATED: 22/8/2010 | ||
1772 | |||
1773 | Encodes LPC energy using an E_LEVELS quantiser. | ||
1774 | |||
1775 | \*---------------------------------------------------------------------------*/ | ||
1776 | |||
1777 | int encode_energy(float e, int bits) | ||
1778 | { | ||
1779 | int index, e_levels = 1<<bits; | ||
1780 | float e_min = E_MIN_DB; | ||
1781 | float e_max = E_MAX_DB; | ||
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 | |||
1803 | float 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 | |||
1811 | step = (e_max - e_min)/e_levels; | ||
1812 | e = e_min + step*(index); | ||
1813 | e = POW10F(e/10.0); | ||
1814 | |||
1815 | return e; | ||
1816 | } | ||
1817 | |||
1818 | #ifdef NOT_USED | ||
1819 | /*---------------------------------------------------------------------------*\ | ||
1820 | |||
1821 | FUNCTION....: decode_amplitudes() | ||
1822 | AUTHOR......: David Rowe | ||
1823 | DATE CREATED: 22/8/2010 | ||
1824 | |||
1825 | Given the amplitude quantiser indexes recovers the harmonic | ||
1826 | amplitudes. | ||
1827 | |||
1828 | \*---------------------------------------------------------------------------*/ | ||
1829 | |||
1830 | float decode_amplitudes(codec2_fft_cfg fft_fwd_cfg, | ||
1831 | MODEL *model, | ||
1832 | float ak[], | ||
1833 | int lsp_indexes[], | ||
1834 | int energy_index, | ||
1835 | float lsps[], | ||
1836 | float *e | ||
1837 | ) | ||
1838 | { | ||
1839 | float snr; | ||
1840 | |||
1841 | decode_lsps_scalar(lsps, lsp_indexes, LPC_ORD); | ||
1842 | bw_expand_lsps(lsps, LPC_ORD); | ||
1843 | lsp_to_lpc(lsps, ak, LPC_ORD); | ||
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 | |||
1848 | return snr; | ||
1849 | } | ||
1850 | #endif | ||
1851 | |||
1852 | static float ge_coeff[2] = {0.8, 0.9}; | ||
1853 | |||
1854 | void compute_weights2(const float *x, const float *xp, float *w) | ||
1855 | { | ||
1856 | w[0] = 30; | ||
1857 | w[1] = 1; | ||
1858 | if (x[1]<0) | ||
1859 | { | ||
1860 | w[0] *= .6; | ||
1861 | w[1] *= .3; | ||
1862 | } | ||
1863 | if (x[1]<-10) | ||
1864 | { | ||
1865 | w[0] *= .3; | ||
1866 | w[1] *= .3; | ||
1867 | } | ||
1868 | /* Higher weight if pitch is stable */ | ||
1869 | if (fabsf(x[0]-xp[0])<.2) | ||
1870 | { | ||
1871 | w[0] *= 2; | ||
1872 | w[1] *= 1.5; | ||
1873 | } else if (fabsf(x[0]-xp[0])>.5) /* Lower if not stable */ | ||
1874 | { | ||
1875 | w[0] *= .5; | ||
1876 | } | ||
1877 | |||
1878 | /* Lower weight for low energy */ | ||
1879 | if (x[1] < xp[1]-10) | ||
1880 | { | ||
1881 | w[1] *= .5; | ||
1882 | } | ||
1883 | if (x[1] < xp[1]-20) | ||
1884 | { | ||
1885 | w[1] *= .5; | ||
1886 | } | ||
1887 | |||
1888 | //w[0] = 30; | ||
1889 | //w[1] = 1; | ||
1890 | |||
1891 | /* Square the weights because it's applied on the squared error */ | ||
1892 | w[0] *= w[0]; | ||
1893 | w[1] *= w[1]; | ||
1894 | |||
1895 | } | ||
1896 | |||
1897 | /*---------------------------------------------------------------------------*\ | ||
1898 | |||
1899 | FUNCTION....: quantise_WoE() | ||
1900 | AUTHOR......: Jean-Marc Valin & David Rowe | ||
1901 | DATE CREATED: 29 Feb 2012 | ||
1902 | |||
1903 | Experimental joint Wo and LPC energy vector quantiser developed by | ||
1904 | Jean-Marc Valin. Exploits correlations between the difference in | ||
1905 | the log pitch and log energy from frame to frame. For example | ||
1906 | both the pitch and energy tend to only change by small amounts | ||
1907 | during voiced speech, however it is important that these changes be | ||
1908 | coded carefully. During unvoiced speech they both change a lot but | ||
1909 | the ear is less sensitve to errors so coarser quantisation is OK. | ||
1910 | |||
1911 | 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 | ||
1913 | values is close to way the ear senses errors. | ||
1914 | |||
1915 | See http://jmspeex.livejournal.com/10446.html | ||
1916 | |||
1917 | \*---------------------------------------------------------------------------*/ | ||
1918 | |||
1919 | void quantise_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[]) | ||
1920 | { | ||
1921 | int i, n1; | ||
1922 | float x[2]; | ||
1923 | float err[2]; | ||
1924 | float w[2]; | ||
1925 | const float *codebook1 = ge_cb[0].cb; | ||
1926 | int nb_entries = ge_cb[0].m; | ||
1927 | int ndim = ge_cb[0].k; | ||
1928 | float Wo_min = c2const->Wo_min; | ||
1929 | float Wo_max = c2const->Wo_max; | ||
1930 | float Fs = c2const->Fs; | ||
1931 | |||
1932 | /* VQ is only trained for Fs = 8000 Hz */ | ||
1933 | |||
1934 | assert(Fs == 8000); | ||
1935 | |||
1936 | x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2); | ||
1937 | x[1] = 10.0*log10f(1e-4 + *e); | ||
1938 | |||
1939 | compute_weights2(x, xq, w); | ||
1940 | for (i=0;i<ndim;i++) | ||
1941 | err[i] = x[i]-ge_coeff[i]*xq[i]; | ||
1942 | n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim); | ||
1943 | |||
1944 | for (i=0;i<ndim;i++) | ||
1945 | { | ||
1946 | xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i]; | ||
1947 | err[i] -= codebook1[ndim*n1+i]; | ||
1948 | } | ||
1949 | |||
1950 | /* | ||
1951 | x = log2(4000*Wo/(PI*50)); | ||
1952 | 2^x = 4000*Wo/(PI*50) | ||
1953 | Wo = (2^x)*(PI*50)/4000; | ||
1954 | */ | ||
1955 | |||
1956 | model->Wo = powf(2.0, xq[0])*(PI*50.0)/4000.0; | ||
1957 | |||
1958 | /* bit errors can make us go out of range leading to all sorts of | ||
1959 | probs like seg faults */ | ||
1960 | |||
1961 | if (model->Wo > Wo_max) model->Wo = Wo_max; | ||
1962 | if (model->Wo < Wo_min) model->Wo = Wo_min; | ||
1963 | |||
1964 | model->L = PI/model->Wo; /* if we quantise Wo re-compute L */ | ||
1965 | |||
1966 | *e = POW10F(xq[1]/10.0); | ||
1967 | } | ||
1968 | |||
1969 | /*---------------------------------------------------------------------------*\ | ||
1970 | |||
1971 | FUNCTION....: encode_WoE() | ||
1972 | AUTHOR......: Jean-Marc Valin & David Rowe | ||
1973 | DATE CREATED: 11 May 2012 | ||
1974 | |||
1975 | Joint Wo and LPC energy vector quantiser developed my Jean-Marc | ||
1976 | Valin. Returns index, and updated states xq[]. | ||
1977 | |||
1978 | \*---------------------------------------------------------------------------*/ | ||
1979 | |||
1980 | int encode_WoE(MODEL *model, float e, float xq[]) | ||
1981 | { | ||
1982 | int i, n1; | ||
1983 | float x[2]; | ||
1984 | float err[2]; | ||
1985 | float w[2]; | ||
1986 | const float *codebook1 = ge_cb[0].cb; | ||
1987 | int nb_entries = ge_cb[0].m; | ||
1988 | int ndim = ge_cb[0].k; | ||
1989 | |||
1990 | assert((1<<WO_E_BITS) == nb_entries); | ||
1991 | |||
1992 | if (e < 0.0) e = 0; /* occasional small negative energies due LPC round off I guess */ | ||
1993 | |||
1994 | x[0] = log10f((model->Wo/PI)*4000.0/50.0)/log10f(2); | ||
1995 | x[1] = 10.0*log10f(1e-4 + e); | ||
1996 | |||
1997 | compute_weights2(x, xq, w); | ||
1998 | for (i=0;i<ndim;i++) | ||
1999 | err[i] = x[i]-ge_coeff[i]*xq[i]; | ||
2000 | n1 = find_nearest_weighted(codebook1, nb_entries, err, w, ndim); | ||
2001 | |||
2002 | for (i=0;i<ndim;i++) | ||
2003 | { | ||
2004 | xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i]; | ||
2005 | err[i] -= codebook1[ndim*n1+i]; | ||
2006 | } | ||
2007 | |||
2008 | //printf("enc: %f %f (%f)(%f) \n", xq[0], xq[1], e, 10.0*log10(1e-4 + e)); | ||
2009 | return n1; | ||
2010 | } | ||
2011 | |||
2012 | |||
2013 | /*---------------------------------------------------------------------------*\ | ||
2014 | |||
2015 | FUNCTION....: decode_WoE() | ||
2016 | AUTHOR......: Jean-Marc Valin & David Rowe | ||
2017 | DATE CREATED: 11 May 2012 | ||
2018 | |||
2019 | Joint Wo and LPC energy vector quantiser developed my Jean-Marc | ||
2020 | Valin. Given index and states xq[], returns Wo & E, and updates | ||
2021 | states xq[]. | ||
2022 | |||
2023 | \*---------------------------------------------------------------------------*/ | ||
2024 | |||
2025 | void decode_WoE(C2CONST *c2const, MODEL *model, float *e, float xq[], int n1) | ||
2026 | { | ||
2027 | int i; | ||
2028 | const float *codebook1 = ge_cb[0].cb; | ||
2029 | int ndim = ge_cb[0].k; | ||
2030 | float Wo_min = c2const->Wo_min; | ||
2031 | float Wo_max = c2const->Wo_max; | ||
2032 | |||
2033 | for (i=0;i<ndim;i++) | ||
2034 | { | ||
2035 | xq[i] = ge_coeff[i]*xq[i] + codebook1[ndim*n1+i]; | ||
2036 | } | ||
2037 | |||
2038 | //printf("dec: %f %f\n", xq[0], xq[1]); | ||
2039 | model->Wo = powf(2.0, xq[0])*(PI*50.0)/4000.0; | ||
2040 | |||
2041 | /* bit errors can make us go out of range leading to all sorts of | ||
2042 | probs like seg faults */ | ||
2043 | |||
2044 | if (model->Wo > Wo_max) model->Wo = Wo_max; | ||
2045 | if (model->Wo < Wo_min) model->Wo = Wo_min; | ||
2046 | |||
2047 | model->L = PI/model->Wo; /* if we quantise Wo re-compute L */ | ||
2048 | |||
2049 | *e = POW10F(xq[1]/10.0); | ||
2050 | } | ||
2051 | |||