FreeCalypso > hg > gsm-codec-lib
view libgsmefr/az_lsp.c @ 242:f081a6850fb5
libgsmfrp: new refined implementation
The previous implementation exhibited the following defects,
which are now fixed:
1) The last received valid SID was cached forever for the purpose of
handling future invalid SIDs - we could have received some valid
SID ages ago, then lots of speech or NO_DATA, and if we then get
an invalid SID, we would resurrect the last valid SID from ancient
history - a bad design. In our new design, we handle invalid SID
based on the current state, much like BFI.
2) GSM 06.11 spec says clearly that after the second lost SID
(received BFI=1 && TAF=1 in CN state) we need to gradually decrease
the output level, rather than jump directly to emitting silence
frames - we previously failed to implement such logic.
3) Per GSM 06.12 section 5.2, Xmaxc should be the same in all 4 subframes
in a SID frame. What should we do if we receive an otherwise valid
SID frame with different Xmaxc? Our previous approach would
replicate this Xmaxc oddity in every subsequent generated CN frame,
which is rather bad. In our new design, the very first CN frame
(which can be seen as a transformation of the SID frame itself)
retains the original 4 distinct Xmaxc, but all subsequent CN frames
are based on the Xmaxc from the last subframe of the most recent SID.
author | Mychaela Falconia <falcon@freecalypso.org> |
---|---|
date | Tue, 09 May 2023 05:16:31 +0000 |
parents | a4d1615e2aa4 |
children |
line wrap: on
line source
/*********************************************************************** * * FUNCTION: Az_lsp * * PURPOSE: Compute the LSPs from the LP coefficients (order=10) * * DESCRIPTION: * - The sum and difference filters are computed and divided by * 1+z^{-1} and 1-z^{-1}, respectively. * * f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5 * f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5 * * - The roots of F1(z) and F2(z) are found using Chebyshev polynomial * evaluation. The polynomials are evaluated at 60 points regularly * spaced in the frequency domain. The sign change interval is * subdivided 4 times to better track the root. * The LSPs are found in the cosine domain [1,-1]. * * - If less than 10 roots are found, the LSPs from the past frame are * used. * ***********************************************************************/ #include "gsm_efr.h" #include "typedef.h" #include "namespace.h" #include "basic_op.h" #include "oper_32b.h" #include "no_count.h" #include "cnst.h" #include "sig_proc.h" #include "grid.tab" /* M = LPC order, NC = M/2 */ #define NC M/2 /* local function */ static Word16 Chebps (Word16 x, Word16 f[], Word16 n); void Az_lsp ( Word16 a[], /* (i) : predictor coefficients */ Word16 lsp[], /* (o) : line spectral pairs */ Word16 old_lsp[] /* (i) : old lsp[] (in case not found 10 roots) */ ) { Word16 i, j, nf, ip; Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint; Word16 x, y, sign, exp; Word16 *coef; Word16 f1[M / 2 + 1], f2[M / 2 + 1]; Word32 t0; /*-------------------------------------------------------------* * find the sum and diff. pol. F1(z) and F2(z) * * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) * * * * f1[0] = 1.0; * * f2[0] = 1.0; * * * * for (i = 0; i< NC; i++) * * { * * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; * * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; * * } * *-------------------------------------------------------------*/ f1[0] = 1024; move16 (); /* f1[0] = 1.0 */ f2[0] = 1024; move16 (); /* f2[0] = 1.0 */ for (i = 0; i < NC; i++) { t0 = L_mult (a[i + 1], 8192); /* x = (a[i+1] + a[M-i]) >> 2 */ t0 = L_mac (t0, a[M - i], 8192); x = extract_h (t0); /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */ f1[i + 1] = sub (x, f1[i]);move16 (); t0 = L_mult (a[i + 1], 8192); /* x = (a[i+1] - a[M-i]) >> 2 */ t0 = L_msu (t0, a[M - i], 8192); x = extract_h (t0); /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */ f2[i + 1] = add (x, f2[i]);move16 (); } /*-------------------------------------------------------------* * find the LSPs using the Chebychev pol. evaluation * *-------------------------------------------------------------*/ nf = 0; move16 (); /* number of found frequencies */ ip = 0; move16 (); /* indicator for f1 or f2 */ coef = f1; move16 (); xlow = grid[0]; move16 (); ylow = Chebps (xlow, coef, NC);move16 (); j = 0; test (); test (); while ( (nf < M) && (j < grid_points) ) /* while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0)) */ { j++; xhigh = xlow; move16 (); yhigh = ylow; move16 (); xlow = grid[j]; move16 (); ylow = Chebps (xlow, coef, NC); move16 (); test (); if (L_mult (ylow, yhigh) <= (Word32) 0L) { /* divide 4 times the interval */ for (i = 0; i < 4; i++) { /* xmid = (xlow + xhigh)/2 */ xmid = add (shr (xlow, 1), shr (xhigh, 1)); ymid = Chebps (xmid, coef, NC); move16 (); test (); if (L_mult (ylow, ymid) <= (Word32) 0L) { yhigh = ymid; move16 (); xhigh = xmid; move16 (); } else { ylow = ymid; move16 (); xlow = xmid; move16 (); } } /*-------------------------------------------------------------* * Linear interpolation * * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); * *-------------------------------------------------------------*/ x = sub (xhigh, xlow); y = sub (yhigh, ylow); test (); if (y == 0) { xint = xlow; move16 (); } else { sign = y; move16 (); y = abs_s (y); exp = norm_s (y); y = shl (y, exp); y = div_s ((Word16) 16383, y); t0 = L_mult (x, y); t0 = L_shr (t0, sub (20, exp)); y = extract_l (t0); /* y= (xhigh-xlow)/(yhigh-ylow) */ test (); if (sign < 0) y = negate (y); t0 = L_mult (ylow, y); t0 = L_shr (t0, 11); xint = sub (xlow, extract_l (t0)); /* xint = xlow - ylow*y */ } lsp[nf] = xint; move16 (); xlow = xint; move16 (); nf++; test (); if (ip == 0) { ip = 1; move16 (); coef = f2; move16 (); } else { ip = 0; move16 (); coef = f1; move16 (); } ylow = Chebps (xlow, coef, NC); move16 (); } test (); test (); } /* Check if M roots found */ test (); if (nf < M) { for (i = 0; i < M; i++) { lsp[i] = old_lsp[i]; move16 (); } } return; } /************************************************************************ * * FUNCTION: Chebps * * PURPOSE: Evaluates the Chebyshev polynomial series * * DESCRIPTION: * - The polynomial order is n = m/2 = 5 * - The polynomial F(z) (F1(z) or F2(z)) is given by * F(w) = 2 exp(-j5w) C(x) * where * C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2 * and T_m(x) = cos(mw) is the mth order Chebyshev polynomial ( x=cos(w) ) * - The function returns the value of C(x) for the input x. * ***********************************************************************/ static Word16 Chebps (Word16 x, Word16 f[], Word16 n) { Word16 i, cheb; Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l; Word32 t0; b2_h = 256; move16 (); /* b2 = 1.0 */ b2_l = 0; move16 (); t0 = L_mult (x, 512); /* 2*x */ t0 = L_mac (t0, f[1], 8192); /* + f[1] */ L_Extract (t0, &b1_h, &b1_l); /* b1 = 2*x + f[1] */ for (i = 2; i < n; i++) { t0 = Mpy_32_16 (b1_h, b1_l, x); /* t0 = 2.0*x*b1 */ t0 = L_shl (t0, 1); t0 = L_mac (t0, b2_h, (Word16) 0x8000); /* t0 = 2.0*x*b1 - b2 */ t0 = L_msu (t0, b2_l, 1); t0 = L_mac (t0, f[i], 8192); /* t0 = 2.0*x*b1 - b2 + f[i] */ L_Extract (t0, &b0_h, &b0_l); /* b0 = 2.0*x*b1 - b2 + f[i]*/ b2_l = b1_l; move16 (); /* b2 = b1; */ b2_h = b1_h; move16 (); b1_l = b0_l; move16 (); /* b1 = b0; */ b1_h = b0_h; move16 (); } t0 = Mpy_32_16 (b1_h, b1_l, x); /* t0 = x*b1; */ t0 = L_mac (t0, b2_h, (Word16) 0x8000); /* t0 = x*b1 - b2 */ t0 = L_msu (t0, b2_l, 1); t0 = L_mac (t0, f[i], 4096); /* t0 = x*b1 - b2 + f[i]/2 */ t0 = L_shl (t0, 6); cheb = extract_h (t0); return (cheb); }