Add SILK floating-point LPC primitives

internal/silk/lpc_float.go

@@ -0,0 +1,121 @@
+// SPDX-FileCopyrightText: 2026 The Pion community <https://pion.ly>
+// SPDX-License-Identifier: MIT
+
+package silk
+
+import "math"
+
+// Floating-point LPC and signal-analysis primitives shared by the pitch and
+// LPC analysis stages (silk/float/*_FLP.c).
+
+// maxLPCOrder matches libopus's SILK_MAX_ORDER_LPC (SigProc_FIX.h), which
+// covers hybrid/SWB orders this port doesn't implement yet — this encoder
+// only ever requests order 10 (NB/MB) or 16 (WB).
+const maxLPCOrder = 24
+
+// innerProductFLP accumulates a dot product in double precision.
+func innerProductFLP(a, b []float32, n int) float64 {
+	var acc float64
+	for i := range n {
+		acc += float64(a[i]) * float64(b[i])
+	}
+
+	return acc
+}
+
+// autocorrelationFLP computes the first correlationCount autocorrelation taps.
+func autocorrelationFLP(results, inputData []float32, inputDataSize, correlationCount int) {
+	if correlationCount > inputDataSize {
+		correlationCount = inputDataSize
+	}
+	for i := range correlationCount {
+		results[i] = float32(innerProductFLP(inputData, inputData[i:], inputDataSize-i))
+	}
+}
+
+// schurFLP computes reflection coefficients from an autocorrelation sequence
+// and returns the residual energy (silk_schur_FLP).
+func schurFLP(reflCoef, autoCorr []float32, order int) float32 {
+	var c [maxLPCOrder + 1][2]float64 //nolint:varnamelen // c is the correlation work matrix, as in the C reference.
+	for k := 0; k <= order; k++ {
+		c[k][0] = float64(autoCorr[k])
+		c[k][1] = float64(autoCorr[k])
+	}
+	for k := range order {
+		rcTmp := -c[k+1][0] / math.Max(c[0][1], 1e-9)
+		reflCoef[k] = float32(rcTmp)
+		for n := range order - k {
+			ctmp1 := c[n+k+1][0]
+			ctmp2 := c[n][1]
+			c[n+k+1][0] = ctmp1 + ctmp2*rcTmp
+			c[n][1] = ctmp2 + ctmp1*rcTmp
+		}
+	}
+
+	return float32(c[0][1])
+}
+
+// k2aFLP converts reflection coefficients to LPC prediction coefficients.
+func k2aFLP(a, rc []float32, order int) {
+	for k := range order {
+		rck := rc[k] //nolint:gosec // G602: k < order <= maxLPCOrder.
+		for n := range (k + 1) >> 1 {
+			tmp1 := a[n] //nolint:gosec // G602: indices bounded by k < order.
+			tmp2 := a[k-n-1]
+			a[n] = tmp1 + tmp2*rck //nolint:gosec // G602: indices bounded by k < order.
+			a[k-n-1] = tmp2 + tmp1*rck
+		}
+		a[k] = -rck //nolint:gosec // G602: k < order <= maxLPCOrder.
+	}
+}
+
+// bwexpanderFLP applies bandwidth expansion (chirp) to an AR filter.
+func bwexpanderFLP(ar []float32, d int, chirp float32) {
+	cfac := chirp
+	for i := range d - 1 {
+		ar[i] *= cfac
+		cfac *= chirp
+	}
+	ar[d-1] *= cfac
+}
+
+// applySineWindowFLP multiplies px by a sine window. winType 1 starts at 0
+// (rising edge); winType 2 starts at 1 (falling edge). length must be a
+// multiple of 4.
+func applySineWindowFLP(pxWin, px []float32, winType, length int) {
+	freq := float32(math.Pi) / float32(length+1)
+	c := 2.0 - freq*freq //nolint:varnamelen // c is the recurrence coefficient, as in the C reference.
+
+	var s0, s1 float32
+	if winType < 2 {
+		s0 = 0.0
+		s1 = freq
+	} else {
+		s0 = 1.0
+		s1 = 0.5 * c
+	}
+
+	for k := 0; k < length; k += 4 {
+		pxWin[k+0] = px[k+0] * 0.5 * (s0 + s1)
+		pxWin[k+1] = px[k+1] * s1 //nolint:gosec // G602: length is a multiple of 4.
+		s0 = c*s1 - s0
+		pxWin[k+2] = px[k+2] * 0.5 * (s1 + s0) //nolint:gosec // G602: length is a multiple of 4.
+		pxWin[k+3] = px[k+3] * s0              //nolint:gosec // G602: length is a multiple of 4.
+		s1 = c*s0 - s1
+	}
+}
+
+// lpcAnalysisFilterFLP computes the LPC residual r = s - predicted(s). The
+// first order samples of r are set to zero (silk_LPC_analysis_filter_FLP).
+func lpcAnalysisFilterFLP(rLPC, predCoef, s []float32, length, order int) {
+	for ix := order; ix < length; ix++ {
+		var pred float32
+		for j := range order {
+			pred += s[ix-1-j] * predCoef[j]
+		}
+		rLPC[ix] = s[ix] - pred
+	}
+	for i := range order {
+		rLPC[i] = 0
+	}
+}

internal/silk/lpc_float_test.go

@@ -0,0 +1,135 @@
+// 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 TestAutocorrelationFLP(t *testing.T) {
+	in := []float32{1, 2, 3, 4}
+	results := make([]float32, 3)
+	autocorrelationFLP(results, in, len(in), 3)
+
+	assert.InDelta(t, 1+4+9+16, results[0], 1e-4)    // sum of squares
+	assert.InDelta(t, 1*2+2*3+3*4, results[1], 1e-4) // lag 1
+	assert.InDelta(t, 1*3+2*4, results[2], 1e-4)     // lag 2
+}
+
+// TestAutocorrelationFLPClampsCorrelationCount exercises the
+// correlationCount > inputDataSize clamp: without it, the lag-i inner
+// product would slice inputData past its end and panic once i exceeds
+// inputDataSize.
+func TestAutocorrelationFLPClampsCorrelationCount(t *testing.T) {
+	in := []float32{1, 2, 3, 4}
+	results := make([]float32, 10)
+
+	require.NotPanics(t, func() {
+		autocorrelationFLP(results, in, len(in), 10)
+	})
+
+	// Clamped to len(in): only the first 4 lags get computed, same values
+	// as requesting exactly that many up front.
+	assert.InDelta(t, 1+4+9+16, results[0], 1e-4)
+	assert.InDelta(t, 1*2+2*3+3*4, results[1], 1e-4)
+	assert.InDelta(t, 1*3+2*4, results[2], 1e-4)
+	assert.InDelta(t, 1*4, results[3], 1e-4)
+	for i := 4; i < len(results); i++ {
+		assert.Equal(t, float32(0), results[i], "lags beyond inputDataSize should be left untouched")
+	}
+}
+
+func TestBwexpanderFLP(t *testing.T) {
+	ar := []float32{1, 1, 1, 1}
+	bwexpanderFLP(ar, len(ar), 0.5)
+	assert.InDelta(t, 0.5, ar[0], 1e-6)
+	assert.InDelta(t, 0.25, ar[1], 1e-6)
+	assert.InDelta(t, 0.125, ar[2], 1e-6)
+	assert.InDelta(t, 0.0625, ar[3], 1e-6)
+}
+
+func TestK2ASingleCoefficient(t *testing.T) {
+	a := make([]float32, 1)
+	k2aFLP(a, []float32{0.5}, 1)
+	assert.InDelta(t, -0.5, a[0], 1e-6)
+}
+
+// TestSchurThenK2A checks the Schur/k2a pair reproduces a known stable AR
+// filter from its autocorrelation.
+func TestSchurThenK2A(t *testing.T) {
+	// AR(1) process x[n] = a*x[n-1] + e: autocorr r[k] = r[0]*a^|k|.
+	const a = 0.8 //nolint:varnamelen // a is the AR(1) coefficient in the test model.
+	order := 4
+	autoCorr := make([]float32, order+1)
+	for k := range autoCorr {
+		autoCorr[k] = float32(math.Pow(a, float64(k)))
+	}
+
+	refl := make([]float32, order)
+	resNrg := schurFLP(refl, autoCorr, order)
+
+	// First reflection coefficient of an AR(1) process is -a.
+	assert.InDelta(t, -a, refl[0], 1e-4)
+	assert.Greater(t, resNrg, float32(0))
+	assert.LessOrEqual(t, resNrg, autoCorr[0])
+
+	coef := make([]float32, order)
+	k2aFLP(coef, refl, order)
+	// Recovered LPC: first tap ~ a, the rest ~ 0.
+	assert.InDelta(t, a, coef[0], 1e-3)
+	for k := 1; k < order; k++ {
+		assert.InDelta(t, 0, coef[k], 1e-3)
+	}
+}
+
+func TestApplySineWindowFLP(t *testing.T) {
+	length := 16
+	px := make([]float32, length)
+	for i := range px {
+		px[i] = 1
+	}
+
+	rising := make([]float32, length)
+	applySineWindowFLP(rising, px, 1, length)
+	// Rising window: starts near 0, increases, ends near the peak.
+	assert.InDelta(t, 0, rising[0], 0.1)
+	assert.Greater(t, rising[length-1], rising[0])
+	for _, v := range rising {
+		assert.GreaterOrEqual(t, v, float32(-1e-3))
+		assert.LessOrEqual(t, v, float32(1.001))
+	}
+
+	falling := make([]float32, length)
+	applySineWindowFLP(falling, px, 2, length)
+	// Falling window: starts near the peak, ends near 0.
+	assert.Greater(t, falling[0], falling[length-1])
+	assert.InDelta(t, 0, falling[length-1], 0.1)
+}
+
+func TestLPCAnalysisFilterFLP(t *testing.T) {
+	s := []float32{1, 2, 3, 4, 5, 6, 7, 8} //nolint:varnamelen // s is the input signal.
+	order := 2
+	r := make([]float32, len(s)) //nolint:varnamelen // r is the residual buffer.
+
+	// Zero predictor: residual equals the input past the first order samples.
+	lpcAnalysisFilterFLP(r, []float32{0, 0}, s, len(s), order)
+	for i := range order {
+		assert.Equal(t, float32(0), r[i])
+	}
+	for i := order; i < len(s); i++ {
+		assert.InDelta(t, s[i], r[i], 1e-6)
+	}
+
+	// A known predictor removes its own prediction.
+	coef := []float32{1.5, -0.5}
+	lpcAnalysisFilterFLP(r, coef, s, len(s), order)
+	for i := order; i < len(s); i++ {
+		want := s[i] - (coef[0]*s[i-1] + coef[1]*s[i-2])
+		require.InDelta(t, want, r[i], 1e-4)
+	}
+}