From fe0e821f42ec9844b31b900d2cf6d82d2ab60972 Mon Sep 17 00:00:00 2001 From: Thomas Vilte Date: Tue, 21 Jul 2026 15:42:08 -0300 Subject: [PATCH] Add SILK Burg LPC analysis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: François Allais --- internal/silk/lpc_analysis.go | 142 +++++++++++++++++++++++++++++ internal/silk/lpc_analysis_test.go | 82 +++++++++++++++++ 2 files changed, 224 insertions(+) create mode 100644 internal/silk/lpc_analysis.go create mode 100644 internal/silk/lpc_analysis_test.go diff --git a/internal/silk/lpc_analysis.go b/internal/silk/lpc_analysis.go new file mode 100644 index 0000000..79143b1 --- /dev/null +++ b/internal/silk/lpc_analysis.go @@ -0,0 +1,142 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +const findLPCCondFac = 1e-5 + +// findLPC estimates LPC coefficients with Burg and converts them to NLSFs +// (silk_find_LPC_FLP). The optional NLSF-interpolation search is deferred; the +// interpolation index stays at 4 (no interpolation). x holds nb_subfr +// subframes of subfrLength samples each (including order preceding samples). +func findLPC(nlsfQ15 []int16, x []float32, minInvGain float32, subfrLength, nbSubfr, order int) { + a := make([]float32, order) + burgModifiedFLP(a, x, minInvGain, subfrLength, nbSubfr, order) + a2nlsfFLP(nlsfQ15, a, order) +} + +// burgModifiedFLP estimates LPC coefficients with Burg's method, stacked over +// nb_subfr subframes, limiting the prediction gain to 1/minInvGain. It writes +// the order coefficients into a and returns the residual energy +// (silk_burg_modified_FLP). +// +//nolint:gocognit,gocyclo,cyclop,maintidx,varnamelen // faithful port of silk_burg_modified_FLP. +func burgModifiedFLP(a, x []float32, minInvGain float32, subfrLength, nbSubfr, order int) float32 { + var cFirstRow, cLastRow [maxLPCOrder]float64 + var cAf, cAb [maxLPCOrder + 1]float64 + var af [maxLPCOrder]float64 + + c0 := energyFLP(x, nbSubfr*subfrLength) + for s := range nbSubfr { + xPtr := x[s*subfrLength:] + for n := 1; n < order+1; n++ { + cFirstRow[n-1] += innerProductFLP(xPtr, xPtr[n:], subfrLength-n) + } + } + copy(cLastRow[:], cFirstRow[:]) + + cAb[0] = c0 + findLPCCondFac*c0 + 1e-9 + cAf[0] = cAb[0] + invGain := 1.0 + reachedMaxGain := false + + for n := range order { + for s := range nbSubfr { + xPtr := x[s*subfrLength:] + tmp1 := float64(xPtr[n]) + tmp2 := float64(xPtr[subfrLength-n-1]) + for k := range n { + cFirstRow[k] -= float64(xPtr[n]) * float64(xPtr[n-k-1]) + cLastRow[k] -= float64(xPtr[subfrLength-n-1]) * float64(xPtr[subfrLength-n+k]) + atmp := af[k] + tmp1 += float64(xPtr[n-k-1]) * atmp + tmp2 += float64(xPtr[subfrLength-n+k]) * atmp + } + for k := 0; k <= n; k++ { + cAf[k] -= tmp1 * float64(xPtr[n-k]) + cAb[k] -= tmp2 * float64(xPtr[subfrLength-n+k-1]) + } + } + + tmp1 := cFirstRow[n] + tmp2 := cLastRow[n] + for k := range n { + atmp := af[k] + tmp1 += cLastRow[n-k-1] * atmp + tmp2 += cFirstRow[n-k-1] * atmp + } + cAf[n+1] = tmp1 //nolint:gosec // G602: n+1 <= order <= maxLPCOrder. + cAb[n+1] = tmp2 //nolint:gosec // G602: n+1 <= order <= maxLPCOrder. + + num := cAb[n+1] //nolint:gosec // G602: n+1 <= order <= maxLPCOrder. + nrgB := cAb[0] + nrgF := cAf[0] + for k := range n { + atmp := af[k] + num += cAb[n-k] * atmp + nrgB += cAb[k+1] * atmp + nrgF += cAf[k+1] * atmp + } + + rc := -2.0 * num / (nrgF + nrgB) + + gain := invGain * (1.0 - rc*rc) + if gain <= float64(minInvGain) { + rc = math.Sqrt(1.0 - float64(minInvGain)/invGain) + if num > 0 { + rc = -rc + } + invGain = float64(minInvGain) + reachedMaxGain = true + } else { + invGain = gain + } + + for k := range (n + 1) >> 1 { + t1 := af[k] + t2 := af[n-k-1] + af[k] = t1 + rc*t2 + af[n-k-1] = t2 + rc*t1 + } + af[n] = rc + + if reachedMaxGain { + for k := n + 1; k < order; k++ { + af[k] = 0.0 + } + + break + } + + for k := 0; k <= n+1; k++ { + t1 := cAf[k] + cAf[k] += rc * cAb[n-k+1] + cAb[n-k+1] += rc * t1 + } + } + + var residualEnergy float64 + if reachedMaxGain { + for k := range order { + a[k] = float32(-af[k]) + } + for s := range nbSubfr { + c0 -= energyFLP(x[s*subfrLength:], order) + } + residualEnergy = c0 * invGain + } else { + residualEnergy = cAf[0] + sumSq := 1.0 + for k := range order { + atmp := af[k] + residualEnergy += cAf[k+1] * atmp + sumSq += atmp * atmp + a[k] = float32(-atmp) + } + residualEnergy -= findLPCCondFac * c0 * sumSq + } + + return float32(residualEnergy) +} diff --git a/internal/silk/lpc_analysis_test.go b/internal/silk/lpc_analysis_test.go new file mode 100644 index 0000000..d0b801e --- /dev/null +++ b/internal/silk/lpc_analysis_test.go @@ -0,0 +1,82 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// arSignal generates a length-n AR(2) process x[i] = a1*x[i-1] + a2*x[i-2] + e. +func arSignal(n int, a1, a2 float32, seed uint32) []float32 { + x := make([]float32, n) + var prev1, prev2 float32 + state := seed + for i := range x { + state = 1664525*state + 1013904223 + e := float32(int32(state>>16)%2000-1000) / 1000 + cur := a1*prev1 + a2*prev2 + e + x[i] = cur + prev2, prev1 = prev1, cur + } + + return x +} + +func TestBurgModifiedFLPRecoversAR(t *testing.T) { + const ( + a1 = 0.5 + a2 = 0.2 + n = 512 + ) + x := arSignal(n, a1, a2, 42) + + a := make([]float32, 2) + resNrg := burgModifiedFLP(a, x, 1e-5, n, 1, 2) + + assert.InDelta(t, a1, a[0], 0.12) + assert.InDelta(t, a2, a[1], 0.12) + assert.Greater(t, resNrg, float32(0)) + assert.Less(t, resNrg, float32(energyFLP(x, n)), "AR signal should be predictable") +} + +func TestBurgModifiedFLPWhiteNoise(t *testing.T) { + const n = 512 + x := make([]float32, n) + state := uint32(777) + for i := range x { + state = 1664525*state + 1013904223 + x[i] = float32(int32(state>>16)%2000-1000) / 1000 + } + + a := make([]float32, 4) + burgModifiedFLP(a, x, 1e-5, n, 1, 4) + + // White noise has no predictable structure; coefficients should be small. + for k, c := range a { + require.InDeltaf(t, 0, c, 0.2, "coefficient %d", k) + } +} + +func TestFindLPC(t *testing.T) { + const ( + order = 16 + nbSubfr = 4 + subfrLength = 40 + order + ) + x := arSignal(nbSubfr*subfrLength, 0.6, 0.2, 99) + + nlsf := make([]int16, order) + findLPC(nlsf, x, 1e-4, subfrLength, nbSubfr, order) + + // A valid NLSF vector: non-decreasing and non-negative (Q15, <= int16 max). + for k := range nlsf { + require.GreaterOrEqual(t, nlsf[k], int16(0)) + } + for k := 1; k < order; k++ { + require.GreaterOrEqual(t, nlsf[k], nlsf[k-1]) + } +}