diff --git a/internal/silk/ltp_analysis.go b/internal/silk/ltp_analysis.go new file mode 100644 index 0000000..b637b5f --- /dev/null +++ b/internal/silk/ltp_analysis.go @@ -0,0 +1,72 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// 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 + } +} diff --git a/internal/silk/ltp_analysis_test.go b/internal/silk/ltp_analysis_test.go new file mode 100644 index 0000000..6db70cb --- /dev/null +++ b/internal/silk/ltp_analysis_test.go @@ -0,0 +1,77 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// 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) + } + } + } +} diff --git a/internal/silk/ltp_quant.go b/internal/silk/ltp_quant.go new file mode 100644 index 0000000..d9e8dee --- /dev/null +++ b/internal/silk/ltp_quant.go @@ -0,0 +1,143 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// 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. diff --git a/internal/silk/ltp_quant_test.go b/internal/silk/ltp_quant_test.go new file mode 100644 index 0000000..84358e7 --- /dev/null +++ b/internal/silk/ltp_quant_test.go @@ -0,0 +1,71 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// 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) + } +} diff --git a/internal/silk/ltp_scale.go b/internal/silk/ltp_scale.go new file mode 100644 index 0000000..c52a084 --- /dev/null +++ b/internal/silk/ltp_scale.go @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// LTP state-scaling gains (silk_LTPScales_table_Q14). +// +//nolint:gochecknoglobals // constant lookup table. +var ltpScalesTableQ14 = [3]int32{15565, 12288, 8192} + +// ltpScaleControl selects the LTP state-scaling index for an independently +// coded voiced frame (silk_LTP_scale_ctrl_FLP). ltpredCodGain is the LTP +// prediction gain in dB and snrDBQ7 the target SNR in Q7. Higher expected loss +// pushes toward stronger scaling so decoded frames recover faster; with no +// configured packet loss the index is always 0. Returns the index and its Q14 +// scale. +// +//nolint:gosec // G115: dB gains and loss percentages are small, in-range values. +func ltpScaleControl( + ltpredCodGain float32, snrDBQ7 int32, packetLossPerc, nFramesPerPacket int, lbrr bool, +) (int, int32) { + roundLoss := packetLossPerc * nFramesPerPacket + if lbrr { + roundLoss = 2 + int(smulbb(int32(roundLoss), int32(roundLoss)))/100 + } + + idx := 0 + if smulbb(int32(ltpredCodGain), int32(roundLoss)) > log2lin(2900-snrDBQ7) { + idx++ + } + if smulbb(int32(ltpredCodGain), int32(roundLoss)) > log2lin(3900-snrDBQ7) { + idx++ + } + + return idx, ltpScalesTableQ14[idx] +} diff --git a/internal/silk/ltp_scale_test.go b/internal/silk/ltp_scale_test.go new file mode 100644 index 0000000..297cc0c --- /dev/null +++ b/internal/silk/ltp_scale_test.go @@ -0,0 +1,35 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" +) + +func TestLTPScaleControl(t *testing.T) { + // No packet loss: always minimum scaling (index 0, 15565), regardless of + // the target SNR — roundLoss is 0, so the product with it never crosses + // either threshold. + idx, q14 := ltpScaleControl(20, 10*128, 0, 1, false) + assert.Equal(t, 0, idx, "no loss: index") + assert.Equal(t, int32(15565), q14, "no loss: scale") + + // High prediction gain and loss push toward stronger scaling. + idx, q14 = ltpScaleControl(40, 25*128, 25, 1, false) + assert.NotEqual(t, 0, idx, "high loss/gain: expected stronger scaling") + assert.Equal(t, ltpScalesTableQ14[idx], q14, "scale matches table entry") + + // LBRR (low bitrate redundancy) present: roundLoss uses the squared-loss + // formula instead of the plain product, still with no plain packet loss. + idx, q14 = ltpScaleControl(40, 18*128, 0, 1, true) + assert.NotEqual(t, 0, idx, "lbrr: expected stronger scaling") + assert.Equal(t, ltpScalesTableQ14[idx], q14, "scale matches table entry") + + // Extreme gain and loss cross both thresholds: strongest scaling (index 2). + idx, q14 = ltpScaleControl(100, 30*128, 100, 3, false) + assert.Equal(t, 2, idx, "extreme loss/gain: index") + assert.Equal(t, int32(8192), q14, "extreme loss/gain: scale") +}