1 /* linearPrediction.cpp - source file for class providing linear prediction capability
2 * written by C. R. Helmrich, last modified in 2021 - see License.htm for legal notices
3 *
4 * The copyright in this software is being made available under the exhale Copyright License
5 * and comes with ABSOLUTELY NO WARRANTY. This software may be subject to other third-
6 * party rights, including patent rights. No such rights are granted under this License.
7 *
8 * Copyright (c) 2018-2021 Christian R. Helmrich, project ecodis. All rights reserved.
9 */
10
11 #include "exhaleLibPch.h"
12 #include "linearPrediction.h"
13
14 // reconstructed 3-bit TNS coefficients
15 static const short tnsQuantCoeff3[9 /*2^3+1*/] = { // = round (2^11 * sin (x * pi / (x < 0 ? 9 : 7)))
16 -2017, -1774, -1316, -700, 0, 889, 1601, 1997, 1997
17 };
18
19 // reconstructed 4-bit TNS coefficients
20 static const short tnsQuantCoeff4[17/*2^4+1*/] = { // = round (2^11 * sin (x * pi / (x < 0 ? 17 : 15)))
21 -2039, -1970, -1833, -1634, -1380, -1078, -740, -376, 0, 426, 833, 1204, 1522, 1774, 1948, 2037, 2037
22 };
23
24 static const short* tnsQuantCoeff[2/*coefRes*/] = {tnsQuantCoeff3, tnsQuantCoeff4};
25
26 // ISO/IEC 14496-3, Sec. 4.6.9.3, 3-bit
27 static const int8_t tnsQuantIndex3[SCHAR_MAX + 1] = { // = round (asin (x / 64) * (x < 0 ? 9 : 7) / pi)
28 -4, -4, -4, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -3, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2,
29 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
30 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
31 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3
32 };
33
34 // ISO/IEC 14496-3, Sec. 4.6.9.3, 4-bit
35 static const int8_t tnsQuantIndex4[SCHAR_MAX + 1] = { // = round (asin (x / 64) * (x < 0 ? 17 : 15) / pi)
36 -8, -7, -7, -7, -6, -6, -6, -6, -6, -5, -5, -5, -5, -5, -5, -5, -4, -4, -4, -4, -4, -4, -4, -4, -4, -3, -3, -3, -3, -3, -3, -3,
37 -3, -3, -3, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0,
38 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,
39 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7
40 };
41
42 static const int8_t* tnsQuantIndex[2/*coefRes*/] = {tnsQuantIndex3, tnsQuantIndex4};
43
44 // static helper functions
quantizeParCorCoeffs(const short * const parCorCoeffs,const uint16_t nCoeffs,const short bitDepth,int8_t * const quantCoeffs,const bool lowRes)45 static int quantizeParCorCoeffs (const short* const parCorCoeffs, const uint16_t nCoeffs, const short bitDepth, int8_t* const quantCoeffs,
46 const bool lowRes)
47 {
48 const short bitShift = bitDepth - 7;
49 const unsigned tabIdx = (lowRes ? 0 : 1);
50 const int8_t tabOffset = 4 << tabIdx;
51 const short* coeffTab = tnsQuantCoeff[tabIdx];
52 const int8_t* indexTab = tnsQuantIndex[tabIdx];
53 int dist0, dist1, distTotal = 0;
54
55 for (uint16_t s = 0; s < nCoeffs; s++)
56 {
57 const short coeff = (bitShift < 0 ? parCorCoeffs[s] << -bitShift : parCorCoeffs[s] >> bitShift);
58 const int8_t coeff1 = indexTab[coeff + 1 + (SCHAR_MAX >> 1)];
59 const int8_t coeff0 = (coeff1 <= -tabOffset ? coeff1 : coeff1 - 1);
60
61 dist0 = (int) coeffTab[coeff0 + tabOffset] - parCorCoeffs[s];
62 dist0 *= dist0;
63 dist1 = (int) coeffTab[coeff1 + tabOffset] - parCorCoeffs[s];
64 dist1 *= dist1;
65 if (dist0 < dist1) // use down-rounded coeff
66 {
67 quantCoeffs[s] = coeff0;
68 distTotal += dist0;
69 }
70 else // dist0 >= dist1, use up-rounded coeff
71 {
72 quantCoeffs[s] = ((dist0 == dist1) && (abs (coeff0) < abs (coeff1)) ? coeff0 : coeff1);
73 distTotal += dist1;
74 }
75 } // for s
76
77 return distTotal; // total quantization error
78 }
79
80 // constructor
LinearPredictor()81 LinearPredictor::LinearPredictor ()
82 {
83 // initialize member variables
84 memset (m_tempBuf, 0, 2 * MAX_PREDICTION_ORDER * sizeof (int64_t));
85 }
86
87 // public functions
calcParCorCoeffs(const int32_t * const anaSignal,const uint16_t nAnaSamples,const uint16_t nCoeffs,short * const parCorCoeffs)88 uint32_t LinearPredictor::calcParCorCoeffs (const int32_t* const anaSignal, const uint16_t nAnaSamples, const uint16_t nCoeffs,
89 short* const parCorCoeffs) // returns 256 - 256 / prediction gain per filter order, or 0
90 {
91 // int64 version of algorithm in Figure 1 of J. LeRoux and C. Gueguen, "A Fixed Point Computation
92 // of Partial Correlation Coefficients," IEEE Trans. ASSP, vol. 27, no. 3, pp. 257-259, June 1977.
93 int64_t* const acf = m_tempBuf; // correlation
94 uint32_t pg[MAX_PREDICTION_ORDER] = {0, 0, 0, 0};
95 int64_t pgDen, pgOff; // for prediction gains
96 short s = nAnaSamples - 1;
97
98 if ((anaSignal == nullptr) || (parCorCoeffs == nullptr) || (nCoeffs == 0) || (nCoeffs > MAX_PREDICTION_ORDER) || (nAnaSamples <= nCoeffs))
99 {
100 return 0;
101 }
102
103 if (nCoeffs >= 2) // high predictor order, 2-4
104 {
105 int64_t* const EN = &m_tempBuf[0];
106 int64_t* const EP = &m_tempBuf[MAX_PREDICTION_ORDER];
107 int64_t sampleHO = anaSignal[s];
108
109 // calculate autocorrelation function values
110 acf[0] = sampleHO * sampleHO;
111 sampleHO = anaSignal[--s];
112 acf[0] += sampleHO * sampleHO;
113 acf[1] = sampleHO * anaSignal[s + 1];
114 sampleHO = anaSignal[--s];
115 acf[0] += sampleHO * sampleHO;
116 acf[1] += sampleHO * anaSignal[s + 1];
117 acf[2] = sampleHO * anaSignal[s + 2];
118
119 if (nCoeffs == 2)
120 {
121 for (s--; s >= 0; s--)
122 {
123 const int64_t sample = anaSignal[s];
124
125 acf[0] += sample * sample;
126 acf[1] += sample * anaSignal[s + 1];
127 acf[2] += sample * anaSignal[s + 2];
128 }
129 }
130 else // order 3-4
131 {
132 sampleHO = anaSignal[--s];
133 acf[0] += sampleHO * sampleHO;
134 acf[1] += sampleHO * anaSignal[s + 1];
135 acf[2] += sampleHO * anaSignal[s + 2];
136 acf[3] = sampleHO * anaSignal[s + 3];
137
138 if (nCoeffs == 3)
139 {
140 for (s--; s >= 0; s--)
141 {
142 const int64_t sample = anaSignal[s];
143
144 acf[0] += sample * sample;
145 acf[1] += sample * anaSignal[s + 1];
146 acf[2] += sample * anaSignal[s + 2];
147 acf[3] += sample * anaSignal[s + 3];
148 }
149 }
150 else // order 4
151 {
152 acf[4] = 0;
153 for (s--; s >= 0; s--)
154 {
155 const int64_t sample = anaSignal[s];
156
157 acf[0] += sample * sample;
158 acf[1] += sample * anaSignal[s + 1];
159 acf[2] += sample * anaSignal[s + 2];
160 acf[3] += sample * anaSignal[s + 3];
161 acf[4] += sample * anaSignal[s + 4];
162 }
163 }
164 }
165
166 // reduce correlation value range to <32 bit
167 acf[0] = (acf[0] - INT_MIN/*eps*/) >> 31;
168
169 EP[3] = (acf[4] - (INT_MIN >> 1)) >> 31;
170 EN[3] = (acf[3] - (INT_MIN >> 1)) >> 31;
171 EP[2] = EN[3];
172 EN[2] = (acf[2] - (INT_MIN >> 1)) >> 31;
173 EP[1] = EN[2];
174 EN[1] = (acf[1] - (INT_MIN >> 1)) >> 31;
175 EP[0] = EN[1]; // finish EP buffer creation
176 EN[0] = acf[0]; // finish EN buffer creation
177
178 pgOff = acf[0];
179 pgDen = pgOff << 1;
180
181 for (s = 0; s < (short) nCoeffs; s++)
182 {
183 uint16_t p;
184 // ParCor coefficient Ks & prediction gain
185 int64_t Ks = (EP[s] * -512) / (EN[0] + LP_EPS);
186
187 Ks = CLIP_PM (Ks, 511); // enforce |Ks|<1
188 parCorCoeffs[s] = (short) Ks;
189
190 for (p = s; p < nCoeffs; p++)
191 {
192 sampleHO = EN[p - s]; // use register?
193 EN[p - s] = (EN[p - s] << 9) + EP[p] * Ks;
194 EP[p] = (EP[p] << 9) + sampleHO * Ks;
195 }
196 if (s > 0 && EN[0] < 0) // EN wrap-around
197 {
198 pg[s] = pg[s - 1]; // "invent" some gain
199 }
200 else
201 {
202 sampleHO = (1 << (s * 9)) >> 1;// offset
203 pg[s] = (pgOff <= LP_EPS ? 0 : (1 << 8) - uint32_t ((((EN[0] + sampleHO) >> (s * 9)) + pgOff) / pgDen));
204 }
205 } // for s
206 }
207 else // nCoeffs == 1, minimum predictor order
208 {
209 int64_t sampleMO = anaSignal[s];
210
211 acf[0] = sampleMO * sampleMO;
212 acf[1] = 0;
213 for (s--; s >= 0; s--)
214 {
215 const int64_t sample = anaSignal[s];
216
217 acf[0] += sample * sample;
218 acf[1] += sample * anaSignal[s + 1];
219 }
220
221 // reduce correlation value range to <32 bit
222 acf[0] = (acf[0] - INT_MIN/*eps*/) >> 31;
223 acf[1] = (acf[1] - (INT_MIN >> 1)) >> 31;
224
225 pgOff = acf[0];
226 pgDen = pgOff << 1;
227
228 // 1st-order coefficient K & prediction gain
229 sampleMO = (acf[1] * -512) / (acf[0] + LP_EPS);
230
231 sampleMO = CLIP_PM (sampleMO, 511); // |K|<1
232 parCorCoeffs[0] = (short) sampleMO;
233
234 sampleMO = (acf[0] << 9) + acf[1] * sampleMO;
235 pg[0] = (pgOff <= LP_EPS ? 0 : (1 << 8) - uint32_t ((sampleMO + pgOff) / pgDen));
236 }
237
238 return (CLIP_UCHAR (pg[3]) << 24) | (CLIP_UCHAR (pg[2]) << 16) | (CLIP_UCHAR (pg[1]) << 8) | CLIP_UCHAR (pg[0]);
239 }
240
calcOptTnsCoeffs(short * const parCorCoeffs,int8_t * const quantCoeffs,bool * const lowCoeffRes,const uint16_t maxOrder,const uint8_t predGain,const uint8_t tonality,const uint16_t parCorCoeffBitDepth)241 uint8_t LinearPredictor::calcOptTnsCoeffs (short* const parCorCoeffs, int8_t* const quantCoeffs, bool* const lowCoeffRes,
242 const uint16_t maxOrder, const uint8_t predGain, const uint8_t tonality /*= 0*/,
243 const uint16_t parCorCoeffBitDepth /*= 10*/) // returns optimized filter order for TNS
244 {
245 const short bitShift = LP_DEPTH - (short) parCorCoeffBitDepth;
246 short shortBuf[MAX_PREDICTION_ORDER];
247 uint16_t s, order = __min (maxOrder, MAX_PREDICTION_ORDER);
248 int d, i;
249
250 if ((parCorCoeffs == nullptr) || (quantCoeffs == nullptr) || (maxOrder == 0) || (maxOrder > MAX_PREDICTION_ORDER) || (parCorCoeffBitDepth < 2) || (bitShift < 0))
251 {
252 if (quantCoeffs) memset (quantCoeffs, 0, order * sizeof (int8_t));
253
254 return 0; // invalid input arguments error
255 }
256
257 // determine direct-form filter damping factor
258 parCorCoeffs[0] <<= bitShift;
259 i = abs (parCorCoeffs[0]);
260 for (s = 1; s < order; s++)
261 {
262 parCorCoeffs[s] <<= bitShift; // scale coeff
263 i = __max (i, abs (parCorCoeffs[s]));
264 }
265 for (/*s*/; s < MAX_PREDICTION_ORDER; s++)
266 {
267 parCorCoeffs[s] = 0; // similarParCorCoeffs
268 }
269
270 if (predGain < 41 + (tonality >> 3)) // 1.5 dB
271 {
272 memset (quantCoeffs, 0, order * sizeof (int8_t));
273
274 return 0; // LPC prediction gain is too low
275 }
276
277 d = (3 << (LP_DEPTH - 1)) >> 2;
278
279 if (i > d) // apply direct-form filter damping
280 {
281 i = ((i >> 1) + d * (1 << LP_SHIFT)) / i;
282 d = i;
283 if (parCorToLpCoeffs (parCorCoeffs, order, shortBuf, LP_DEPTH) > 0)
284 {
285 return 0; // coefficient conversion error
286 }
287 for (s = 0; s < order; s++)
288 {
289 shortBuf[s] = short ((LP_OFFSET + d * shortBuf[s]) >> LP_SHIFT);
290 d = (LP_OFFSET + d * i) >> LP_SHIFT;
291 }
292 // verify order and go back to ParCor domain
293 if ((s > 0) && (shortBuf[s - 1] == 0))
294 {
295 order--;
296 }
297 if (lpToParCorCoeffs (shortBuf, order, parCorCoeffs, LP_DEPTH) > 0)
298 {
299 return 0; // coefficient conversion error
300 }
301 }
302
303 // high-res quantizer, obtain coeff distortion
304 d = quantizeParCorCoeffs (parCorCoeffs, order, LP_DEPTH, quantCoeffs, false);
305
306 if ((lowCoeffRes != nullptr) && (quantizeParCorCoeffs (parCorCoeffs, order, LP_DEPTH, (int8_t* const) m_tempBuf, true) < d))
307 {
308 // low-res quantizer yields lower distortion
309 *lowCoeffRes = true;
310 memcpy (quantCoeffs, m_tempBuf, order * sizeof (int8_t));
311 }
312
313 for (; order > 0; order--) // return opt order
314 {
315 if (quantCoeffs[order - 1] != 0) return (uint8_t) order;
316 }
317 return 0;
318 }
319
lpToParCorCoeffs(short * const lpCoeffs,const uint16_t nCoeffs,short * const parCorCoeffs,const uint16_t parCorCoeffBitDepth)320 unsigned LinearPredictor::lpToParCorCoeffs (/*mod!*/short* const lpCoeffs, const uint16_t nCoeffs, short* const parCorCoeffs,
321 const uint16_t parCorCoeffBitDepth /*= 10*/)
322 {
323 int* const intBuf = (int* const) m_tempBuf;
324 const int shift = parCorCoeffBitDepth - 1;
325 const int offset = 1 << (shift - 1);
326
327 if ((lpCoeffs == nullptr) || (parCorCoeffs == nullptr) || (nCoeffs == 0) || (nCoeffs > MAX_PREDICTION_ORDER) || (parCorCoeffBitDepth < 2))
328 {
329 return 1; // error
330 }
331
332 for (uint16_t p, s = nCoeffs - 1; s > 0; s--)
333 {
334 const int i = (parCorCoeffs[s] = lpCoeffs[s]);
335 const int d = (1 << shift) - ((offset + i * i) >> shift);
336 const int o = d >> 1;
337
338 if (d <= 0) return s; // invalid coefficient
339
340 for (p = 0; p < s; p++)
341 {
342 intBuf[p] = lpCoeffs[s - 1 - p];
343 }
344 for (p = 0; p < s; p++)
345 {
346 lpCoeffs[p] = short ((o + ((int) lpCoeffs[p] << shift) - intBuf[p] * i) / d);
347 }
348 }
349 parCorCoeffs[0] = lpCoeffs[0];
350
351 return 0; // no error
352 }
353
parCorToLpCoeffs(const short * const parCorCoeffs,const uint16_t nCoeffs,short * const lpCoeffs,const uint16_t parCorCoeffBitDepth)354 unsigned LinearPredictor::parCorToLpCoeffs (const short* const parCorCoeffs, const uint16_t nCoeffs, short* const lpCoeffs,
355 const uint16_t parCorCoeffBitDepth /*= 10*/)
356 {
357 int* const intBuf = (int* const) m_tempBuf;
358 const int shift = parCorCoeffBitDepth - 1;
359 const int offset = 1 << (shift - 1);
360
361 if ((parCorCoeffs == nullptr) || (lpCoeffs == nullptr) || (nCoeffs == 0) || (nCoeffs > MAX_PREDICTION_ORDER) || (parCorCoeffBitDepth < 2))
362 {
363 return 1; // error
364 }
365
366 lpCoeffs[0] = parCorCoeffs[0];
367 for (uint16_t p, s = 1; s < nCoeffs; s++)
368 {
369 const int i = (lpCoeffs[s] = parCorCoeffs[s]);
370
371 if (abs (i) > (1 << shift)) return s; // > 1
372
373 for (p = 0; p < s; p++)
374 {
375 intBuf[p] = lpCoeffs[s - 1 - p];
376 }
377 for (p = 0; p < s; p++)
378 {
379 lpCoeffs[p] += short ((offset + intBuf[p] * i) >> shift);
380 }
381 }
382 return 0; // no error
383 }
384
quantTnsToLpCoeffs(const int8_t * const quantTnsCoeffs,const uint16_t nCoeffs,const bool lowCoeffRes,short * const parCorCoeffs,short * const lpCoeffs)385 unsigned LinearPredictor::quantTnsToLpCoeffs (const int8_t* const quantTnsCoeffs, const uint16_t nCoeffs, const bool lowCoeffRes,
386 short* const parCorCoeffs, short* const lpCoeffs)
387 {
388 const unsigned tabIdx = (lowCoeffRes ? 0 : 1);
389 const int8_t tabOffset = 4 << tabIdx;
390 const short* coeffTab = tnsQuantCoeff[tabIdx];
391
392 if ((quantTnsCoeffs == nullptr) || (parCorCoeffs == nullptr) || (lpCoeffs == nullptr) || (nCoeffs == 0) || (nCoeffs > MAX_PREDICTION_ORDER))
393 {
394 return 1; // error
395 }
396
397 for (uint16_t s = 0; s < nCoeffs; s++)
398 {
399 parCorCoeffs[s] = coeffTab[CLIP_PM (quantTnsCoeffs[s], tabOffset) + tabOffset];
400 }
401
402 return parCorToLpCoeffs (parCorCoeffs, nCoeffs, lpCoeffs, LP_DEPTH);
403 }
404
similarParCorCoeffs(const short * const parCorCoeffs1,const short * const parCorCoeffs2,const uint16_t nCoeffs,const uint16_t parCorCoeffBitDepth)405 bool LinearPredictor::similarParCorCoeffs (const short* const parCorCoeffs1, const short* const parCorCoeffs2, const uint16_t nCoeffs,
406 const uint16_t parCorCoeffBitDepth /*= 10*/)
407 {
408 unsigned sumAbsDiff = 0;
409
410 if ((parCorCoeffs1 == nullptr) || (parCorCoeffs2 == nullptr) || (nCoeffs == 0) || (nCoeffs > MAX_PREDICTION_ORDER) || (parCorCoeffBitDepth < 2))
411 {
412 return false; // error
413 }
414
415 for (uint16_t s = 0; s < nCoeffs; s++)
416 {
417 sumAbsDiff += abs (parCorCoeffs1[s] - parCorCoeffs2[s]);
418 }
419
420 return (sumAbsDiff + 12u * nCoeffs < ((4u * nCoeffs) << (parCorCoeffBitDepth >> 1)));
421 }
422