Add SILK LTP analysis and quantization

internal/silk/ltp_analysis.go

@@ -0,0 +1,72 @@
+// SPDX-FileCopyrightText: 2026 The Pion community <https://pion.ly>
+// SPDX-License-Identifier: MIT
+
+package silk
+
+const (
+	ltpOrder       = 5
+	ltpCorrInvMax  = 0.03
+	ltpMatrixSize  = ltpOrder * ltpOrder
+	ltpHalfOrder   = ltpOrder / 2
+	ltpLastDiagIdx = ltpMatrixSize - 1
+)
+
+// corrVectorFLP computes the cross-correlation X'*t (silk_corrVector_FLP).
+func corrVectorFLP(x, t []float32, l, order int, xt []float32) {
+	for lag := range order {
+		xt[lag] = float32(innerProductFLP(x[order-1-lag:], t, l))
+	}
+}
+
+// corrMatrixFLP computes the symmetric correlation matrix X'*X, stored row-major
+// (silk_corrMatrix_FLP).
+func corrMatrixFLP(x []float32, l, order int, xx []float32) { //nolint:varnamelen // l is the vector length.
+	p1 := order - 1
+	energy := energyFLP(x[p1:], l)
+	xx[0] = float32(energy)
+	for j := 1; j < order; j++ {
+		energy += float64(x[p1-j])*float64(x[p1-j]) - float64(x[p1+l-j])*float64(x[p1+l-j])
+		xx[j*order+j] = float32(energy)
+	}
+
+	p2 := order - 2
+	for lag := 1; lag < order; lag++ {
+		energy = innerProductFLP(x[p1:], x[p2:], l)
+		xx[lag*order] = float32(energy)
+		xx[lag] = float32(energy)
+		for j := 1; j < order-lag; j++ {
+			energy += float64(x[p1-j])*float64(x[p2-j]) - float64(x[p1+l-j])*float64(x[p2+l-j])
+			xx[(lag+j)*order+j] = float32(energy)
+			xx[j*order+(lag+j)] = float32(energy)
+		}
+		p2--
+	}
+}
+
+// scaleVectorFLP multiplies each element by gain.
+func scaleVectorFLP(data []float32, gain float32, size int) {
+	for i := range size {
+		data[i] *= gain
+	}
+}
+
+// findLTPFLP builds the normalized LTP correlation matrices (XX) and vectors
+// (xX) for each subframe (silk_find_LTP_FLP). r is the LPC residual and rOffset
+// points at the start of the current frame within it (with history before).
+func findLTPFLP(xx, xX, r []float32, rOffset int, lag []int, subfrLength, nbSubfr int) {
+	xxPtr, xXPtr, rPtr := 0, 0, rOffset
+	for k := range nbSubfr {
+		lagPtr := rPtr - (lag[k] + ltpHalfOrder)
+		corrMatrixFLP(r[lagPtr:], subfrLength, ltpOrder, xx[xxPtr:])
+		corrVectorFLP(r[lagPtr:], r[rPtr:], subfrLength, ltpOrder, xX[xXPtr:])
+
+		energy := float32(energyFLP(r[rPtr:], subfrLength+ltpOrder))
+		temp := 1.0 / max(energy, ltpCorrInvMax*0.5*(xx[xxPtr]+xx[xxPtr+ltpLastDiagIdx])+1.0)
+		scaleVectorFLP(xx[xxPtr:], temp, ltpMatrixSize)
+		scaleVectorFLP(xX[xXPtr:], temp, ltpOrder)
+
+		rPtr += subfrLength
+		xxPtr += ltpMatrixSize
+		xXPtr += ltpOrder
+	}
+}

internal/silk/ltp_analysis_test.go

@@ -0,0 +1,77 @@
+// SPDX-FileCopyrightText: 2026 The Pion community <https://pion.ly>
+// SPDX-License-Identifier: MIT
+
+package silk
+
+import (
+	"math"
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+	"github.com/stretchr/testify/require"
+)
+
+func TestCorrVectorFLP(t *testing.T) {
+	order := 3
+	l := 4
+	x := []float32{1, 2, 3, 4, 5, 6} // length l+order-1
+	target := []float32{1, 1, 1, 1}
+	xt := make([]float32, order)
+	corrVectorFLP(x, target, l, order, xt)
+
+	// lag 0 sums x[order-1 .. order-1+l) = x[2..6) = 3+4+5+6
+	assert.InDelta(t, 3+4+5+6, xt[0], 1e-4)
+	assert.InDelta(t, 2+3+4+5, xt[1], 1e-4)
+	assert.InDelta(t, 1+2+3+4, xt[2], 1e-4)
+}
+
+func TestCorrMatrixFLP(t *testing.T) {
+	order := 3
+	l := 4
+	x := []float32{1, 2, 3, 4, 5, 6}
+	xx := make([]float32, order*order)
+	corrMatrixFLP(x, l, order, xx)
+
+	// Diagonal [0][0] = energy of x[2..6) = 9+16+25+36 = 86.
+	assert.InDelta(t, 86, xx[0], 1e-3)
+	// Symmetric.
+	for i := range order {
+		for j := range order {
+			require.InDeltaf(t, xx[i*order+j], xx[j*order+i], 1e-3, "asymmetry at %d,%d", i, j)
+		}
+	}
+}
+
+func TestFindLTPFLP(t *testing.T) {
+	const (
+		nbSubfr     = 4
+		subfrLength = 40
+		maxLag      = 120
+	)
+	// Residual buffer with history ahead of the frame start.
+	total := maxLag + ltpOrder + nbSubfr*subfrLength + ltpOrder
+	r := make([]float32, total)
+	state := uint32(2024)
+	for i := range r {
+		state = 1664525*state + 1013904223
+		r[i] = float32(int32(state>>16)%2000-1000) / 1000
+	}
+	rOffset := maxLag + ltpOrder
+	lag := []int{80, 82, 78, 81}
+
+	xx := make([]float32, nbSubfr*ltpMatrixSize)
+	xX := make([]float32, nbSubfr*ltpOrder)
+	findLTPFLP(xx, xX, r, rOffset, lag, subfrLength, nbSubfr)
+
+	// Each subframe matrix is finite and symmetric.
+	for k := range nbSubfr {
+		base := k * ltpMatrixSize
+		for i := range ltpOrder {
+			for j := range ltpOrder {
+				v := xx[base+i*ltpOrder+j]
+				require.Falsef(t, math.IsNaN(float64(v)) || math.IsInf(float64(v), 0), "non-finite at subfr %d", k)
+				require.InDeltaf(t, xx[base+i*ltpOrder+j], xx[base+j*ltpOrder+i], 1e-5, "asymmetry subfr %d", k)
+			}
+		}
+	}
+}

internal/silk/ltp_quant.go

@@ -0,0 +1,143 @@
+// SPDX-FileCopyrightText: 2026 The Pion community <https://pion.ly>
+// SPDX-License-Identifier: MIT
+
+package silk
+
+import "math"
+
+// Long-term-predictor gain quantization (silk_quant_LTP_gains + silk_VQ_WMat_EC).
+// Given the find_LTP correlation matrices it searches the three LTP codebooks
+// for the periodicity index and per-subframe filter indices that minimize a
+// weighted quantization error plus rate cost.
+
+const nLTPCodebooks = 3
+
+//nolint:gochecknoglobals // LTP codebook effective gains and bit costs (tables_LTP.c).
+var (
+	ltpVQSizes = [nLTPCodebooks]int{8, 16, 32}
+
+	ltpGainVQ0 = []uint8{46, 2, 90, 87, 93, 91, 82, 98}
+	ltpGainVQ1 = []uint8{109, 120, 118, 12, 113, 115, 117, 119, 99, 59, 87, 111, 63, 111, 112, 80}
+	ltpGainVQ2 = []uint8{
+		126, 124, 125, 124, 129, 121, 126, 23, 132, 127, 127, 127, 126, 127, 122, 133,
+		130, 134, 101, 118, 119, 145, 126, 86, 124, 120, 123, 119, 170, 173, 107, 109,
+	}
+
+	ltpBitsQ50 = []uint8{15, 131, 138, 138, 155, 155, 173, 173}
+	ltpBitsQ51 = []uint8{69, 93, 115, 118, 131, 138, 141, 138, 150, 150, 155, 150, 155, 160, 166, 160}
+	ltpBitsQ52 = []uint8{
+		131, 128, 134, 141, 141, 141, 145, 145, 145, 150, 155, 155, 155, 155, 160, 160,
+		160, 160, 166, 166, 173, 173, 182, 192, 182, 192, 192, 192, 205, 192, 205, 224,
+	}
+)
+
+func ltpCodebook(k int) [][]int8 {
+	switch k {
+	case 1:
+		return codebookLTPFilterPeriodicityIndex1
+	case 2:
+		return codebookLTPFilterPeriodicityIndex2
+	default:
+		return codebookLTPFilterPeriodicityIndex0
+	}
+}
+
+func ltpGainTable(k int) []uint8 {
+	switch k {
+	case 1:
+		return ltpGainVQ1
+	case 2:
+		return ltpGainVQ2
+	default:
+		return ltpGainVQ0
+	}
+}
+
+func ltpBitsTable(k int) []uint8 {
+	switch k {
+	case 1:
+		return ltpBitsQ51
+	case 2:
+		return ltpBitsQ52
+	default:
+		return ltpBitsQ50
+	}
+}
+
+// log2lin approximates 2^(inLogQ7/128) (silk_log2lin).
+func log2lin(inLogQ7 int32) int32 {
+	if inLogQ7 < 0 {
+		return 0
+	}
+
+	if inLogQ7 >= 3967 {
+		return math.MaxInt32
+	}
+	out := int32(1) << (inLogQ7 >> 7)
+	frac := inLogQ7 & 0x7F
+	adj := smlawb(frac, smulbb(frac, 128-frac), -174)
+	if inLogQ7 < 2048 {
+		return out + ((out * adj) >> 7)
+	}
+
+	return out + (out>>7)*adj
+}
+
+// vqWMatEC finds the codebook vector minimizing the weighted quantization error
+// plus rate (silk_VQ_WMat_EC). XX_Q17 is the 5x5 correlation matrix (row-major),
+// xX_Q17 the correlation vector, both Q17.
+func vqWMatEC(
+	xxQ17, xXQ17 []int32, cb [][]int8, cbGain, clQ5 []uint8, subfrLen int, maxGainQ7 int32, l int,
+) (ind int, resNrgQ15, rateDistQ8, gainQ7 int32) {
+	var negXX [ltpOrder]int32
+	for i := range ltpOrder {
+		negXX[i] = -(xXQ17[i] << 7)
+	}
+
+	rateDistQ8 = math.MaxInt32
+	resNrgQ15 = math.MaxInt32
+	for k := range l { //nolint:varnamelen // k indexes the codebook, as in the C reference.
+		row := cb[k]
+		gainTmp := int32(cbGain[k])
+		sum1 := int32(32801) // 1.001 in Q15
+		penalty := max(gainTmp-maxGainQ7, 0) << 11
+
+		sum2 := negXX[0] + xxQ17[1]*int32(row[1]) + xxQ17[2]*int32(row[2]) + xxQ17[3]*int32(row[3]) + xxQ17[4]*int32(row[4])
+		sum2 = (sum2 << 1) + xxQ17[0]*int32(row[0])
+		sum1 = smlawb(sum1, sum2, int32(row[0]))
+
+		sum2 = negXX[1] + xxQ17[7]*int32(row[2]) + xxQ17[8]*int32(row[3]) + xxQ17[9]*int32(row[4])
+		sum2 = (sum2 << 1) + xxQ17[6]*int32(row[1])
+		sum1 = smlawb(sum1, sum2, int32(row[1]))
+
+		sum2 = negXX[2] + xxQ17[13]*int32(row[3]) + xxQ17[14]*int32(row[4])
+		sum2 = (sum2 << 1) + xxQ17[12]*int32(row[2])
+		sum1 = smlawb(sum1, sum2, int32(row[2]))
+
+		sum2 = negXX[3] + xxQ17[19]*int32(row[4])
+		sum2 = (sum2 << 1) + xxQ17[18]*int32(row[3])
+		sum1 = smlawb(sum1, sum2, int32(row[3]))
+
+		sum2 = (negXX[4] << 1) + xxQ17[24]*int32(row[4])
+		sum1 = smlawb(sum1, sum2, int32(row[4]))
+
+		if sum1 >= 0 {
+			bitsResQ8 := smulbb(int32(subfrLen), lin2log(sum1+penalty)-(15<<7)) //nolint:gosec // G115
+			bitsTotQ8 := addLShift32(bitsResQ8, int32(clQ5[k]), 2)
+			if bitsTotQ8 <= rateDistQ8 {
+				rateDistQ8 = bitsTotQ8
+				resNrgQ15 = sum1 + penalty
+				ind = k
+				gainQ7 = gainTmp
+			}
+		}
+	}
+
+	return ind, resNrgQ15, rateDistQ8, gainQ7
+}
+
+// quantLTPGains (silk_quant_LTP_gains) is deferred to the orchestration piece:
+// it's a method on *Encoder (threads sumLogGainQ7 across frames), and that
+// type doesn't exist in this package yet — same reasoning as
+// findPitchLags/encodeNLSF. What's here is the per-subframe codebook search
+// (vqWMatEC) it calls in a loop over the 3 LTP codebooks.

internal/silk/ltp_quant_test.go

@@ -0,0 +1,71 @@
+// SPDX-FileCopyrightText: 2026 The Pion community <https://pion.ly>
+// SPDX-License-Identifier: MIT
+
+package silk
+
+import (
+	"math"
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+	"github.com/stretchr/testify/require"
+)
+
+// TestLog2Lin checks log2lin against its float inverse, math.Log2, over a
+// representative Q7 range (silk_log2lin approximates 2^(x/128)).
+func TestLog2Lin(t *testing.T) {
+	assert.Zero(t, log2lin(-1))
+	assert.Equal(t, int32(math.MaxInt32), log2lin(3967))
+
+	for _, inLogQ7 := range []int32{0, 128, 256, 1000, 2048, 3000} {
+		got := float64(log2lin(inLogQ7))
+		want := math.Pow(2, float64(inLogQ7)/128.0)
+		assert.InEpsilonf(t, want, got, 0.03, "log2lin(%d)", inLogQ7)
+	}
+}
+
+// TestVQWMatEC feeds a real (not synthetic) correlation matrix/vector — built
+// from an actual signal via corrMatrixFLP/corrVectorFLP, the same way
+// quantLTPGains will once it's wired to *Encoder — and checks the returned
+// index, energy, and gain are self-consistent, exercising the fixed-point
+// codebook search on physically meaningful input rather than arbitrary noise.
+func TestVQWMatEC(t *testing.T) {
+	const (
+		subfrLength = 40
+		order       = ltpOrder
+	)
+	// A short periodic-ish signal so the correlation matrix isn't degenerate.
+	x := make([]float32, subfrLength+order-1)
+	for i := range x {
+		x[i] = float32(10 * math.Sin(2*math.Pi*float64(i)/17))
+	}
+
+	xxFLP := make([]float32, order*order)
+	xXFLP := make([]float32, order)
+	corrMatrixFLP(x, subfrLength, order, xxFLP)
+	corrVectorFLP(x, x[order-1:], subfrLength, order, xXFLP)
+
+	xxQ17 := make([]int32, order*order)
+	for i, v := range xxFLP {
+		xxQ17[i] = int32(math.RoundToEven(float64(v) * 131072.0))
+	}
+	xXQ17 := make([]int32, order)
+	for i, v := range xXFLP {
+		xXQ17[i] = int32(math.RoundToEven(float64(v) * 131072.0))
+	}
+
+	for k := range nLTPCodebooks {
+		cb := ltpCodebook(k)
+		cbGain := ltpGainTable(k)
+		clQ5 := ltpBitsTable(k)
+		size := ltpVQSizes[k]
+
+		ind, resNrgQ15, rateDistQ8, gainQ7 := vqWMatEC(xxQ17, xXQ17, cb, cbGain, clQ5, subfrLength, 1<<20, size)
+
+		require.GreaterOrEqualf(t, ind, 0, "codebook %d index", k)
+		require.Lessf(t, ind, size, "codebook %d index", k)
+		assert.NotEqual(t, int32(math.MaxInt32), resNrgQ15, "codebook %d: no candidate scored", k)
+		assert.NotEqual(t, int32(math.MaxInt32), rateDistQ8, "codebook %d: no candidate scored", k)
+		assert.Equal(t, int32(cbGain[ind]), gainQ7, "codebook %d: gain matches winning row", k)
+	}
+}