From 30eb2f854ce70315a9d483b327069069ef334291 Mon Sep 17 00:00:00 2001 From: Thomas Vilte Date: Mon, 20 Jul 2026 18:34:14 -0300 Subject: [PATCH] Add SILK LPC-to-NLSF conversion (a2nlsf) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: François Allais --- internal/silk/a2nlsf.go | 196 +++++++++++++++++++++++++++++++++++ internal/silk/a2nlsf_test.go | 151 +++++++++++++++++++++++++++ 2 files changed, 347 insertions(+) create mode 100644 internal/silk/a2nlsf.go create mode 100644 internal/silk/a2nlsf_test.go diff --git a/internal/silk/a2nlsf.go b/internal/silk/a2nlsf.go new file mode 100644 index 0000000..9512c59 --- /dev/null +++ b/internal/silk/a2nlsf.go @@ -0,0 +1,196 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// LPC-to-NLSF conversion (silk_A2NLSF, A2NLSF.c). The whitening filter's even +// and odd polynomials are root-found on a cosine grid to recover the NLSFs. + +const ( + binDivStepsA2NLSF = 3 + maxIterA2NLSF = 16 + lsfCosTabSz = 128 +) + +// silkLSFCosTabFIXQ12 is the Q12 cosine table used for the root search. +// +//nolint:gochecknoglobals +var silkLSFCosTabFIXQ12 = [lsfCosTabSz + 1]int16{ + 8192, 8190, 8182, 8170, 8152, 8130, 8104, 8072, 8034, 7994, 7946, 7896, + 7840, 7778, 7714, 7644, 7568, 7490, 7406, 7318, 7226, 7128, 7026, 6922, + 6812, 6698, 6580, 6458, 6332, 6204, 6070, 5934, 5792, 5648, 5502, 5352, + 5198, 5040, 4880, 4718, 4552, 4382, 4212, 4038, 3862, 3684, 3502, 3320, + 3136, 2948, 2760, 2570, 2378, 2186, 1990, 1794, 1598, 1400, 1202, 1002, + 802, 602, 402, 202, 0, -202, -402, -602, -802, -1002, -1202, -1400, + -1598, -1794, -1990, -2186, -2378, -2570, -2760, -2948, -3136, -3320, + -3502, -3684, -3862, -4038, -4212, -4382, -4552, -4718, -4880, -5040, + -5198, -5352, -5502, -5648, -5792, -5934, -6070, -6204, -6332, -6458, + -6580, -6698, -6812, -6922, -7026, -7128, -7226, -7318, -7406, -7490, + -7568, -7644, -7714, -7778, -7840, -7896, -7946, -7994, -8034, -8072, + -8104, -8130, -8152, -8170, -8182, -8190, -8192, +} + +// smlaww returns a + ((b * c) >> 16). +func smlaww(a, b, c int32) int32 { + return a + smulww(b, c) +} + +// a2nlsfTransPoly transforms a polynomial from cos(n*f) to cos(f)^n. +func a2nlsfTransPoly(p []int32, dd int) { + for k := 2; k <= dd; k++ { + for n := dd; n > k; n-- { + p[n-2] -= p[n] + } + p[k-2] -= p[k] << 1 + } +} + +// a2nlsfEvalPoly evaluates a Q16 polynomial at the Q12 point x. +func a2nlsfEvalPoly(p []int32, x int32, dd int) int32 { + y32 := p[dd] + xQ16 := x << 4 + for n := dd - 1; n >= 0; n-- { + y32 = smlaww(p[n], y32, xQ16) + } + + return y32 +} + +// a2nlsfInit builds the even (P) and odd (Q) polynomials from the LPC coefs. +func a2nlsfInit(aQ16, p, q []int32, dd int) { + p[dd] = 1 << 16 + q[dd] = 1 << 16 + for k := range dd { + p[k] = -aQ16[dd-k-1] - aQ16[dd+k] + q[k] = -aQ16[dd-k-1] + aQ16[dd+k] + } + for k := dd; k > 0; k-- { + p[k-1] -= p[k] + q[k-1] += q[k] + } + a2nlsfTransPoly(p, dd) + a2nlsfTransPoly(q, dd) +} + +// bwexpander32 applies a Q16 chirp to int32 AR coefficients. +func bwexpander32(ar []int32, d int, chirpQ16 int32) { + chirpMinusOneQ16 := chirpQ16 - 65536 + for i := range d - 1 { + ar[i] = smulww(chirpQ16, ar[i]) + chirpQ16 += rshiftRound32(chirpQ16*chirpMinusOneQ16, 16) + } + ar[d-1] = smulww(chirpQ16, ar[d-1]) +} + +// a2nlsf converts monic whitening filter coefficients (Q16, modified in place) +// to Q15 NLSFs, bandwidth-expanding until the roots converge. +// +//nolint:gocognit,gocyclo,cyclop,maintidx,varnamelen // faithful port of silk_A2NLSF. +func a2nlsf(nlsf []int16, aQ16 []int32, d int) { + dd := d >> 1 + pPoly := make([]int32, dd+1) + qPoly := make([]int32, dd+1) + a2nlsfInit(aQ16, pPoly, qPoly, dd) + polys := [2][]int32{pPoly, qPoly} + + p := pPoly //nolint:varnamelen // p tracks the active P/Q polynomial, as in the C reference. + xlo := int32(silkLSFCosTabFIXQ12[0]) + ylo := a2nlsfEvalPoly(p, xlo, dd) + + var rootIx int + if ylo < 0 { + nlsf[0] = 0 + p = qPoly + ylo = a2nlsfEvalPoly(p, xlo, dd) + rootIx = 1 + } + + k := 1 //nolint:varnamelen // k indexes the cosine grid, as in the C reference. + iterations := 0 + thr := int32(0) + for { + xhi := int32(silkLSFCosTabFIXQ12[k]) + yhi := a2nlsfEvalPoly(p, xhi, dd) + + if (ylo <= 0 && yhi >= thr) || (ylo >= 0 && yhi <= -thr) { //nolint:nestif // faithful port of silk_A2NLSF. + if yhi == 0 { + thr = 1 + } else { + thr = 0 + } + ffrac := int32(-256) + for m := range binDivStepsA2NLSF { + xmid := rshiftRound32(xlo+xhi, 1) + ymid := a2nlsfEvalPoly(p, xmid, dd) + if (ylo <= 0 && ymid >= 0) || (ylo >= 0 && ymid <= 0) { + xhi = xmid + yhi = ymid + } else { + xlo = xmid + ylo = ymid + ffrac += 128 >> m + } + } + if absInt32(ylo) < 65536 { + den := ylo - yhi + nom := (ylo << (8 - binDivStepsA2NLSF)) + (den >> 1) + if den != 0 { + ffrac += nom / den + } + } else { + ffrac += ylo / ((ylo - yhi) >> (8 - binDivStepsA2NLSF)) + } + nlsf[rootIx] = int16(min((int32(k)<<8)+ffrac, math.MaxInt16)) + + rootIx++ + if rootIx >= d { + return + } + p = polys[rootIx&1] + xlo = int32(silkLSFCosTabFIXQ12[k-1]) + ylo = int32(1-(rootIx&2)) << 12 + } else { + k++ + xlo = xhi + ylo = yhi + thr = 0 + + if k > lsfCosTabSz { + iterations++ + if iterations > maxIterA2NLSF { + nlsf[0] = int16((1 << 15) / int32(d+1)) //nolint:gosec // G115: quotient fits int16. + for k = 1; k < d; k++ { + nlsf[k] = nlsf[k-1] + nlsf[0] + } + + return + } + bwexpander32(aQ16, d, 65536-(1< +// SPDX-License-Identifier: MIT + +package silk + +import ( + "fmt" + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// genNLSF builds a deterministic, strictly increasing NLSF vector in Q15. +func genNLSF(order int, seed uint32) []int16 { + nlsf := make([]int16, order) + state := seed + cur := int32(400) + for k := range nlsf { + state = 1664525*state + 1013904223 + cur += int32(state>>20)%2000 + 700 //nolint:gosec // G115: bounded arithmetic. + if cur > 32200 { + cur = 32200 + } + nlsf[k] = int16(cur) + } + + return nlsf +} + +// TestA2NLSFRoundTrip converts a valid NLSF vector to LPC via the decoder and +// back with A2NLSF, expecting to recover the input within the fixed-point +// tolerance. +func TestA2NLSFRoundTrip(t *testing.T) { + cases := []struct { + bandwidth Bandwidth + order int + }{ + {BandwidthNarrowband, 10}, + {BandwidthWideband, 16}, + } + + for _, tc := range cases { + for seed := range 8 { + t.Run(fmt.Sprintf("bw%d_seed%d", tc.bandwidth, seed), func(t *testing.T) { + nlsf := genNLSF(tc.order, uint32(seed*131+tc.order)) //nolint:gosec // G115: non-negative test seed. + + dec := NewDecoder() + dec.normalizeLSFStabilization(nlsf, tc.order, tc.bandwidth) + a32Q17 := dec.convertNormalizedLSFsToLPCCoefficients(nlsf, tc.bandwidth) + aFloat := make([]float32, tc.order) + for i := range aFloat { + aFloat[i] = float32(a32Q17[i]) / (1 << 17) + } + + recovered := make([]int16, tc.order) + a2nlsfFLP(recovered, aFloat, tc.order) + + // Round-trip through two different fixed-point approximations + // (decoder NLSF->LPC and A2NLSF LPC->NLSF); band-edge roots are + // the least precise. ~1% of full scale is expected. + for k := range nlsf { + assert.InDeltaf(t, nlsf[k], recovered[k], 300, "coefficient %d", k) + } + // Result must be a valid increasing NLSF vector. + for k := 1; k < tc.order; k++ { + require.Greaterf(t, recovered[k], recovered[k-1], "not increasing at %d", k) + } + }) + } + } +} + +// TestA2NLSFInitialRootBelowZero exercises the branch where the even +// polynomial's value at cos(0) is already negative on the very first +// evaluation, forcing nlsf[0] to 0 and switching to the odd polynomial to +// search for the first root. The branch is only observable indirectly (via +// nlsf[0]==0, since a2nlsf keeps no other externally visible state) — these +// aQ16 vectors were found by brute-force search over random coefficients, +// not derived analytically. +// +// Unlike TestA2NLSFRoundTrip, these aQ16 values are NOT derived from a valid +// whitening filter (no real encoder would produce them), so the classic +// interlacing-root guarantee that makes NLSFs monotonic doesn't apply here — +// only boundedness does. +func TestA2NLSFInitialRootBelowZero(t *testing.T) { + cases := []struct { + name string + order int + aQ16 []int32 + }{ + {"order10", 10, []int32{479952, -135678, 36600, -1536, 84504, -11100, 118244, 21727, -7612, -30372}}, + {"order16", 16, []int32{82040, -77114, 55740, -3567, 27247, -762, 3861, 1517, -119, 1306, 213, 3, -42, 80, 24, 66}}, + } + for _, tc := range cases { + t.Run(tc.name, func(t *testing.T) { + aQ16 := append([]int32(nil), tc.aQ16...) + nlsf := make([]int16, tc.order) + a2nlsf(nlsf, aQ16, tc.order) + + assert.Equal(t, int16(0), nlsf[0], "expected the initial-root-negative branch to fire") + for k := 1; k < tc.order; k++ { + require.GreaterOrEqualf(t, nlsf[k], int16(0), "coefficient %d negative", k) + } + }) + } +} + +// TestA2NLSFDegenerate exercises the bandwidth-expansion retry path with +// pathological LPC coefficients. Zeros and huge values force the algorithm +// through multiple expansion passes; the output must still be a valid +// increasing NLSF vector. +func TestA2NLSFDegenerate(t *testing.T) { + cases := []struct { + name string + order int + fill func(aQ16 []int32) + }{ + {"zeros_order10", 10, func(aQ16 []int32) {}}, + {"zeros_order16", 16, func(aQ16 []int32) {}}, + {"huge_order10", 10, func(aQ16 []int32) { + for i := range aQ16 { + aQ16[i] = -(1 << 22) + } + }}, + {"huge_order16", 16, func(aQ16 []int32) { + for i := range aQ16 { + aQ16[i] = -(1 << 22) + } + }}, + } + + for _, tc := range cases { + t.Run(tc.name, func(t *testing.T) { + nlsf := make([]int16, tc.order) + aQ16 := make([]int32, tc.order) + tc.fill(aQ16) + + a2nlsf(nlsf, aQ16, tc.order) + + // Whatever path was taken, the result must be a valid increasing + // NLSF vector — never NaN-equivalent garbage or negative values. + for k := range nlsf { + require.GreaterOrEqual(t, nlsf[k], int16(0), "coefficient %d negative", k) + } + for k := 1; k < tc.order; k++ { + require.Greaterf(t, nlsf[k], nlsf[k-1], "not increasing at %d", k) + } + }) + } +}