summaryrefslogtreecommitdiff
path: root/newamp1.c
diff options
context:
space:
mode:
Diffstat (limited to 'newamp1.c')
-rw-r--r--newamp1.c770
1 files changed, 394 insertions, 376 deletions
diff --git a/newamp1.c b/newamp1.c
index 8980ac6..3ba2de0 100644
--- a/newamp1.c
+++ b/newamp1.c
@@ -28,19 +28,18 @@
28 28
29*/ 29*/
30 30
31#include "newamp1.h"
32
31#include <assert.h> 33#include <assert.h>
34#include <math.h>
32#include <stdio.h> 35#include <stdio.h>
33#include <stdlib.h> 36#include <stdlib.h>
34#include <string.h> 37#include <string.h>
35#include <math.h>
36 38
37#include "defines.h" 39#include "defines.h"
40#include "mbest.h"
38#include "phase.h" 41#include "phase.h"
39#include "quantise.h" 42#include "quantise.h"
40#include "mbest.h"
41#include "newamp1.h"
42
43#define NEWAMP1_VQ_MBEST_DEPTH 5 /* how many candidates we keep for each stage of mbest search */
44 43
45/*---------------------------------------------------------------------------*\ 44/*---------------------------------------------------------------------------*\
46 45
@@ -48,39 +47,44 @@
48 AUTHOR......: David Rowe 47 AUTHOR......: David Rowe
49 DATE CREATED: Jan 2017 48 DATE CREATED: Jan 2017
50 49
51 General 2nd order parabolic interpolator. Used splines orginally, 50 General 2nd order parabolic interpolator. Used splines originally,
52 but this is much simpler and we don't need much accuracy. Given two 51 but this is much simpler and we don't need much accuracy. Given two
53 vectors of points xp and yp, find interpolated values y at points x. 52 vectors of points xp and yp, find interpolated values y at points x.
54 53
55\*---------------------------------------------------------------------------*/ 54\*---------------------------------------------------------------------------*/
56 55
57void interp_para(float y[], float xp[], float yp[], int np, float x[], int n) 56void interp_para(float y[], float xp[], float yp[], int np, float x[], int n) {
58{ 57 assert(np >= 3);
59 assert(np >= 3);
60 58
61 int k,i; 59 int k, i;
62 float xi, x1, y1, x2, y2, x3, y3, a, b; 60 float xi, x1, y1, x2, y2, x3, y3, a, b;
63 61
64 k = 0; 62 k = 0;
65 for (i=0; i<n; i++) { 63 for (i = 0; i < n; i++) {
66 xi = x[i]; 64 xi = x[i];
67 65
68 /* k is index into xp of where we start 3 points used to form parabola */ 66 /* k is index into xp of where we start 3 points used to form parabola */
69 67
70 while ((xp[k+1] < xi) && (k < (np-3))) 68 while ((xp[k + 1] < xi) && (k < (np - 3))) k++;
71 k++;
72
73 x1 = xp[k]; y1 = yp[k]; x2 = xp[k+1]; y2 = yp[k+1]; x3 = xp[k+2]; y3 = yp[k+2];
74 69
75 //printf("k: %d np: %d i: %d xi: %f x1: %f y1: %f\n", k, np, i, xi, x1, y1); 70 x1 = xp[k];
71 y1 = yp[k];
72 x2 = xp[k + 1];
73 y2 = yp[k + 1];
74 x3 = xp[k + 2];
75 y3 = yp[k + 2];
76 76
77 a = ((y3-y2)/(x3-x2)-(y2-y1)/(x2-x1))/(x3-x1); 77 // printf("k: %d np: %d i: %d xi: %f x1: %f y1: %f\n", k, np, i, xi, x1,
78 b = ((y3-y2)/(x3-x2)*(x2-x1)+(y2-y1)/(x2-x1)*(x3-x2))/(x3-x1); 78 // y1);
79
80 y[i] = a*(xi-x2)*(xi-x2) + b*(xi-x2) + y2;
81 }
82}
83 79
80 a = ((y3 - y2) / (x3 - x2) - (y2 - y1) / (x2 - x1)) / (x3 - x1);
81 b = ((y3 - y2) / (x3 - x2) * (x2 - x1) +
82 (y2 - y1) / (x2 - x1) * (x3 - x2)) /
83 (x3 - x1);
84
85 y[i] = a * (xi - x2) * (xi - x2) + b * (xi - x2) + y2;
86 }
87}
84 88
85/*---------------------------------------------------------------------------*\ 89/*---------------------------------------------------------------------------*\
86 90
@@ -94,24 +98,23 @@ void interp_para(float y[], float xp[], float yp[], int np, float x[], int n)
94\*---------------------------------------------------------------------------*/ 98\*---------------------------------------------------------------------------*/
95 99
96float ftomel(float fHz) { 100float ftomel(float fHz) {
97 float mel = floorf(2595.0*log10f(1.0 + fHz/700.0)+0.5); 101 float mel = floorf(2595.0 * log10f(1.0 + fHz / 700.0) + 0.5);
98 return mel; 102 return mel;
99} 103}
100 104
101void mel_sample_freqs_kHz(float rate_K_sample_freqs_kHz[], int K, float mel_start, float mel_end) 105void mel_sample_freqs_kHz(float rate_K_sample_freqs_kHz[], int K,
102{ 106 float mel_start, float mel_end) {
103 float step = (mel_end-mel_start)/(K-1); 107 float step = (mel_end - mel_start) / (K - 1);
104 float mel; 108 float mel;
105 int k; 109 int k;
106 110
107 mel = mel_start; 111 mel = mel_start;
108 for (k=0; k<K; k++) { 112 for (k = 0; k < K; k++) {
109 rate_K_sample_freqs_kHz[k] = 0.7*(POW10F(mel/2595.0) - 1.0); 113 rate_K_sample_freqs_kHz[k] = 0.7 * (POW10F(mel / 2595.0) - 1.0);
110 mel += step; 114 mel += step;
111 } 115 }
112} 116}
113 117
114
115/*---------------------------------------------------------------------------*\ 118/*---------------------------------------------------------------------------*\
116 119
117 FUNCTION....: resample_const_rate_f() 120 FUNCTION....: resample_const_rate_f()
@@ -122,35 +125,36 @@ void mel_sample_freqs_kHz(float rate_K_sample_freqs_kHz[], int K, float mel_star
122 125
123\*---------------------------------------------------------------------------*/ 126\*---------------------------------------------------------------------------*/
124 127
125void resample_const_rate_f(C2CONST *c2const, MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], int K) 128void resample_const_rate_f(C2CONST *c2const, MODEL *model, float rate_K_vec[],
126{ 129 float rate_K_sample_freqs_kHz[], int K) {
127 int m; 130 int m;
128 float AmdB[MAX_AMP+1], rate_L_sample_freqs_kHz[MAX_AMP+1], AmdB_peak; 131 float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1], AmdB_peak;
129 132
130 /* convert rate L=pi/Wo amplitude samples to fixed rate K */ 133 /* convert rate L=pi/Wo amplitude samples to fixed rate K */
131 134
132 AmdB_peak = -100.0; 135 AmdB_peak = -100.0;
133 for(m=1; m<=model->L; m++) { 136 for (m = 1; m <= model->L; m++) {
134 AmdB[m] = 20.0*log10f(model->A[m]+1E-16); 137 AmdB[m] = 20.0 * log10f(model->A[m] + 1E-16);
135 if (AmdB[m] > AmdB_peak) { 138 if (AmdB[m] > AmdB_peak) {
136 AmdB_peak = AmdB[m]; 139 AmdB_peak = AmdB[m];
137 }
138 rate_L_sample_freqs_kHz[m] = m*model->Wo*(c2const->Fs/2000.0)/M_PI;
139 //printf("m: %d AmdB: %f AmdB_peak: %f sf: %f\n", m, AmdB[m], AmdB_peak, rate_L_sample_freqs_kHz[m]);
140 } 140 }
141 141 rate_L_sample_freqs_kHz[m] = m * model->Wo * (c2const->Fs / 2000.0) / M_PI;
142 /* clip between peak and peak -50dB, to reduce dynamic range */ 142 // printf("m: %d AmdB: %f AmdB_peak: %f sf: %f\n", m, AmdB[m], AmdB_peak,
143 // rate_L_sample_freqs_kHz[m]);
144 }
145
146 /* clip between peak and peak -50dB, to reduce dynamic range */
143 147
144 for(m=1; m<=model->L; m++) { 148 for (m = 1; m <= model->L; m++) {
145 if (AmdB[m] < (AmdB_peak-50.0)) { 149 if (AmdB[m] < (AmdB_peak - 50.0)) {
146 AmdB[m] = AmdB_peak-50.0; 150 AmdB[m] = AmdB_peak - 50.0;
147 }
148 } 151 }
152 }
149 153
150 interp_para(rate_K_vec, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, rate_K_sample_freqs_kHz, K); 154 interp_para(rate_K_vec, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L,
155 rate_K_sample_freqs_kHz, K);
151} 156}
152 157
153
154/*---------------------------------------------------------------------------*\ 158/*---------------------------------------------------------------------------*\
155 159
156 FUNCTION....: rate_K_mbest_encode 160 FUNCTION....: rate_K_mbest_encode
@@ -161,64 +165,55 @@ void resample_const_rate_f(C2CONST *c2const, MODEL *model, float rate_K_vec[], f
161 165
162\*---------------------------------------------------------------------------*/ 166\*---------------------------------------------------------------------------*/
163 167
164float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest_entries) 168float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim,
165{ 169 int mbest_entries) {
166 int i, j, n1, n2; 170 int i, j, n1, n2;
167 const float *codebook1 = newamp1vq_cb[0].cb; 171 const float *codebook1 = newamp1vq_cb[0].cb;
168 const float *codebook2 = newamp1vq_cb[1].cb; 172 const float *codebook2 = newamp1vq_cb[1].cb;
169 struct MBEST *mbest_stage1, *mbest_stage2; 173 struct MBEST *mbest_stage1, *mbest_stage2;
170 float target[ndim]; 174 float target[ndim];
171 float w[ndim]; 175 int index[MBEST_STAGES];
172 int index[MBEST_STAGES];
173 float mse, tmp; 176 float mse, tmp;
174 177
175 /* codebook is compiled for a fixed K */ 178 /* codebook is compiled for a fixed K */
176 179
177 assert(ndim == newamp1vq_cb[0].k); 180 assert(ndim == newamp1vq_cb[0].k);
178 181
179 /* equal weights, could be argued mel freq axis gives freq dep weighting */
180
181 for(i=0; i<ndim; i++)
182 w[i] = 1.0;
183
184 mbest_stage1 = mbest_create(mbest_entries); 182 mbest_stage1 = mbest_create(mbest_entries);
185 mbest_stage2 = mbest_create(mbest_entries); 183 mbest_stage2 = mbest_create(mbest_entries);
186 for(i=0; i<MBEST_STAGES; i++) 184 for (i = 0; i < MBEST_STAGES; i++) index[i] = 0;
187 index[i] = 0;
188 185
189 /* Stage 1 */ 186 /* Stage 1 */
190 187
191 mbest_search(codebook1, x, w, ndim, newamp1vq_cb[0].m, mbest_stage1, index); 188 mbest_search(codebook1, x, ndim, newamp1vq_cb[0].m, mbest_stage1, index);
192 MBEST_PRINT("Stage 1:", mbest_stage1);
193 189
194 /* Stage 2 */ 190 /* Stage 2 */
195 191
196 for (j=0; j<mbest_entries; j++) { 192 for (j = 0; j < mbest_entries; j++) {
197 index[1] = n1 = mbest_stage1->list[j].index[0]; 193 index[1] = n1 = mbest_stage1->list[j].index[0];
198 for(i=0; i<ndim; i++) 194 for (i = 0; i < ndim; i++) target[i] = x[i] - codebook1[ndim * n1 + i];
199 target[i] = x[i] - codebook1[ndim*n1+i]; 195 mbest_search(codebook2, target, ndim, newamp1vq_cb[1].m, mbest_stage2,
200 mbest_search(codebook2, target, w, ndim, newamp1vq_cb[1].m, mbest_stage2, index); 196 index);
201 } 197 }
202 MBEST_PRINT("Stage 2:", mbest_stage2);
203 198
204 n1 = mbest_stage2->list[0].index[1]; 199 n1 = mbest_stage2->list[0].index[1];
205 n2 = mbest_stage2->list[0].index[0]; 200 n2 = mbest_stage2->list[0].index[0];
206 mse = 0.0; 201 mse = 0.0;
207 for (i=0;i<ndim;i++) { 202 for (i = 0; i < ndim; i++) {
208 tmp = codebook1[ndim*n1+i] + codebook2[ndim*n2+i]; 203 tmp = codebook1[ndim * n1 + i] + codebook2[ndim * n2 + i];
209 mse += (x[i]-tmp)*(x[i]-tmp); 204 mse += (x[i] - tmp) * (x[i] - tmp);
210 xq[i] = tmp; 205 xq[i] = tmp;
211 } 206 }
212 207
213 mbest_destroy(mbest_stage1); 208 mbest_destroy(mbest_stage1);
214 mbest_destroy(mbest_stage2); 209 mbest_destroy(mbest_stage2);
215 210
216 indexes[0] = n1; indexes[1] = n2; 211 indexes[0] = n1;
212 indexes[1] = n2;
217 213
218 return mse; 214 return mse;
219} 215}
220 216
221
222/*---------------------------------------------------------------------------*\ 217/*---------------------------------------------------------------------------*\
223 218
224 FUNCTION....: post_filter 219 FUNCTION....: post_filter
@@ -226,7 +221,7 @@ float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest
226 DATE CREATED: Jan 2017 221 DATE CREATED: Jan 2017
227 222
228 Post Filter, has a big impact on speech quality after VQ. When used 223 Post Filter, has a big impact on speech quality after VQ. When used
229 on a mean removed rate K vector, it raises formants, and supresses 224 on a mean removed rate K vector, it raises formants, and suppresses
230 anti-formants. As it manipulates amplitudes, we normalise energy to 225 anti-formants. As it manipulates amplitudes, we normalise energy to
231 prevent clipping or large level variations. pf_gain of 1.2 to 1.5 226 prevent clipping or large level variations. pf_gain of 1.2 to 1.5
232 (dB) seems to work OK. Good area for further investigations and 227 (dB) seems to work OK. Good area for further investigations and
@@ -234,37 +229,36 @@ float rate_K_mbest_encode(int *indexes, float *x, float *xq, int ndim, int mbest
234 229
235\*---------------------------------------------------------------------------*/ 230\*---------------------------------------------------------------------------*/
236 231
237void post_filter_newamp1(float vec[], float sample_freq_kHz[], int K, float pf_gain) 232void post_filter_newamp1(float vec[], float sample_freq_kHz[], int K,
238{ 233 float pf_gain) {
239 int k; 234 int k;
240 235
241 /* 236 /*
242 vec is rate K vector describing spectrum of current frame lets 237 vec is rate K vector describing spectrum of current frame lets
243 pre-emp before applying PF. 20dB/dec over 300Hz. Postfilter 238 pre-emp before applying PF. 20dB/dec over 300Hz. Postfilter
244 affects energy of frame so we measure energy before and after 239 affects energy of frame so we measure energy before and after
245 and normalise. Plenty of room for experiment here as well. 240 and normalise. Plenty of room for experimentation here.
246 */ 241 */
247 242
248 float pre[K]; 243 float pre[K];
249 float e_before = 0.0; 244 float e_before = 0.0;
250 float e_after = 0.0; 245 float e_after = 0.0;
251 for(k=0; k<K; k++) { 246 for (k = 0; k < K; k++) {
252 pre[k] = 20.0*log10f(sample_freq_kHz[k]/0.3); 247 pre[k] = 20.0 * log10f(sample_freq_kHz[k] / 0.3);
253 vec[k] += pre[k]; 248 vec[k] += pre[k];
254 e_before += POW10F(vec[k]/10.0); 249 e_before += POW10F(vec[k] / 10.0);
255 vec[k] *= pf_gain; 250 vec[k] *= pf_gain;
256 e_after += POW10F(vec[k]/10.0); 251 e_after += POW10F(vec[k] / 10.0);
257 } 252 }
258 253
259 float gain = e_after/e_before; 254 float gain = e_after / e_before;
260 float gaindB = 10*log10f(gain); 255 float gaindB = 10 * log10f(gain);
261
262 for(k=0; k<K; k++) {
263 vec[k] -= gaindB;
264 vec[k] -= pre[k];
265 }
266}
267 256
257 for (k = 0; k < K; k++) {
258 vec[k] -= gaindB;
259 vec[k] -= pre[k];
260 }
261}
268 262
269/*---------------------------------------------------------------------------*\ 263/*---------------------------------------------------------------------------*\
270 264
@@ -273,49 +267,46 @@ void post_filter_newamp1(float vec[], float sample_freq_kHz[], int K, float pf_g
273 DATE CREATED: Jan 2017 267 DATE CREATED: Jan 2017
274 268
275 Decoder side interpolation of Wo and voicing, to go from 25 Hz 269 Decoder side interpolation of Wo and voicing, to go from 25 Hz
276 sample rate used over channle to 100Hz internal sample rate of Codec 2. 270 sample rate used over channel to 100Hz internal sample rate of Codec 2.
277 271
278\*---------------------------------------------------------------------------*/ 272\*---------------------------------------------------------------------------*/
279 273
280void interp_Wo_v(float Wo_[], int L_[], int voicing_[], float Wo1, float Wo2, int voicing1, int voicing2) 274void interp_Wo_v(float Wo_[], int L_[], int voicing_[], float Wo1, float Wo2,
281{ 275 int voicing1, int voicing2) {
282 int i; 276 int i;
283 int M = 4; /* interpolation rate */ 277 int M = 4; /* interpolation rate */
284 278
285 for(i=0; i<M; i++) 279 for (i = 0; i < M; i++) voicing_[i] = 0;
286 voicing_[i] = 0;
287 280
288 if (!voicing1 && !voicing2) { 281 if (!voicing1 && !voicing2) {
289 for(i=0; i<M; i++) 282 for (i = 0; i < M; i++) Wo_[i] = 2.0 * M_PI / 100.0;
290 Wo_[i] = 2.0*M_PI/100.0; 283 }
291 }
292 284
293 if (voicing1 && !voicing2) { 285 if (voicing1 && !voicing2) {
294 Wo_[0] = Wo_[1] = Wo1; 286 Wo_[0] = Wo_[1] = Wo1;
295 Wo_[2] = Wo_[3] = 2.0*M_PI/100.0; 287 Wo_[2] = Wo_[3] = 2.0 * M_PI / 100.0;
296 voicing_[0] = voicing_[1] = 1; 288 voicing_[0] = voicing_[1] = 1;
297 } 289 }
298 290
299 if (!voicing1 && voicing2) { 291 if (!voicing1 && voicing2) {
300 Wo_[0] = Wo_[1] = 2.0*M_PI/100.0; 292 Wo_[0] = Wo_[1] = 2.0 * M_PI / 100.0;
301 Wo_[2] = Wo_[3] = Wo2; 293 Wo_[2] = Wo_[3] = Wo2;
302 voicing_[2] = voicing_[3] = 1; 294 voicing_[2] = voicing_[3] = 1;
303 } 295 }
304 296
305 if (voicing1 && voicing2) { 297 if (voicing1 && voicing2) {
306 float c; 298 float c;
307 for(i=0,c=1.0; i<M; i++,c-=1.0/M) { 299 for (i = 0, c = 1.0; i < M; i++, c -= 1.0 / M) {
308 Wo_[i] = Wo1*c + Wo2*(1.0-c); 300 Wo_[i] = Wo1 * c + Wo2 * (1.0 - c);
309 voicing_[i] = 1; 301 voicing_[i] = 1;
310 }
311 } 302 }
303 }
312 304
313 for(i=0; i<M; i++) { 305 for (i = 0; i < M; i++) {
314 L_[i] = floorf(M_PI/Wo_[i]); 306 L_[i] = floorf(M_PI / Wo_[i]);
315 } 307 }
316} 308}
317 309
318
319/*---------------------------------------------------------------------------*\ 310/*---------------------------------------------------------------------------*\
320 311
321 FUNCTION....: resample_rate_L 312 FUNCTION....: resample_rate_L
@@ -326,36 +317,37 @@ void interp_Wo_v(float Wo_[], int L_[], int voicing_[], float Wo1, float Wo2, in
326 317
327\*---------------------------------------------------------------------------*/ 318\*---------------------------------------------------------------------------*/
328 319
329void resample_rate_L(C2CONST *c2const, MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], int K) 320void resample_rate_L(C2CONST *c2const, MODEL *model, float rate_K_vec[],
330{ 321 float rate_K_sample_freqs_kHz[], int K) {
331 float rate_K_vec_term[K+2], rate_K_sample_freqs_kHz_term[K+2]; 322 float rate_K_vec_term[K + 2], rate_K_sample_freqs_kHz_term[K + 2];
332 float AmdB[MAX_AMP+1], rate_L_sample_freqs_kHz[MAX_AMP+1]; 323 float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1];
333 int m,k; 324 int m, k;
334
335 /* terminate either end of the rate K vecs with 0dB points */
336
337 rate_K_vec_term[0] = rate_K_vec_term[K+1] = 0.0;
338 rate_K_sample_freqs_kHz_term[0] = 0.0;
339 rate_K_sample_freqs_kHz_term[K+1] = 4.0;
340
341 for(k=0; k<K; k++) {
342 rate_K_vec_term[k+1] = rate_K_vec[k];
343 rate_K_sample_freqs_kHz_term[k+1] = rate_K_sample_freqs_kHz[k];
344
345 //printf("k: %d f: %f rate_K: %f\n", k, rate_K_sample_freqs_kHz[k], rate_K_vec[k]);
346 }
347
348 for(m=1; m<=model->L; m++) {
349 rate_L_sample_freqs_kHz[m] = m*model->Wo*(c2const->Fs/2000.0)/M_PI;
350 }
351
352 interp_para(&AmdB[1], rate_K_sample_freqs_kHz_term, rate_K_vec_term, K+2, &rate_L_sample_freqs_kHz[1], model->L);
353 for(m=1; m<=model->L; m++) {
354 model->A[m] = POW10F(AmdB[m]/20.0);
355 // printf("m: %d f: %f AdB: %f A: %f\n", m, rate_L_sample_freqs_kHz[m], AmdB[m], model->A[m]);
356 }
357}
358 325
326 /* terminate either end of the rate K vecs with 0dB points */
327
328 rate_K_vec_term[0] = rate_K_vec_term[K + 1] = 0.0;
329 rate_K_sample_freqs_kHz_term[0] = 0.0;
330 rate_K_sample_freqs_kHz_term[K + 1] = 4.0;
331
332 for (k = 0; k < K; k++) {
333 rate_K_vec_term[k + 1] = rate_K_vec[k];
334 rate_K_sample_freqs_kHz_term[k + 1] = rate_K_sample_freqs_kHz[k];
335 // printf("k: %d f: %f rate_K: %f\n", k, rate_K_sample_freqs_kHz[k],
336 // rate_K_vec[k]);
337 }
338
339 for (m = 1; m <= model->L; m++) {
340 rate_L_sample_freqs_kHz[m] = m * model->Wo * (c2const->Fs / 2000.0) / M_PI;
341 }
342
343 interp_para(&AmdB[1], rate_K_sample_freqs_kHz_term, rate_K_vec_term, K + 2,
344 &rate_L_sample_freqs_kHz[1], model->L);
345 for (m = 1; m <= model->L; m++) {
346 model->A[m] = POW10F(AmdB[m] / 20.0);
347 // printf("m: %d f: %f AdB: %f A: %f\n", m, rate_L_sample_freqs_kHz[m],
348 // AmdB[m], model->A[m]);
349 }
350}
359 351
360/*---------------------------------------------------------------------------*\ 352/*---------------------------------------------------------------------------*\
361 353
@@ -368,34 +360,100 @@ void resample_rate_L(C2CONST *c2const, MODEL *model, float rate_K_vec[], float r
368 360
369\*---------------------------------------------------------------------------*/ 361\*---------------------------------------------------------------------------*/
370 362
371void determine_phase(C2CONST *c2const, COMP H[], MODEL *model, int Nfft, codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg) 363void determine_phase(C2CONST *c2const, COMP H[], MODEL *model, int Nfft,
372{ 364 codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg) {
373 int i,m,b; 365 int i, m, b;
374 int Ns = Nfft/2+1; 366 int Ns = Nfft / 2 + 1;
375 float Gdbfk[Ns], sample_freqs_kHz[Ns], phase[Ns]; 367 float Gdbfk[Ns], sample_freqs_kHz[Ns], phase[Ns];
376 float AmdB[MAX_AMP+1], rate_L_sample_freqs_kHz[MAX_AMP+1]; 368 float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1];
377 369
378 for(m=1; m<=model->L; m++) { 370 for (m = 1; m <= model->L; m++) {
379 assert(model->A[m] != 0.0); 371 assert(model->A[m] != 0.0);
380 AmdB[m] = 20.0*log10f(model->A[m]); 372 AmdB[m] = 20.0 * log10f(model->A[m]);
381 rate_L_sample_freqs_kHz[m] = (float)m*model->Wo*(c2const->Fs/2000.0)/M_PI; 373 rate_L_sample_freqs_kHz[m] =
382 } 374 (float)m * model->Wo * (c2const->Fs / 2000.0) / M_PI;
383 375 }
384 for(i=0; i<Ns; i++) {
385 sample_freqs_kHz[i] = (c2const->Fs/1000.0)*(float)i/Nfft;
386 }
387 376
388 interp_para(Gdbfk, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, sample_freqs_kHz, Ns); 377 for (i = 0; i < Ns; i++) {
389 mag_to_phase(phase, Gdbfk, Nfft, fwd_cfg, inv_cfg); 378 sample_freqs_kHz[i] = (c2const->Fs / 1000.0) * (float)i / Nfft;
379 }
390 380
391 for(m=1; m<=model->L; m++) { 381 interp_para(Gdbfk, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L,
392 b = floorf(0.5+m*model->Wo*Nfft/(2.0*M_PI)); 382 sample_freqs_kHz, Ns);
393 H[m].real = cosf(phase[b]); H[m].imag = sinf(phase[b]); 383 mag_to_phase(phase, Gdbfk, Nfft, fwd_cfg, inv_cfg);
394 } 384
385 for (m = 1; m <= model->L; m++) {
386 b = floorf(0.5 + m * model->Wo * Nfft / (2.0 * M_PI));
387 H[m].real = cosf(phase[b]);
388 H[m].imag = sinf(phase[b]);
389 }
395} 390}
396 391
392/*---------------------------------------------------------------------------* \
397 393
398/*---------------------------------------------------------------------------*\ 394 FUNCTION....: determine_autoc
395 AUTHOR......: David Rowe
396 DATE CREATED: April 2020
397
398 Determine autocorrelation coefficients from model params, for machine
399 learning experiments.
400
401\*---------------------------------------------------------------------------*/
402
403void determine_autoc(C2CONST *c2const, float Rk[], int order, MODEL *model,
404 int Nfft, codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg) {
405 int i, m;
406 int Ns = Nfft / 2 + 1;
407 float Gdbfk[Ns], sample_freqs_kHz[Ns];
408 float AmdB[MAX_AMP + 1], rate_L_sample_freqs_kHz[MAX_AMP + 1];
409
410 /* interpolate in the log domain */
411 for (m = 1; m <= model->L; m++) {
412 assert(model->A[m] != 0.0);
413 AmdB[m] = 20.0 * log10f(model->A[m]);
414 rate_L_sample_freqs_kHz[m] =
415 (float)m * model->Wo * (c2const->Fs / 2000.0) / M_PI;
416 }
417
418 for (i = 0; i < Ns; i++) {
419 sample_freqs_kHz[i] = (c2const->Fs / 1000.0) * (float)i / Nfft;
420 }
421
422 interp_para(Gdbfk, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L,
423 sample_freqs_kHz, Ns);
424
425 COMP S[Nfft], R[Nfft];
426
427 /* install negative frequency components, convert to mag squared of spectrum
428 */
429 S[0].real = pow(10.0, Gdbfk[0] / 10.0);
430 S[0].imag = 0.0;
431 for (i = 1; i < Ns; i++) {
432 S[i].real = S[Nfft - i].real = pow(10.0, Gdbfk[i] / 10.0);
433 S[i].imag = S[Nfft - i].imag = 0.0;
434 }
435
436 /* IDFT of mag squared is autocorrelation function */
437 codec2_fft(inv_cfg, S, R);
438 for (int k = 0; k < order + 1; k++) Rk[k] = R[k].real;
439}
440
441/* update and optionally run "front eq" equaliser on before VQ */
442void newamp1_eq(float rate_K_vec_no_mean[], float eq[], int K, int eq_en) {
443 static float ideal[] = {8, 10, 12, 14, 14, 14, 14, 14, 14, 14,
444 14, 14, 14, 14, 14, 14, 14, 14, 14, -20};
445 float gain = 0.02;
446 float update;
447
448 for (int k = 0; k < K; k++) {
449 update = rate_K_vec_no_mean[k] - ideal[k];
450 eq[k] = (1.0 - gain) * eq[k] + gain * update;
451 if (eq[k] < 0.0) eq[k] = 0.0;
452 if (eq_en) rate_K_vec_no_mean[k] -= eq[k];
453 }
454}
455
456/*---------------------------------------------------------------------------* \
399 457
400 FUNCTION....: newamp1_model_to_indexes 458 FUNCTION....: newamp1_model_to_indexes
401 AUTHOR......: David Rowe 459 AUTHOR......: David Rowe
@@ -406,78 +464,53 @@ void determine_phase(C2CONST *c2const, COMP H[], MODEL *model, int Nfft, codec2_
406 464
407\*---------------------------------------------------------------------------*/ 465\*---------------------------------------------------------------------------*/
408 466
409void newamp1_model_to_indexes(C2CONST *c2const, 467void newamp1_model_to_indexes(C2CONST *c2const, int indexes[], MODEL *model,
410 int indexes[], 468 float rate_K_vec[],
411 MODEL *model, 469 float rate_K_sample_freqs_kHz[], int K,
412 float rate_K_vec[], 470 float *mean, float rate_K_vec_no_mean[],
413 float rate_K_sample_freqs_kHz[], 471 float rate_K_vec_no_mean_[], float *se, float *eq,
414 int K, 472 int eq_en) {
415 float *mean, 473 int k;
416 float rate_K_vec_no_mean[], 474
417 float rate_K_vec_no_mean_[], 475 /* convert variable rate L to fixed rate K */
418 float *se, 476 resample_const_rate_f(c2const, model, rate_K_vec, rate_K_sample_freqs_kHz, K);
419 float *eq, 477
420 int eq_en 478 /* remove mean */
421 ) 479 float sum = 0.0;
422{ 480 for (k = 0; k < K; k++) sum += rate_K_vec[k];
423 int k; 481 *mean = sum / K;
424 482 for (k = 0; k < K; k++) rate_K_vec_no_mean[k] = rate_K_vec[k] - *mean;
425 /* convert variable rate L to fixed rate K */ 483
426 resample_const_rate_f(c2const, model, rate_K_vec, rate_K_sample_freqs_kHz, K); 484 /* update and optionally run "front eq" equaliser on before VQ */
427 485 newamp1_eq(rate_K_vec_no_mean, eq, K, eq_en);
428 /* remove mean */ 486
429 float sum = 0.0; 487 /* two stage VQ */
430 for(k=0; k<K; k++) 488 rate_K_mbest_encode(indexes, rate_K_vec_no_mean, rate_K_vec_no_mean_, K,
431 sum += rate_K_vec[k]; 489 NEWAMP1_VQ_MBEST_DEPTH);
432 *mean = sum/K; 490
433 for(k=0; k<K; k++) 491 /* running sum of squared error for variance calculation */
434 rate_K_vec_no_mean[k] = rate_K_vec[k] - *mean; 492 for (k = 0; k < K; k++)
435 493 *se += (float)pow(rate_K_vec_no_mean[k] - rate_K_vec_no_mean_[k], 2.0);
436 /* update and optionally run "front eq" equaliser on before VQ */ 494
437 static float ideal[] = {8,10,12,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,14,-20}; 495 /* scalar quantise mean (effectively the frame energy) */
438 float gain = 0.02; 496 float w[1] = {1.0};
439 float update; 497 float se_mean;
440 498 indexes[2] =
441 for(k=0; k<K; k++) { 499 quantise(newamp1_energy_cb[0].cb, mean, w, newamp1_energy_cb[0].k,
442 update = rate_K_vec_no_mean[k] - ideal[k]; 500 newamp1_energy_cb[0].m, &se_mean);
443 eq[k] = (1.0-gain)*eq[k] + gain*update; 501
444 if (eq[k] < 0.0) eq[k] = 0.0; 502 /* scalar quantise Wo. We steal the smallest Wo index to signal
445 if (eq_en) 503 an unvoiced frame */
446 rate_K_vec_no_mean[k] -= eq[k]; 504 if (model->voiced) {
447 } 505 int index = encode_log_Wo(c2const, model->Wo, 6);
448 506 if (index == 0) {
449 /* two stage VQ */ 507 index = 1;
450 rate_K_mbest_encode(indexes, rate_K_vec_no_mean, rate_K_vec_no_mean_, K, NEWAMP1_VQ_MBEST_DEPTH);
451
452 /* running sum of squared error for variance calculation */
453 for(k=0; k<K; k++)
454 *se += pow(rate_K_vec_no_mean[k]-rate_K_vec_no_mean_[k],2.0);
455
456 /* scalar quantise mean (effectively the frame energy) */
457 float w[1] = {1.0};
458 float se_mean;
459 indexes[2] = quantise(newamp1_energy_cb[0].cb,
460 mean,
461 w,
462 newamp1_energy_cb[0].k,
463 newamp1_energy_cb[0].m,
464 &se_mean);
465
466 /* scalar quantise Wo. We steal the smallest Wo index to signal
467 an unvoiced frame */
468 if (model->voiced) {
469 int index = encode_log_Wo(c2const, model->Wo, 6);
470 if (index == 0) {
471 index = 1;
472 }
473 indexes[3] = index;
474 }
475 else {
476 indexes[3] = 0;
477 } 508 }
478 509 indexes[3] = index;
479 } 510 } else {
480 511 indexes[3] = 0;
512 }
513}
481 514
482/*---------------------------------------------------------------------------*\ 515/*---------------------------------------------------------------------------*\
483 516
@@ -487,22 +520,22 @@ void newamp1_model_to_indexes(C2CONST *c2const,
487 520
488\*---------------------------------------------------------------------------*/ 521\*---------------------------------------------------------------------------*/
489 522
490void newamp1_interpolate(float interpolated_surface_[], float left_vec[], float right_vec[], int K) 523void newamp1_interpolate(float interpolated_surface_[], float left_vec[],
491{ 524 float right_vec[], int K) {
492 int i, k; 525 int i, k;
493 int M = 4; 526 int M = 4;
494 float c; 527 float c;
495 528
496 /* (linearly) interpolate 25Hz amplitude vectors back to 100Hz */ 529 /* (linearly) interpolate 25Hz amplitude vectors back to 100Hz */
497 530
498 for(i=0,c=1.0; i<M; i++,c-=1.0/M) { 531 for (i = 0, c = 1.0; i < M; i++, c -= 1.0 / M) {
499 for(k=0; k<K; k++) { 532 for (k = 0; k < K; k++) {
500 interpolated_surface_[i*K+k] = left_vec[k]*c + right_vec[k]*(1.0-c); 533 interpolated_surface_[i * K + k] =
501 } 534 left_vec[k] * c + right_vec[k] * (1.0 - c);
502 } 535 }
536 }
503} 537}
504 538
505
506/*---------------------------------------------------------------------------*\ 539/*---------------------------------------------------------------------------*\
507 540
508 FUNCTION....: newamp1_indexes_to_rate_K_vec 541 FUNCTION....: newamp1_indexes_to_rate_K_vec
@@ -514,42 +547,39 @@ void newamp1_interpolate(float interpolated_surface_[], float left_vec[], float
514 547
515\*---------------------------------------------------------------------------*/ 548\*---------------------------------------------------------------------------*/
516 549
517void newamp1_indexes_to_rate_K_vec(float rate_K_vec_[], 550void newamp1_indexes_to_rate_K_vec(float rate_K_vec_[],
518 float rate_K_vec_no_mean_[], 551 float rate_K_vec_no_mean_[],
519 float rate_K_sample_freqs_kHz[], 552 float rate_K_sample_freqs_kHz[], int K,
520 int K, 553 float *mean_, int indexes[],
521 float *mean_,
522 int indexes[],
523 float user_rate_K_vec_no_mean_[], 554 float user_rate_K_vec_no_mean_[],
524 int post_filter_en) 555 int post_filter_en) {
525{ 556 int k;
526 int k; 557 const float *codebook1 = newamp1vq_cb[0].cb;
527 const float *codebook1 = newamp1vq_cb[0].cb; 558 const float *codebook2 = newamp1vq_cb[1].cb;
528 const float *codebook2 = newamp1vq_cb[1].cb; 559 int n1 = indexes[0];
529 int n1 = indexes[0]; 560 int n2 = indexes[1];
530 int n2 = indexes[1]; 561
531 562 if (user_rate_K_vec_no_mean_ == NULL) {
532 if (user_rate_K_vec_no_mean_ == NULL) { 563 /* normal operation */
533 /* normal operation */ 564 for (k = 0; k < K; k++) {
534 for(k=0; k<K; k++) { 565 rate_K_vec_no_mean_[k] = codebook1[K * n1 + k] + codebook2[K * n2 + k];
535 rate_K_vec_no_mean_[k] = codebook1[K*n1+k] + codebook2[K*n2+k];
536 }
537 } else {
538 /* for development we can optionally inject the quantised rate K vector here */
539 for(k=0; k<K; k++)
540 rate_K_vec_no_mean_[k] = user_rate_K_vec_no_mean_[k];
541 } 566 }
542 567 } else {
543 if (post_filter_en) 568 /* for development we can optionally inject the quantised rate K vector here
544 post_filter_newamp1(rate_K_vec_no_mean_, rate_K_sample_freqs_kHz, K, 1.5); 569 */
570 for (k = 0; k < K; k++)
571 rate_K_vec_no_mean_[k] = user_rate_K_vec_no_mean_[k];
572 }
545 573
546 *mean_ = newamp1_energy_cb[0].cb[indexes[2]]; 574 if (post_filter_en)
575 post_filter_newamp1(rate_K_vec_no_mean_, rate_K_sample_freqs_kHz, K, 1.5);
547 576
548 for(k=0; k<K; k++) { 577 *mean_ = newamp1_energy_cb[0].cb[indexes[2]];
549 rate_K_vec_[k] = rate_K_vec_no_mean_[k] + *mean_;
550 }
551}
552 578
579 for (k = 0; k < K; k++) {
580 rate_K_vec_[k] = rate_K_vec_no_mean_[k] + *mean_;
581 }
582}
553 583
554/*---------------------------------------------------------------------------*\ 584/*---------------------------------------------------------------------------*\
555 585
@@ -561,78 +591,66 @@ void newamp1_indexes_to_rate_K_vec(float rate_K_vec_[],
561 591
562\*---------------------------------------------------------------------------*/ 592\*---------------------------------------------------------------------------*/
563 593
564void newamp1_indexes_to_model(C2CONST *c2const, 594void newamp1_indexes_to_model(C2CONST *c2const, MODEL model_[], COMP H[],
565 MODEL model_[],
566 COMP H[],
567 float *interpolated_surface_, 595 float *interpolated_surface_,
568 float prev_rate_K_vec_[], 596 float prev_rate_K_vec_[], float *Wo_left,
569 float *Wo_left, 597 int *voicing_left,
570 int *voicing_left, 598 float rate_K_sample_freqs_kHz[], int K,
571 float rate_K_sample_freqs_kHz[], 599 codec2_fft_cfg fwd_cfg, codec2_fft_cfg inv_cfg,
572 int K, 600 int indexes[], float user_rate_K_vec_no_mean_[],
573 codec2_fft_cfg fwd_cfg, 601 int post_filter_en) {
574 codec2_fft_cfg inv_cfg, 602 float rate_K_vec_[K], rate_K_vec_no_mean_[K], mean_, Wo_right;
575 int indexes[], 603 int voicing_right, k;
576 float user_rate_K_vec_no_mean_[], 604 int M = 4;
577 int post_filter_en) 605
578{ 606 /* extract latest rate K vector */
579 float rate_K_vec_[K], rate_K_vec_no_mean_[K], mean_, Wo_right; 607
580 int voicing_right, k; 608 newamp1_indexes_to_rate_K_vec(rate_K_vec_, rate_K_vec_no_mean_,
581 int M = 4; 609 rate_K_sample_freqs_kHz, K, &mean_, indexes,
582 610 user_rate_K_vec_no_mean_, post_filter_en);
583 /* extract latest rate K vector */ 611
584 612 /* decode latest Wo and voicing */
585 newamp1_indexes_to_rate_K_vec(rate_K_vec_, 613
586 rate_K_vec_no_mean_, 614 if (indexes[3]) {
587 rate_K_sample_freqs_kHz, 615 Wo_right = decode_log_Wo(c2const, indexes[3], 6);
588 K, 616 voicing_right = 1;
589 &mean_, 617 } else {
590 indexes, 618 Wo_right = 2.0 * M_PI / 100.0;
591 user_rate_K_vec_no_mean_, 619 voicing_right = 0;
592 post_filter_en); 620 }
593
594 /* decode latest Wo and voicing */
595
596 if (indexes[3]) {
597 Wo_right = decode_log_Wo(c2const, indexes[3], 6);
598 voicing_right = 1;
599 }
600 else {
601 Wo_right = 2.0*M_PI/100.0;
602 voicing_right = 0;
603 }
604
605 /* interpolate 25Hz rate K vec back to 100Hz */
606 621
607 float *left_vec = prev_rate_K_vec_; 622 /* interpolate 25Hz rate K vec back to 100Hz */
608 float *right_vec = rate_K_vec_;
609 newamp1_interpolate(interpolated_surface_, left_vec, right_vec, K);
610 623
611 /* interpolate 25Hz v and Wo back to 100Hz */ 624 float *left_vec = prev_rate_K_vec_;
625 float *right_vec = rate_K_vec_;
626 newamp1_interpolate(interpolated_surface_, left_vec, right_vec, K);
612 627
613 float aWo_[M]; 628 /* interpolate 25Hz v and Wo back to 100Hz */
614 int avoicing_[M], aL_[M], i;
615 629
616 interp_Wo_v(aWo_, aL_, avoicing_, *Wo_left, Wo_right, *voicing_left, voicing_right); 630 float aWo_[M];
631 int avoicing_[M], aL_[M], i;
617 632
618 /* back to rate L amplitudes, synthesis phase for each frame */ 633 interp_Wo_v(aWo_, aL_, avoicing_, *Wo_left, Wo_right, *voicing_left,
634 voicing_right);
619 635
620 for(i=0; i<M; i++) { 636 /* back to rate L amplitudes, synthesise phase for each frame */
621 model_[i].Wo = aWo_[i];
622 model_[i].L = aL_[i];
623 model_[i].voiced = avoicing_[i];
624 637
625 resample_rate_L(c2const, &model_[i], &interpolated_surface_[K*i], rate_K_sample_freqs_kHz, K); 638 for (i = 0; i < M; i++) {
626 determine_phase(c2const, &H[(MAX_AMP+1)*i], &model_[i], NEWAMP1_PHASE_NFFT, fwd_cfg, inv_cfg); 639 model_[i].Wo = aWo_[i];
627 } 640 model_[i].L = aL_[i];
641 model_[i].voiced = avoicing_[i];
628 642
629 /* update memories for next time */ 643 resample_rate_L(c2const, &model_[i], &interpolated_surface_[K * i],
644 rate_K_sample_freqs_kHz, K);
645 determine_phase(c2const, &H[(MAX_AMP + 1) * i], &model_[i],
646 NEWAMP1_PHASE_NFFT, fwd_cfg, inv_cfg);
647 }
630 648
631 for(k=0; k<K; k++) { 649 /* update memories for next time */
632 prev_rate_K_vec_[k] = rate_K_vec_[k];
633 }
634 *Wo_left = Wo_right;
635 *voicing_left = voicing_right;
636 650
651 for (k = 0; k < K; k++) {
652 prev_rate_K_vec_[k] = rate_K_vec_[k];
653 }
654 *Wo_left = Wo_right;
655 *voicing_left = voicing_right;
637} 656}
638
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-11 20:24:00 +0000