From 7fa92ea3463b4d11c4bd64361e42b040b374fb6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Allais?= Date: Fri, 17 Jul 2026 18:10:09 -0300 Subject: [PATCH 1/2] Add SILK floating-point LPC primitives --- internal/silk/lpc_float.go | 118 ++++++++++++++++++++++++++++++++ internal/silk/lpc_float_test.go | 112 ++++++++++++++++++++++++++++++ 2 files changed, 230 insertions(+) create mode 100644 internal/silk/lpc_float.go create mode 100644 internal/silk/lpc_float_test.go diff --git a/internal/silk/lpc_float.go b/internal/silk/lpc_float.go new file mode 100644 index 0000000..63cf1f9 --- /dev/null +++ b/internal/silk/lpc_float.go @@ -0,0 +1,118 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// 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). + +const maxLPCOrder = 16 // SILK_MAX_ORDER_LPC + +// 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 + } +} diff --git a/internal/silk/lpc_float_test.go b/internal/silk/lpc_float_test.go new file mode 100644 index 0000000..379a7e9 --- /dev/null +++ b/internal/silk/lpc_float_test.go @@ -0,0 +1,112 @@ +// 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 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 +} + +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) + } +} From 58af8356911b060b673021165e490474c245a424 Mon Sep 17 00:00:00 2001 From: Thomas Vilte Date: Sun, 19 Jul 2026 15:18:12 -0300 Subject: [PATCH 2/2] Fix maxLPCOrder, add autocorrelation clamp test MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit maxLPCOrder was 16 with a comment claiming it matched libopus's SILK_MAX_ORDER_LPC, but that constant is actually 24 — harmless today since this port only requests order 10/16, but a latent trap if a future hybrid/SWB sub-PR reuses schurFLP at a higher order. Also covers the correlationCount > inputDataSize clamp in autocorrelationFLP, untested before (75% coverage on that function) — without it, the lag loop slices inputData past its end once the requested count exceeds what's available. Co-authored-by: François Allais --- internal/silk/lpc_float.go | 5 ++++- internal/silk/lpc_float_test.go | 23 +++++++++++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/internal/silk/lpc_float.go b/internal/silk/lpc_float.go index 63cf1f9..7f212af 100644 --- a/internal/silk/lpc_float.go +++ b/internal/silk/lpc_float.go @@ -8,7 +8,10 @@ import "math" // Floating-point LPC and signal-analysis primitives shared by the pitch and // LPC analysis stages (silk/float/*_FLP.c). -const maxLPCOrder = 16 // SILK_MAX_ORDER_LPC +// 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 { diff --git a/internal/silk/lpc_float_test.go b/internal/silk/lpc_float_test.go index 379a7e9..2cdf02c 100644 --- a/internal/silk/lpc_float_test.go +++ b/internal/silk/lpc_float_test.go @@ -21,6 +21,29 @@ func TestAutocorrelationFLP(t *testing.T) { 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)