diff --git a/internal/silk/encoder.go b/internal/silk/encoder.go new file mode 100644 index 0000000..acd6f62 --- /dev/null +++ b/internal/silk/encoder.go @@ -0,0 +1,68 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "github.com/pion/opus/internal/rangecoding" + +// Encoder quantizes and range-encodes a single SILK channel. It is the +// counterpart to Decoder and is built up one stage at a time. +// +// State fields match the corresponding Decoder fields and use the same reset +// values so the encoder and decoder stay in sync across an uncoded frame. +type Encoder struct { + rangeEncoder rangecoding.Encoder + + // haveEncoded reports whether a frame has been encoded yet; it selects + // independent gain coding for the first subframe of the first frame. + haveEncoded bool + + // previousLogGain is the running quantized log-gain index carried across + // subframes and frames. + previousLogGain int32 + + // previousLag and isPreviousFrameVoiced carry pitch state across frames, + // selecting relative vs absolute primary-lag coding. + previousLag int + isPreviousFrameVoiced bool + + // firstFrameAfterReset caps the predictor more aggressively on the first + // frame after a reset (find_pred_coefs). + firstFrameAfterReset bool + + // prevNLSFq holds the previous frame's quantized NLSFs (Q15) for LSF + // interpolation. + prevNLSFq []int16 + + // Analysis state for the frame encoder. + vad vadState + nsq *nsqState + targetBitrate int // target bitrate in bps (drives control_SNR) + sumLogGainQ7 int32 // cumulative LTP gain limit (quant_LTP_gains) + ltpCorr float32 // normalized correlation carried across frames + tiltSmth float32 // smoothed spectral tilt (shape state) + harmShapeGainSmth float32 // smoothed harmonic shaping gain (shape state) +} + +// NewEncoder creates a SILK Encoder with its prediction state reset. +func NewEncoder() Encoder { + e := Encoder{vad: newVADState(), nsq: newNSQState()} + e.resetPredictionState() + + return e +} + +// resetPredictionState resets the encoder prediction state. The values must +// match Decoder.resetPredictionState. +func (e *Encoder) resetPredictionState() { + e.haveEncoded = false + e.previousLogGain = 10 + e.previousLag = 100 + e.isPreviousFrameVoiced = false + e.firstFrameAfterReset = true + e.sumLogGainQ7 = 0 + e.prevNLSFq = make([]int16, maxLPCOrder) + if e.targetBitrate == 0 { + e.targetBitrate = 24000 + } +} diff --git a/internal/silk/encoder_test.go b/internal/silk/encoder_test.go new file mode 100644 index 0000000..cfd2471 --- /dev/null +++ b/internal/silk/encoder_test.go @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +func TestNewEncoderDefaults(t *testing.T) { + enc := NewEncoder() + + assert.False(t, enc.haveEncoded) + assert.Equal(t, int32(10), enc.previousLogGain) + assert.Equal(t, 100, enc.previousLag) + assert.False(t, enc.isPreviousFrameVoiced) + assert.True(t, enc.firstFrameAfterReset) + assert.Equal(t, int32(0), enc.sumLogGainQ7) + assert.Len(t, enc.prevNLSFq, maxLPCOrder) + assert.Equal(t, 24000, enc.targetBitrate) + require.NotNil(t, enc.nsq) +} + +func TestResetPredictionStateKeepsExplicitBitrate(t *testing.T) { + enc := NewEncoder() + enc.targetBitrate = 32000 + enc.haveEncoded = true + enc.previousLogGain = 40 + enc.firstFrameAfterReset = false + + enc.resetPredictionState() + + assert.False(t, enc.haveEncoded) + assert.Equal(t, int32(10), enc.previousLogGain) + assert.True(t, enc.firstFrameAfterReset) + // A non-zero bitrate the caller already set is left alone — only the + // zero-value default (uninitialized Encoder) gets the 24000 fallback. + assert.Equal(t, 32000, enc.targetBitrate) +} diff --git a/internal/silk/find_pitch_lags.go b/internal/silk/find_pitch_lags.go new file mode 100644 index 0000000..886ac8e --- /dev/null +++ b/internal/silk/find_pitch_lags.go @@ -0,0 +1,70 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +const ( + findPitchWhiteNoiseFraction = 1e-3 + findPitchBandwidthExpansion = 0.99 + laPitchMS = 2 + findPitchLPCWinMS = 20 + laPitchMS<<1 // FIND_PITCH_LPC_WIN_MS + pitchEstLPCOrder = 16 + pitchEstComplexity = 2 // 0..2, highest = best + pitchSearchThreshold = 0.3 // first-stage candidate threshold +) + +// findPitchLags whitens the signal and runs the pitch estimator (a port of +// silk_find_pitch_lags_FLP). analysisBuf holds the LTP-memory history followed +// by the current frame (length ltp_mem_length + frame_length). It returns +// whether the frame is voiced, the per-subframe pitch lags, the lag/contour +// indices, and the whitening residual (reused by the LTP analysis). +func (e *Encoder) findPitchLags( + analysisBuf []float32, fsKHz, nbSubfr, speechActivityQ8 int, inputTiltQ15 int32, +) (voiced bool, pitchL []int, lagIndex int16, contourIndex int8, res []float32, predGain float32) { + bufLen := len(analysisBuf) + laPitch := laPitchMS * fsKHz + winLength := findPitchLPCWinMS * fsKHz + + // Windowed signal: rising edge, flat middle, falling edge. + wsig := make([]float32, winLength) + start := bufLen - winLength + applySineWindowFLP(wsig, analysisBuf[start:], 1, laPitch) + copy(wsig[laPitch:winLength-laPitch], analysisBuf[start+laPitch:start+winLength-laPitch]) + applySineWindowFLP(wsig[winLength-laPitch:], analysisBuf[start+winLength-laPitch:], 2, laPitch) + + // Whitening LPC via autocorrelation + Schur. + autoCorr := make([]float32, pitchEstLPCOrder+1) + autocorrelationFLP(autoCorr, wsig, winLength, pitchEstLPCOrder+1) + autoCorr[0] += autoCorr[0]*findPitchWhiteNoiseFraction + 1 + + refl := make([]float32, pitchEstLPCOrder) + resNrg := schurFLP(refl, autoCorr, pitchEstLPCOrder) + predGain = autoCorr[0] / max(resNrg, 1.0) + a := make([]float32, pitchEstLPCOrder) + k2aFLP(a, refl, pitchEstLPCOrder) + bwexpanderFLP(a, pitchEstLPCOrder, findPitchBandwidthExpansion) + + res = make([]float32, bufLen) + lpcAnalysisFilterFLP(res, a, analysisBuf, bufLen, pitchEstLPCOrder) + + // Voicing threshold (search_thres2). + prevVoiced := float32(0) + if e.isPreviousFrameVoiced { + prevVoiced = 1 + } + thrhld := float32(0.6) - + 0.004*pitchEstLPCOrder - + 0.1*float32(speechActivityQ8)*(1.0/256.0) - + 0.15*prevVoiced - + 0.1*float32(inputTiltQ15)*(1.0/32768.0) + + pitchL = make([]int, nbSubfr) + prevLag := 0 + if e.isPreviousFrameVoiced { + prevLag = e.previousLag + } + lagIndex, contourIndex, voiced = pitchAnalysisCore( + res, pitchL, &e.ltpCorr, prevLag, pitchSearchThreshold, thrhld, fsKHz, pitchEstComplexity, nbSubfr) + + return voiced, pitchL, lagIndex, contourIndex, res, predGain +} diff --git a/internal/silk/find_pitch_lags_test.go b/internal/silk/find_pitch_lags_test.go new file mode 100644 index 0000000..fb66ccd --- /dev/null +++ b/internal/silk/find_pitch_lags_test.go @@ -0,0 +1,72 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// TestFindPitchLagsVoiced checks that a clearly periodic signal is detected as +// voiced with a pitch lag near its true period. +func TestFindPitchLagsVoiced(t *testing.T) { + const ( + fsKHz = 16 + nbSubfr = 4 + ltpMemLength = 20 * fsKHz + frameLength = 20 * fsKHz + period = 96 // 16 kHz / 96 ~ 167 Hz, inside the pitch range + ) + + // Impulse train: a strong pulse every `period` samples. + buf := make([]float32, ltpMemLength+frameLength) + for i := range buf { + if i%period == 0 { + buf[i] = 8000 + } + } + + enc := NewEncoder() + voiced, pitchL, lagIndex, contourIndex, res, predGain := enc.findPitchLags(buf, fsKHz, nbSubfr, 200, 0) + + require.True(t, voiced, "periodic signal should be voiced") + require.Len(t, res, ltpMemLength+frameLength) + primaryLag := int(lagIndex) + peMinLagMS*fsKHz + assert.InDeltaf(t, period, primaryLag, 8, "primary lag") + for _, p := range pitchL { + assert.InDeltaf(t, period, p, 12, "subframe lag") + } + // contourIndex selects a row of the stage-3 lag codebook (peNBCbksStage3Max + // entries at complexity 2, the level findPitchLags always requests). + assert.GreaterOrEqualf(t, contourIndex, int8(0), "contour index") + assert.Lessf(t, contourIndex, int8(peNBCbksStage3Max), "contour index") + // predGain is the whitening filter's prediction gain (autoCorr[0]/resNrg). + // This is a short-term (16-tap) LPC gain, not the long-term pitch gain — + // an impulse train with a 96-sample period has no structure a 16-tap + // predictor can see, so the meaningful invariant here is just positivity + // (autoCorr[0] and the clamped resNrg are both always > 0). + assert.Greaterf(t, predGain, float32(0), "prediction gain") +} + +// TestFindPitchLagsUnvoiced checks noise is not classified as strongly voiced. +func TestFindPitchLagsUnvoiced(t *testing.T) { + const ( + fsKHz = 16 + ltpMemLength = 20 * fsKHz + frameLength = 20 * fsKHz + ) + buf := make([]float32, ltpMemLength+frameLength) + state := uint32(1) + for i := range buf { + state = 1664525*state + 1013904223 + buf[i] = float32(int32(state>>16)%1000 - 500) + } + + enc := NewEncoder() + // Not asserting hard (noise can occasionally correlate) — just that it runs. + _, _, _, _, res, _ := enc.findPitchLags(buf, fsKHz, 4, 50, 0) //nolint:dogsled // only the residual is under test + require.Len(t, res, ltpMemLength+frameLength) +} diff --git a/internal/silk/find_pred_coefs.go b/internal/silk/find_pred_coefs.go new file mode 100644 index 0000000..7f74ff7 --- /dev/null +++ b/internal/silk/find_pred_coefs.go @@ -0,0 +1,122 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +const ( + maxPredictionPowerGain = 1e4 + maxPredictionPowerGainAfterReset = 1e2 +) + +// ltpAnalysisFilterFLP forms the LTP residual (silk_LTP_analysis_filter_FLP). +// x is the speech buffer and xBase the index of the first sample to process, +// which must have at least predictLPCOrder + max(pitchL) samples of history +// before it. The output holds nbSubfr blocks of preLength+subfrLength samples, +// each scaled by the subframe inverse gain. +func ltpAnalysisFilterFLP( + ltpRes, x []float32, xBase int, b []float32, pitchL []int, invGains []float32, + subfrLength, nbSubfr, preLength int, +) { + xPtr := xBase + resIdx := 0 + for k := range nbSubfr { + lagPtr := xPtr - pitchL[k] //nolint:gosec // G602: pitchL has nbSubfr entries by caller contract. + invGain := invGains[k] //nolint:gosec // G602: invGains has nbSubfr entries by caller contract. + for i := range subfrLength + preLength { + v := x[xPtr+i] + for j := range ltpOrder { + v -= b[k*ltpOrder+j] * x[lagPtr+ltpOrder/2-j+i] + } + ltpRes[resIdx+i] = v * invGain + } + resIdx += subfrLength + preLength + xPtr += subfrLength + } +} + +// residualEnergyFLP measures the per-subframe LPC residual energy of the +// gain-normalized signal, rescaled by the squared gains +// (silk_residual_energy_FLP). lpcInPre holds nbSubfr blocks of order+subfrLength +// samples; a0 and a1 are the quantized LPC for the first and second frame +// halves (equal when NLSF interpolation is inactive). +func residualEnergyFLP(nrgs, lpcInPre, a0, a1, gains []float32, subfrLength, nbSubfr, order int) { + shift := order + subfrLength + lpcRes := make([]float32, 2*shift) + halves := [][]float32{a0, a1} + for half := 0; half < nbSubfr; half += 2 { + lpcAnalysisFilterFLP(lpcRes, halves[half/2], lpcInPre[half*shift:], 2*shift, order) + nrgs[half] = gains[half] * gains[half] * float32(energyFLP(lpcRes[order:], subfrLength)) + //nolint:gosec // G602: gains/nrgs have nbSubfr entries by caller contract. + nrgs[half+1] = gains[half+1] * gains[half+1] * float32(energyFLP(lpcRes[order+shift:], subfrLength)) + } +} + +// interpolateNLSF linearly interpolates two NLSF vectors (silk_interpolate). +// ifactQ2 (0..4) is the weight on x1. +func interpolateNLSF(xi, x0, x1 []int16, ifactQ2, d int) { + for i := range d { + //nolint:gosec // G115: NLSFs and their interpolation stay within int16. + xi[i] = int16(int32(x0[i]) + (smulbb(int32(x1[i])-int32(x0[i]), int32(ifactQ2)) >> 2)) + } +} + +// findLPCNLSF runs Burg over LPC_in_pre and, for 20 ms frames, searches the +// NLSF interpolation factor with the lowest first-half residual energy +// (silk_find_LPC_FLP). It returns the interpolation index (4 = no +// interpolation) and the NLSFs to quantize. +func (e *Encoder) findLPCNLSF( + lpcInPre []float32, minInvGain float32, bandwidth Bandwidth, order, nbSubfr, subfrLength int, +) (int, []int16) { + blockLen := subfrLength + order + interpQ2 := 4 + aFull := make([]float32, order) + resNrg := burgModifiedFLP(aFull, lpcInPre, minInvGain, blockLen, nbSubfr, order) + + nlsf := make([]int16, order) + if !e.firstFrameAfterReset && nbSubfr == maxSubframeCount { + half := maxSubframeCount / 2 + aTmp := make([]float32, order) + resNrg -= burgModifiedFLP(aTmp, lpcInPre[half*blockLen:], minInvGain, blockLen, half, order) + a2nlsfFLP(nlsf, aTmp, order) // NLSFs of the last 10 ms + + resNrg2nd := float32(math.MaxFloat32) + nlsf0 := make([]int16, order) + lpcRes := make([]float32, 2*blockLen) + aInterp := make([]float32, order) + for k := 3; k >= 0; k-- { + interpolateNLSF(nlsf0, e.prevNLSFq, nlsf, k, order) + aQ12 := nlsfToLPCQ12(nlsf0, bandwidth) + for i := range order { + aInterp[i] = float32(aQ12[i]) * (1.0 / 4096.0) + } + lpcAnalysisFilterFLP(lpcRes, aInterp, lpcInPre, 2*blockLen, order) + resInterp := float32(energyFLP(lpcRes[order:], subfrLength)) + + float32(energyFLP(lpcRes[order+blockLen:], subfrLength)) + if resInterp < resNrg { + resNrg = resInterp + interpQ2 = k + } else if resInterp > resNrg2nd { + break + } + resNrg2nd = resInterp + } + } + if interpQ2 == 4 { + a2nlsfFLP(nlsf, aFull, order) + } + + return interpQ2, nlsf +} + +// predCoefsMinInvGain returns the minimum inverse prediction gain used to cap +// the short-term predictor (silk_find_pred_coefs_FLP). +func predCoefsMinInvGain(firstFrameAfterReset bool, ltpredCodGain, codingQuality float32) float32 { + if firstFrameAfterReset { + return 1.0 / maxPredictionPowerGainAfterReset + } + minInvGain := float32(math.Pow(2, float64(ltpredCodGain)/3)) / maxPredictionPowerGain + + return minInvGain / (0.25 + 0.75*codingQuality) +} diff --git a/internal/silk/find_pred_coefs_test.go b/internal/silk/find_pred_coefs_test.go new file mode 100644 index 0000000..05b226b --- /dev/null +++ b/internal/silk/find_pred_coefs_test.go @@ -0,0 +1,143 @@ +// 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 TestLtpAnalysisFilterFLPZeroCoefsPassesThrough(t *testing.T) { + const ( + nbSubfr = 2 + subfrLength = 8 + preLength = 2 + lag = 10 + // x must have at least lag+ltpOrder/2 samples of history before + // xBase, since the filter reads back to lagPtr-ltpOrder/2. + xBase = lag + ltpOrder/2 + ) + signal := make([]float32, xBase+nbSubfr*subfrLength+preLength) + for i := range signal { + signal[i] = float32(i + 1) + } + b := make([]float32, ltpOrder*nbSubfr) // zero LTP taps: residual == input. + pitchL := []int{lag, lag} + invGains := []float32{1, 2} + + ltpRes := make([]float32, nbSubfr*(subfrLength+preLength)) + ltpAnalysisFilterFLP(ltpRes, signal, xBase, b, pitchL, invGains, subfrLength, nbSubfr, preLength) + + for k := range nbSubfr { + for i := range subfrLength + preLength { + want := signal[xBase+k*subfrLength+i] * invGains[k] + assert.InDeltaf(t, want, ltpRes[k*(subfrLength+preLength)+i], 1e-4, "subfr %d sample %d", k, i) + } + } +} + +func TestResidualEnergyFLPZeroCoefsMatchesRawEnergy(t *testing.T) { + const ( + order = 4 + subfrLength = 6 + nbSubfr = 4 + shift = order + subfrLength + ) + lpcInPre := make([]float32, nbSubfr*shift) + state := uint32(11) + for i := range lpcInPre { + state = 1664525*state + 1013904223 + lpcInPre[i] = float32(int32(state>>16)%2000-1000) / 100 + } + a0 := make([]float32, order) // zero AR coefs: residual == input, unfiltered. + a1 := make([]float32, order) + gains := []float32{1, 1, 1, 1} + + nrgs := make([]float32, nbSubfr) + residualEnergyFLP(nrgs, lpcInPre, a0, a1, gains, subfrLength, nbSubfr, order) + + for half := 0; half < nbSubfr; half += 2 { + base := half * shift + want0 := float32(energyFLP(lpcInPre[base+order:], subfrLength)) + want1 := float32(energyFLP(lpcInPre[base+order+shift:], subfrLength)) + assert.InDeltaf(t, want0, nrgs[half], 1e-2, "nrgs[%d]", half) + //nolint:gosec // G602: half> 1)) + assert.Equal(t, want, xi[i]) + } +} + +func TestPredCoefsMinInvGain(t *testing.T) { + assert.InDelta(t, 1.0/maxPredictionPowerGainAfterReset, predCoefsMinInvGain(true, 5, 0.5), 1e-9) + + got := predCoefsMinInvGain(false, 6.0, 0.5) + want := float32(math.Pow(2, 6.0/3)) / maxPredictionPowerGain / (0.25 + 0.75*0.5) + assert.InDelta(t, want, got, 1e-6) +} + +// TestFindLPCNLSF exercises both the interpolation-search path (4 subframes, +// not the first frame after reset — matching find_pitch_lags_FLP's caller +// convention) and the no-interpolation path (2 subframes, always skips the +// search since it requires MAX_NB_SUBFR), including nlsfToLPCQ12 (only +// reachable from the interpolation search). +func TestFindLPCNLSF(t *testing.T) { + const ( + order = 16 + nbSubfr = 4 + subfrLength = 80 + blockLen = subfrLength + order + ) + lpcInPre := make([]float32, nbSubfr*blockLen) + for i := range lpcInPre { + lpcInPre[i] = float32(3000 * math.Sin(2*math.Pi*float64(i)/37)) + } + + enc := NewEncoder() + enc.firstFrameAfterReset = false + enc.prevNLSFq = genNLSF(order, 42) + + interpQ2, nlsf := enc.findLPCNLSF(lpcInPre, 1e-4, BandwidthWideband, order, nbSubfr, subfrLength) + + require.Len(t, nlsf, order) + require.GreaterOrEqual(t, interpQ2, 0) + require.LessOrEqual(t, interpQ2, 4) + for k := range nlsf { + assert.GreaterOrEqualf(t, nlsf[k], int16(0), "coefficient %d", k) + } + for k := 1; k < order; k++ { + assert.Greaterf(t, nlsf[k], nlsf[k-1], "not increasing at %d", k) + } + + // nbSubfr==2 always takes the no-interpolation path (interpQ2 stays 4) + // regardless of firstFrameAfterReset/prevNLSFq, since the search requires + // exactly MAX_NB_SUBFR subframes. + interpQ2NB, nlsfNB := enc.findLPCNLSF(lpcInPre[:2*blockLen], 1e-4, BandwidthNarrowband, order, 2, subfrLength) + assert.Equal(t, 4, interpQ2NB) + require.Len(t, nlsfNB, order) +} diff --git a/internal/silk/gains.go b/internal/silk/gains.go index b9a1625..9bdff16 100644 --- a/internal/silk/gains.go +++ b/internal/silk/gains.go @@ -15,10 +15,22 @@ const ( gainMaxLogQ7 = 3967 // 31 in Q7, the upper clamp used by silk_log2lin() ) -// encodeSubframeGains (the *Encoder wrapper that also range-encodes the -// indices via emitGainIndices) is deferred to the orchestration piece — same -// reasoning as findPitchLags/encodeNLSF/quantLTPGains, since it needs -// e.haveEncoded and e.rangeEncoder, and that type doesn't exist here yet. +// encodeSubframeGains quantizes the per-subframe target gains (Q16), emits the +// indices, and returns the dequantized gains. It is the counterpart of +// Decoder.decodeSubframeQuantizations and produces bit-identical gains. +func (e *Encoder) encodeSubframeGains( + gainsTargetQ16 []int32, + signalType frameSignalType, + subframeCount int, + isFirstSilkFrameInOpusFrame bool, +) (gainQ16 []float32) { + conditional := !isFirstSilkFrameInOpusFrame && e.haveEncoded + indices, gainQ16, _ := quantizeGains(gainsTargetQ16, &e.previousLogGain, subframeCount, conditional) + e.emitGainIndices(indices, signalType, conditional) + e.haveEncoded = true + + return gainQ16 +} // quantizeGains is silk_gains_quant: it turns target gains into transmit // indices and the dequantized gains, updating the running previousLogGain @@ -79,5 +91,23 @@ func quantizeGains( return indices, gainQ16, gainQ16Int } -// emitGainIndices (the range-coder wiring for quantizeGains's output) is -// deferred alongside encodeSubframeGains, for the same *Encoder reason. +// emitGainIndices range-encodes the gain indices produced by quantizeGains. +func (e *Encoder) emitGainIndices(indices []int8, signalType frameSignalType, conditional bool) { + for subframeIndex, index := range indices { + if subframeIndex == 0 && !conditional { + msb := uint32(index >> 3) //nolint:gosec // G115: index is in [0,63]. + lsb := uint32(index & 0x7) //nolint:gosec // G115 + switch signalType { + case frameSignalTypeInactive: + e.rangeEncoder.EncodeSymbolWithICDF(icdfIndependentQuantizationGainMSBInactive, msb) + case frameSignalTypeVoiced: + e.rangeEncoder.EncodeSymbolWithICDF(icdfIndependentQuantizationGainMSBVoiced, msb) + case frameSignalTypeUnvoiced: + e.rangeEncoder.EncodeSymbolWithICDF(icdfIndependentQuantizationGainMSBUnvoiced, msb) + } + e.rangeEncoder.EncodeSymbolWithICDF(icdfIndependentQuantizationGainLSB, lsb) + } else { + e.rangeEncoder.EncodeSymbolWithICDF(icdfDeltaQuantizationGain, uint32(index)) //nolint:gosec // G115 + } + } +} diff --git a/internal/silk/gains_test.go b/internal/silk/gains_test.go index e408eef..3d71516 100644 --- a/internal/silk/gains_test.go +++ b/internal/silk/gains_test.go @@ -123,6 +123,102 @@ func TestQuantizeGainsRoundTrip(t *testing.T) { } } +// TestEncodeSubframeGainsRoundTrip checks that encoded gains decode back to +// the same values through the decoder, with the range coder and gain state +// in sync — the same check nlsf_encode_test.go's TestQuantizeNLSFRoundTrip +// does for NLSFs, but this one goes through the real range coder since +// *Encoder exists now. +func TestEncodeSubframeGainsRoundTrip(t *testing.T) { + cases := []struct { + name string + signalType frameSignalType + subframeCount int + isFirst bool + haveState bool + prevLogGain int32 + targetsQ16 []int32 + }{ + { + name: "independent_first_voiced", + signalType: frameSignalTypeVoiced, + subframeCount: 4, + isFirst: true, + haveState: false, + prevLogGain: 10, + targetsQ16: []int32{200000, 500000, 1500000, 800000}, + }, + { + name: "independent_first_unvoiced", + signalType: frameSignalTypeUnvoiced, + subframeCount: 4, + isFirst: true, + haveState: false, + prevLogGain: 10, + targetsQ16: []int32{81920, 120000, 90000, 300000}, + }, + { + name: "independent_first_inactive", + signalType: frameSignalTypeInactive, + subframeCount: 2, + isFirst: true, + haveState: false, + prevLogGain: 10, + targetsQ16: []int32{400000, 250000}, + }, + { + name: "conditional_first_delta", + signalType: frameSignalTypeVoiced, + subframeCount: 4, + isFirst: false, + haveState: true, + prevLogGain: 30, + targetsQ16: []int32{600000, 650000, 500000, 700000}, + }, + { + name: "large_upward_jump_double_step", + signalType: frameSignalTypeVoiced, + subframeCount: 4, + isFirst: true, + haveState: false, + prevLogGain: 10, + targetsQ16: []int32{81920, 1000000000, 1500000000, 90000}, + }, + { + name: "monotone_decrease", + signalType: frameSignalTypeUnvoiced, + subframeCount: 4, + isFirst: false, + haveState: true, + prevLogGain: 45, + targetsQ16: []int32{900000, 400000, 150000, 82000}, + }, + } + + for _, tc := range cases { + t.Run(tc.name, func(t *testing.T) { + enc := NewEncoder() + enc.haveEncoded = tc.haveState + enc.previousLogGain = tc.prevLogGain + enc.rangeEncoder.Init() + + gains := enc.encodeSubframeGains(tc.targetsQ16, tc.signalType, tc.subframeCount, tc.isFirst) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.haveDecoded = tc.haveState + dec.previousLogGain = tc.prevLogGain + dec.rangeDecoder.Init(data) + + decGains := dec.decodeSubframeQuantizations(tc.signalType, tc.subframeCount, tc.isFirst) + + require.Equal(t, gains, decGains, "reconstructed gains differ") + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange(), "range coder desync") + assert.Equal(t, enc.previousLogGain, dec.previousLogGain, "previousLogGain desync") + }) + } +} + // TestLin2LogRoundTripAgainstGainDequant checks that lin2log is a near-inverse // of the silk_log2lin() spelled out inline in the gain code, matching the // reference tolerance (they are only approximate inverses). diff --git a/internal/silk/ltp_quant.go b/internal/silk/ltp_quant.go index d9e8dee..3fc5e0f 100644 --- a/internal/silk/ltp_quant.go +++ b/internal/silk/ltp_quant.go @@ -10,7 +10,11 @@ import "math" // for the periodicity index and per-subframe filter indices that minimize a // weighted quantization error plus rate cost. -const nLTPCodebooks = 3 +const ( + maxSumLogGainDBQ7 = 5333 // round(MAX_SUM_LOG_GAIN_DB/6 * 128), MAX_SUM_LOG_GAIN_DB=250 + ltpGainSafetyQ7 = 51 // round(0.4 * 128) + nLTPCodebooks = 3 +) //nolint:gochecknoglobals // LTP codebook effective gains and bit costs (tables_LTP.c). var ( @@ -136,8 +140,71 @@ func vqWMatEC( 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. +// quantLTPGains quantizes the LTP taps for all subframes, returning the Q14 +// filter coefficients, per-subframe codebook indices, the periodicity index, +// and the LTP prediction gain in dB. sumLogGainQ7 carries the cumulative gain +// limit across frames. +func (e *Encoder) quantLTPGains( + xx, xX []float32, subfrLen, nbSubfr int, +) (ltpCoefQ14 []int16, cbkIndex []int8, periodicityIndex int, predGainDB float32) { + xxQ17 := make([]int32, nbSubfr*ltpMatrixSize) + xXQ17 := make([]int32, nbSubfr*ltpOrder) + for i := range xxQ17 { + xxQ17[i] = int32(math.RoundToEven(float64(xx[i]) * 131072.0)) + } + for i := range xXQ17 { + xXQ17[i] = int32(math.RoundToEven(float64(xX[i]) * 131072.0)) + } + + ltpCoefQ14 = make([]int16, nbSubfr*ltpOrder) + cbkIndex = make([]int8, nbSubfr) + tempIdx := make([]int8, nbSubfr) + + minRateDist := int32(math.MaxInt32) + bestSumLogGain := int32(0) + var lastResNrg int32 + for k := range nLTPCodebooks { //nolint:varnamelen // k indexes the codebook, as in the C reference. + cb := ltpCodebook(k) + cbGain := ltpGainTable(k) + clQ5 := ltpBitsTable(k) + size := ltpVQSizes[k] + + resNrg := int32(0) + rateDist := int32(0) + sumLogGainTmp := e.sumLogGainQ7 + for j := range nbSubfr { + maxGainQ7 := log2lin((maxSumLogGainDBQ7-sumLogGainTmp)+(7<<7)) - ltpGainSafetyQ7 + ind, resNrgSubfr, rateDistSubfr, gainQ7 := vqWMatEC( + xxQ17[j*ltpMatrixSize:], xXQ17[j*ltpOrder:], cb, cbGain, clQ5, subfrLen, maxGainQ7, size) + tempIdx[j] = int8(ind) //nolint:gosec // G115 + resNrg = addPosSat32(resNrg, resNrgSubfr) + rateDist = addPosSat32(rateDist, rateDistSubfr) + sumLogGainTmp = max(0, sumLogGainTmp+lin2log(ltpGainSafetyQ7+gainQ7)-(7<<7)) + } + if rateDist <= minRateDist { + minRateDist = rateDist + periodicityIndex = k + copy(cbkIndex, tempIdx) + bestSumLogGain = sumLogGainTmp + } + lastResNrg = resNrg + } + + cb := ltpCodebook(periodicityIndex) + for j := range nbSubfr { + for k := range ltpOrder { + //nolint:gosec // G115: int8 codebook value, well within int16 range after <<7. + ltpCoefQ14[j*ltpOrder+k] = int16(int32(cb[cbkIndex[j]][k]) << 7) + } + } + + if nbSubfr == 2 { + lastResNrg >>= 1 + } else { + lastResNrg >>= 2 + } + e.sumLogGainQ7 = bestSumLogGain + predGainDB = float32(smulbb(-3, lin2log(lastResNrg)-(15<<7))) / 128.0 + + return ltpCoefQ14, cbkIndex, periodicityIndex, predGainDB +} diff --git a/internal/silk/ltp_quant_test.go b/internal/silk/ltp_quant_test.go index 84358e7..3cdcd76 100644 --- a/internal/silk/ltp_quant_test.go +++ b/internal/silk/ltp_quant_test.go @@ -11,6 +11,45 @@ import ( "github.com/stretchr/testify/require" ) +// TestQuantLTPGains feeds quantLTPGains the correlation matrices findLTPFLP +// actually produces (the same wiring find_pred_coefs_FLP.c uses: find_LTP's +// output goes straight into quant_LTP_gains), across multiple calls to also +// exercise the cumulative sumLogGainQ7 state threading across frames. +func TestQuantLTPGains(t *testing.T) { + const ( + nbSubfr = 4 + subfrLength = 80 + maxLag = 120 + ) + total := maxLag + ltpOrder + nbSubfr*subfrLength + ltpOrder + r := make([]float32, total) + for i := range r { + r[i] = float32(500 * math.Sin(2*math.Pi*float64(i)/41)) + } + 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) + + enc := NewEncoder() + for i := range 3 { // repeat: also checks sumLogGainQ7 keeps accumulating sanely. + ltpCoefQ14, cbkIndex, periodicityIndex, predGainDB := enc.quantLTPGains(xx, xX, subfrLength, nbSubfr) + + require.Len(t, ltpCoefQ14, nbSubfr*ltpOrder) + require.Len(t, cbkIndex, nbSubfr) + require.GreaterOrEqualf(t, periodicityIndex, 0, "call %d", i) + require.Lessf(t, periodicityIndex, nLTPCodebooks, "call %d", i) + for k, idx := range cbkIndex { + assert.GreaterOrEqualf(t, idx, int8(0), "call %d cbkIndex[%d]", i, k) + } + assert.Falsef(t, math.IsNaN(float64(predGainDB)) || math.IsInf(float64(predGainDB), 0), "call %d predGainDB", i) + assert.GreaterOrEqualf(t, enc.sumLogGainQ7, int32(0), "call %d sumLogGainQ7", i) + assert.LessOrEqualf(t, enc.sumLogGainQ7, int32(maxSumLogGainDBQ7), "call %d sumLogGainQ7", i) + } +} + // 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) { diff --git a/internal/silk/nlsf_encode.go b/internal/silk/nlsf_encode.go index 27048d8..e4e6342 100644 --- a/internal/silk/nlsf_encode.go +++ b/internal/silk/nlsf_encode.go @@ -9,10 +9,26 @@ import ( ) const ( + nlsfQuantMaxAmplitude = 4 nlsfQuantMaxAmplitudeExt = 10 nlsfLevelAdjQ10 = 102 // round(NLSF_QUANT_LEVEL_ADJ * 1024) ) +// nlsfToLPCQ12 converts NLSFs to Q12 LPC coefficients, limited and gain +// checked as the decoder would (silk_NLSF2A + the limiting steps applied when +// decoding a frame). Used by findLPCNLSF's interpolation search. +func nlsfToLPCQ12(nlsfQ15 []int16, bandwidth Bandwidth) []int16 { + d := NewDecoder() + a32Q17 := d.convertNormalizedLSFsToLPCCoefficients(nlsfQ15, bandwidth) + d.limitLPCCoefficientsRange(a32Q17) + d.limitLPCFilterPredictionGainInto(a32Q17, 0) + + out := make([]int16, len(d.aQ12Int[0])) + copy(out, d.aQ12Int[0]) + + return out +} + // nlsfStepSizes returns the Q16 quantization step and its Q6 inverse. func nlsfStepSizes(bandwidth Bandwidth) (qstepQ16, invQstepQ6 int32) { if bandwidth == BandwidthWideband { @@ -51,6 +67,18 @@ func nlsfSecondOperand(ind, qstepQ16 int32) int32 { return (((ind << 10) - int32(sign(int(ind)))*nlsfLevelAdjQ10) * qstepQ16) >> 16 //nolint:gosec // G115 } +// encodeNLSF quantizes and range-encodes the input NLSF vector, returning the +// quantized NLSFs the decoder will reconstruct. It searches every stage-1 +// codebook vector and greedily quantizes the stage-2 residual for each, +// keeping the lowest weighted distortion. +func (e *Encoder) encodeNLSF(nlsfQ15 []int16, bandwidth Bandwidth, voiced bool) []int16 { + stabilizeNLSF(nlsfQ15, len(nlsfQ15), bandwidth) + index1, indices2, quantized := quantizeNLSF(nlsfQ15, bandwidth) + e.emitNLSFIndices(index1, indices2, bandwidth, voiced) + + return quantized +} + // quantizeNLSF searches the two-stage NLSF codebooks (silk_NLSF_encode) and // returns the stage-1 index, stage-2 indices, and reconstructed NLSF vector. // The input must already be stabilized. @@ -126,6 +154,46 @@ func quantizeNLSF(nlsfQ15 []int16, bandwidth Bandwidth) (int, []int8, []int16) { return bestIndex1, bestIndices2, bestNLSF } +// emitNLSFIndices range-encodes the NLSF codebook indices. +func (e *Encoder) emitNLSFIndices(index1 int, indices2 []int8, bandwidth Bandwidth, voiced bool) { + cb2Select := codebookNormalizedLSFStageTwoIndexNarrowbandOrMediumband + stageOnePDF := icdfNormalizedLSFStageOneIndexNarrowbandOrMediumbandUnvoiced + if bandwidth == BandwidthWideband { + cb2Select = codebookNormalizedLSFStageTwoIndexWideband + stageOnePDF = icdfNormalizedLSFStageOneIndexWidebandUnvoiced + if voiced { + stageOnePDF = icdfNormalizedLSFStageOneIndexWidebandVoiced + } + } else if voiced { + stageOnePDF = icdfNormalizedLSFStageOneIndexNarrowbandOrMediumbandVoiced + } + + e.rangeEncoder.EncodeSymbolWithICDF(stageOnePDF, uint32(index1)) //nolint:gosec // G115 + cb2 := cb2Select[index1] + for k := range indices2 { //nolint:varnamelen // k indexes the coefficient. + v := int(indices2[k]) //nolint:varnamelen // v is the stage-2 index value. + switch { + case v <= -nlsfQuantMaxAmplitude: + e.rangeEncoder.EncodeSymbolWithICDF(icdfNormalizedLSFStageTwoIndex[cb2[k]], 0) + e.rangeEncoder.EncodeSymbolWithICDF( + icdfNormalizedLSFStageTwoIndexExtension, + uint32(-nlsfQuantMaxAmplitude-v), //nolint:gosec // G115 + ) + case v >= nlsfQuantMaxAmplitude: + e.rangeEncoder.EncodeSymbolWithICDF(icdfNormalizedLSFStageTwoIndex[cb2[k]], 2*nlsfQuantMaxAmplitude) + e.rangeEncoder.EncodeSymbolWithICDF( + icdfNormalizedLSFStageTwoIndexExtension, + uint32(v-nlsfQuantMaxAmplitude), //nolint:gosec // G115 + ) + default: + e.rangeEncoder.EncodeSymbolWithICDF( + icdfNormalizedLSFStageTwoIndex[cb2[k]], + uint32(v+nlsfQuantMaxAmplitude), //nolint:gosec // G115 + ) + } + } +} + // stabilizeNLSF enforces the minimum spacing between consecutive NLSF // coefficients (RFC 6716 Section 4.2.7.5.4). // diff --git a/internal/silk/nlsf_encode_test.go b/internal/silk/nlsf_encode_test.go index 6d8954c..009bfdb 100644 --- a/internal/silk/nlsf_encode_test.go +++ b/internal/silk/nlsf_encode_test.go @@ -94,6 +94,64 @@ func reconstructNLSF(index1 int, indices2 []int8, bandwidth Bandwidth) []int16 { return nlsf } +// TestEncodeNLSFRoundTrip encodes an NLSF vector through the real range coder +// and decodes it back, checking the reconstructed vector and range-coder +// state match — now possible since *Encoder exists, unlike +// TestQuantizeNLSFRoundTrip above (which predates it and stays as a +// non-bitstream cross-check of the same math). +func TestEncodeNLSFRoundTrip(t *testing.T) { + cases := []struct { + bandwidth Bandwidth + order int + }{ + {BandwidthNarrowband, 10}, + {BandwidthMediumband, 10}, + {BandwidthWideband, 16}, + } + + for _, tc := range cases { + for _, voiced := range []bool{false, true} { + for seed := range 12 { + name := fmt.Sprintf("bw%d_voiced%t_seed%d", tc.bandwidth, voiced, seed) + t.Run(name, func(t *testing.T) { + input := genNLSF(tc.order, uint32(seed*97+tc.order+boolSeed(voiced))) //nolint:gosec // G115 + + enc := NewEncoder() + enc.rangeEncoder.Init() + quant := enc.encodeNLSF(append([]int16(nil), input...), tc.bandwidth, voiced) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.rangeDecoder.Init(data) + index1 := dec.normalizeLineSpectralFrequencyStageOne(voiced, tc.bandwidth) + dLPC, resQ10 := dec.normalizeLineSpectralFrequencyStageTwo(tc.bandwidth, index1) + nlsf := dec.normalizeLineSpectralFrequencyCoefficients(dLPC, tc.bandwidth, resQ10, index1) + dec.normalizeLSFStabilization(nlsf, dLPC, tc.bandwidth) + + require.Len(t, nlsf, len(quant)) + for k := range quant { + require.Equalf(t, quant[k], nlsf[k], "coefficient %d", k) + } + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange(), "range coder desync") + + for k := 1; k < len(quant); k++ { + require.Greaterf(t, quant[k], quant[k-1], "not increasing at %d", k) + } + }) + } + } + } +} + +func boolSeed(b bool) int { + if b { + return 1 + } + + return 0 +} + // genClusteredNLSF builds a tightly clustered NLSF vector: all coefficients sit // close together around a seed-dependent center, so stabilization pushes them to // minimum spacing and the stage-2 residuals grow large. These are the vectors @@ -157,6 +215,29 @@ func TestQuantizeNLSFClusteredStaysInRange(t *testing.T) { reconstructed := reconstructNLSF(index1, indices2, tc.bandwidth) stabilizeNLSF(reconstructed, len(reconstructed), tc.bandwidth) assert.Equal(t, quant, reconstructed) + + // The original form of this regression: the full encode->decode + // round trip through the real range coder must stay in sync. + for _, voiced := range []bool{false, true} { + enc := NewEncoder() + enc.rangeEncoder.Init() + encQuant := enc.encodeNLSF(append([]int16(nil), input...), tc.bandwidth, voiced) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.rangeDecoder.Init(data) + decIndex1 := dec.normalizeLineSpectralFrequencyStageOne(voiced, tc.bandwidth) + decDLPC, decResQ10 := dec.normalizeLineSpectralFrequencyStageTwo(tc.bandwidth, decIndex1) + decNLSF := dec.normalizeLineSpectralFrequencyCoefficients(decDLPC, tc.bandwidth, decResQ10, decIndex1) + dec.normalizeLSFStabilization(decNLSF, decDLPC, tc.bandwidth) + + require.Len(t, decNLSF, len(encQuant)) + for k := range encQuant { + require.Equalf(t, encQuant[k], decNLSF[k], "voiced=%v coefficient %d", voiced, k) + } + assert.Equalf(t, encRange, dec.rangeDecoder.FinalRange(), "voiced=%v range coder desync", voiced) + } }) } } diff --git a/internal/silk/noise_shape_analysis.go b/internal/silk/noise_shape_analysis.go new file mode 100644 index 0000000..8eadb63 --- /dev/null +++ b/internal/silk/noise_shape_analysis.go @@ -0,0 +1,311 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// Faithful port of silk_noise_shape_analysis_FLP and silk_process_gains_FLP +// (warping disabled, matching the low-complexity / non-delayed-decision path). + +const ( + bgSNRDecrDB = 2.0 + harmSNRIncrDB = 2.0 + shapeWhiteNoiseFraction = 3e-5 + bandwidthExpansion = 0.94 + harmonicShapingConst = 0.3 + highRateHarmonicShaping = 0.2 + hpNoiseCoefConst = 0.25 + harmHPNoiseCoefConst = 0.35 + lowFreqShapingConst = 4.0 + lowQualityLFShapingDecr = 0.5 + subframeSmoothCoef = 0.4 + energyVarThresholdQntOff = 0.6 + minQGainDB = 2 + subFrameLengthMS = 5 + + lambdaOffset = 1.2 + lambdaSpeechAct = -0.2 + lambdaDelayedDecisions = -0.05 + lambdaInputQuality = -0.1 + lambdaCodingQuality = -0.2 + lambdaQuantOffsetConst = 0.8 + nStatesDelayedDecision = 1 // non-delayed NSQ + shapeLPCOrderLowComplex = 12 + laShapeMSLowComplex = 3 +) + +// f2iQ rounds x to the nearest integer (silk_float2int). +func f2iQ(x float32) int32 { + return int32(math.RoundToEven(float64(x))) +} + +// limitCoefsFLP bandwidth-expands an AR filter until every coefficient's +// magnitude is within limit (silk_noise_shape_analysis's limit_coefs). +func limitCoefsFLP(coefs []float32, limit float32, order int) { + for iter := range 10 { + maxAbs := float32(-1) + ind := 0 + for i := range order { + if a := absFloat32(coefs[i]); a > maxAbs { + maxAbs = a + ind = i + } + } + if maxAbs <= limit { + return + } + chirp := 0.99 - (0.8+0.1*float32(iter))*(maxAbs-limit)/(maxAbs*float32(ind+1)) + bwexpanderFLP(coefs, order, chirp) + } +} + +func absFloat32(x float32) float32 { + if x < 0 { + return -x + } + + return x +} + +// shapeResult holds the noise-shaping outputs for one frame. +type shapeResult struct { + gains []float32 + arQ13 []int16 + tiltQ14 []int32 + lfShpQ14 []int32 + harmShapeQ14 []int32 + inputQuality float32 + codingQuality float32 + snrAdjDB float32 + quantOffset frameQuantizationOffsetType // unvoiced choice; voiced set in processGains +} + +// noiseShapeAnalysis computes the noise-shaping AR filters, initial gains, +// spectral tilt, low-frequency and harmonic shaping (silk_noise_shape_analysis_FLP). +// shapeBuf holds la_shape samples of history, the frame, then la_shape samples +// of look-ahead padding. pitchRes is the whitened LPC residual for the current +// frame (findPitchLags's res, offset past its ltpMemLength history) — the +// reference uses this, not the raw signal, for the unvoiced sparseness measure +// that picks quantOffsetType; reusing shapeBuf there was a bug. +// +//nolint:cyclop // faithful port of the noise-shape analysis stage. +func (e *Encoder) noiseShapeAnalysis( + shapeBuf, pitchRes []float32, + signalType frameSignalType, + pitchL []int, + predGain float32, + snrDBQ7 int32, + speechActQ8 int, + qualityBands [vadNBands]int32, + fsKHz, nbSubfr, subfrLength int, +) *shapeResult { + order := shapeLPCOrderLowComplex + laShape := laShapeMSLowComplex * fsKHz + shapeWinLength := subFrameLengthMS*fsKHz + 2*laShape + + sr := &shapeResult{ + gains: make([]float32, nbSubfr), + arQ13: make([]int16, nbSubfr*maxShapeLPCOrder), + tiltQ14: make([]int32, nbSubfr), + lfShpQ14: make([]int32, nbSubfr), + harmShapeQ14: make([]int32, nbSubfr), + } + + // Gain control. + snrAdjDB := float32(snrDBQ7) * (1.0 / 128.0) + sr.inputQuality = 0.5 * (float32(qualityBands[0]) + float32(qualityBands[1])) * (1.0 / 32768.0) + sr.codingQuality = sigmoid(0.25 * (snrAdjDB - 20.0)) + speechAct := float32(speechActQ8) * (1.0 / 256.0) + + // useCBR == 0: reduce coding SNR during low speech activity. + b := 1.0 - speechAct + snrAdjDB -= bgSNRDecrDB * sr.codingQuality * (0.5 + 0.5*sr.inputQuality) * b * b + if signalType == frameSignalTypeVoiced { + snrAdjDB += harmSNRIncrDB * e.ltpCorr + } else { + snrAdjDB += (-0.4*float32(snrDBQ7)*(1.0/128.0) + 6.0) * (1.0 - sr.inputQuality) + } + sr.snrAdjDB = snrAdjDB + + // Quantizer offset for unvoiced from a sparseness measure. + sr.quantOffset = frameQuantizationOffsetTypeLow + if signalType == frameSignalTypeVoiced { + sr.quantOffset = frameQuantizationOffsetTypeLow // may be overridden in processGains + } else { + nSamples := 2 * fsKHz + nSegs := subFrameLengthMS * nbSubfr / 2 + var energyVariation, logEnergyPrev float32 + ptr := 0 + for k := range nSegs { + nrg := float32(nSamples) + float32(energyFLP(pitchRes[ptr:], nSamples)) + logEnergy := silkLog2(float64(nrg)) + if k > 0 { + energyVariation += absFloat32(logEnergy - logEnergyPrev) + } + logEnergyPrev = logEnergy + ptr += nSamples + } + if energyVariation > energyVarThresholdQntOff*float32(nSegs-1) { + sr.quantOffset = frameQuantizationOffsetTypeLow + } else { + sr.quantOffset = frameQuantizationOffsetTypeHigh + } + } + + // Bandwidth expansion for the shaping filter. + strength := float32(findPitchWhiteNoiseFraction) * predGain + bwExp := float32(bandwidthExpansion) / (1.0 + strength*strength) + + // Per-subframe shaping AR coefficients and gains. + xWindowed := make([]float32, shapeWinLength) + autoCorr := make([]float32, order+1) + rc := make([]float32, order) + flatPart := fsKHz * 3 + slopePart := (shapeWinLength - flatPart) / 2 + xPtr := 0 + for k := range nbSubfr { //nolint:varnamelen // k is the subframe index throughout. + applySineWindowFLP(xWindowed, shapeBuf[xPtr:], 1, slopePart) + copy(xWindowed[slopePart:slopePart+flatPart], shapeBuf[xPtr+slopePart:]) + applySineWindowFLP(xWindowed[slopePart+flatPart:], shapeBuf[xPtr+slopePart+flatPart:], 2, slopePart) + xPtr += subfrLength + + autocorrelationFLP(autoCorr, xWindowed, shapeWinLength, order+1) + autoCorr[0] += autoCorr[0]*shapeWhiteNoiseFraction + 1.0 + + nrg := schurFLP(rc, autoCorr, order) + arSub := make([]float32, order) + k2aFLP(arSub, rc, order) + sr.gains[k] = float32(math.Sqrt(float64(nrg))) + bwexpanderFLP(arSub, order, bwExp) + limitCoefsFLP(arSub, 3.999, order) + for j := range order { + //nolint:gosec // G115: bandwidth-expanded, limited shaping coefs fit int16. + sr.arQ13[k*maxShapeLPCOrder+j] = int16(f2iQ(arSub[j] * 8192.0)) + } + } + + // Gain tweaking: higher gains during low speech activity. + gainMult := float32(math.Pow(2.0, float64(-0.16*snrAdjDB))) + gainAdd := float32(math.Pow(2.0, 0.16*minQGainDB)) + for k := range nbSubfr { + sr.gains[k] = sr.gains[k]*gainMult + gainAdd + } + + // Low-frequency shaping and spectral tilt. + lfStrength := lowFreqShapingConst * (1.0 + lowQualityLFShapingDecr*(float32(qualityBands[0])*(1.0/32768.0)-1.0)) + lfStrength *= speechAct + var tilt float32 + lfMA := make([]float32, nbSubfr) + lfAR := make([]float32, nbSubfr) + if signalType == frameSignalTypeVoiced { + for k := range nbSubfr { + bb := 0.2/float32(fsKHz) + 3.0/float32(pitchL[k]) + lfMA[k] = -1.0 + bb + lfAR[k] = 1.0 - bb - bb*lfStrength + } + tilt = -hpNoiseCoefConst - (1.0-hpNoiseCoefConst)*harmHPNoiseCoefConst*speechAct + } else { + bb := 1.3 / float32(fsKHz) + lfMA[0] = -1.0 + bb + lfAR[0] = 1.0 - bb - bb*lfStrength*0.6 + for k := 1; k < nbSubfr; k++ { + lfMA[k] = lfMA[0] + lfAR[k] = lfAR[0] + } + tilt = -hpNoiseCoefConst + } + + // Harmonic shaping (voiced). + var harmShapeGain float32 + if signalType == frameSignalTypeVoiced { + harmShapeGain = harmonicShapingConst + harmShapeGain += highRateHarmonicShaping * (1.0 - (1.0-sr.codingQuality)*sr.inputQuality) + harmShapeGain *= float32(math.Sqrt(float64(e.ltpCorr))) + } + + // Smooth over subframes. + for k := range nbSubfr { + e.harmShapeGainSmth += subframeSmoothCoef * (harmShapeGain - e.harmShapeGainSmth) + e.tiltSmth += subframeSmoothCoef * (tilt - e.tiltSmth) + sr.harmShapeQ14[k] = f2iQ(e.harmShapeGainSmth * 16384.0) + sr.tiltQ14[k] = f2iQ(e.tiltSmth * 16384.0) + sr.lfShpQ14[k] = (f2iQ(lfAR[k]*16384.0) << 16) | int32(uint16(int16(f2iQ(lfMA[k]*16384.0)))) //nolint:gosec + } + + return sr +} + +// processGains reduces gains for high LTP gain, soft-limits them against the +// residual energy, quantizes them, and computes Lambda (silk_process_gains_FLP). +func (e *Encoder) processGains( + sr *shapeResult, + resNrg []float32, + signalType frameSignalType, + ltpredCodGain float32, + snrDBQ7 int32, + speechActQ8 int, + inputTiltQ15 int32, + subfrLength, nbSubfr int, + conditional bool, +) (gainsQ16Int []int32, gainIndices []int8, lambdaQ10 int32, quantOffset frameQuantizationOffsetType) { + quantOffset = sr.quantOffset + + // Gain reduction when LTP coding gain is high. + if signalType == frameSignalTypeVoiced { + s := 1.0 - 0.5*sigmoid(0.25*(ltpredCodGain-12.0)) + for k := range nbSubfr { + sr.gains[k] *= s + } + } + + // Soft limit on the ratio of residual energy to squared gains. + invMaxSqrVal := float32(math.Pow(2.0, float64(0.33*(21.0-float32(snrDBQ7)*(1.0/128.0))))) / float32(subfrLength) + gainsTargetQ16 := make([]int32, nbSubfr) + for k := range nbSubfr { + gain := sr.gains[k] + gain = float32(math.Sqrt(float64(gain*gain + resNrg[k]*invMaxSqrVal))) + sr.gains[k] = float32(math.Min(float64(gain), 32767.0)) + gainsTargetQ16[k] = int32(sr.gains[k] * 65536.0) + } + + // Quantize. + gainIndices, gainsFloat, gainsQ16Int := quantizeGains(gainsTargetQ16, &e.previousLogGain, nbSubfr, conditional) + for k := range nbSubfr { + sr.gains[k] = gainsFloat[k] * (1.0 / 65536.0) + } + + // Quantizer offset for voiced. + if signalType == frameSignalTypeVoiced { + if ltpredCodGain+float32(inputTiltQ15)*(1.0/32768.0) > 1.0 { + quantOffset = frameQuantizationOffsetTypeLow + } else { + quantOffset = frameQuantizationOffsetTypeHigh + } + } + + // Rate/distortion tradeoff. + offsetIndex := int(quantOffset) - int(frameQuantizationOffsetTypeLow) + signalIndex := int(signalType) - int(frameSignalTypeInactive) + quantOffsetVal := float32(silkSignQuantOffsetQ10(signalIndex, offsetIndex)) * (1.0 / 1024.0) + lambda := float32(lambdaOffset) + + lambdaDelayedDecisions*nStatesDelayedDecision + + lambdaSpeechAct*float32(speechActQ8)*(1.0/256.0) + + lambdaInputQuality*sr.inputQuality + + lambdaCodingQuality*sr.codingQuality + + lambdaQuantOffsetConst*quantOffsetVal + lambdaQ10 = f2iQ(lambda * 1024.0) + + return gainsQ16Int, gainIndices, lambdaQ10, quantOffset +} + +// silkSignQuantOffsetQ10 returns the Q10 quantization offset for the signal and +// offset type (silk_Quantization_Offsets_Q10), indexed as [signalType>>1][offset]. +func silkSignQuantOffsetQ10(signalIndex, offsetIndex int) int32 { + sig := 0 + if signalIndex == 2 { // voiced + sig = 1 + } + + return quantizationOffsetsQ10[sig][offsetIndex] +} diff --git a/internal/silk/noise_shape_analysis_test.go b/internal/silk/noise_shape_analysis_test.go new file mode 100644 index 0000000..e473120 --- /dev/null +++ b/internal/silk/noise_shape_analysis_test.go @@ -0,0 +1,118 @@ +// 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" +) + +// TestNoiseShapeAnalysisAndProcessGainsVoiced runs the full analysis chain +// (findPitchLags -> noiseShapeAnalysis -> processGains) for a periodic +// (voiced) signal, using findPitchLags's own residual as noiseShapeAnalysis's +// pitchRes — the same wiring encode_frame_FLP.c uses, and a regression check +// for the pitchRes/shapeBuf mixup this test also exists to catch. +func TestNoiseShapeAnalysisAndProcessGainsVoiced(t *testing.T) { + const ( + fsKHz = 16 + nbSubfr = 4 + subfrLength = 5 * fsKHz + laShape = 3 * fsKHz + frameLength = nbSubfr * subfrLength + ltpMemLength = 20 * fsKHz + period = 100 + ) + + analysisBuf := make([]float32, ltpMemLength+frameLength) + for i := range analysisBuf { + analysisBuf[i] = float32(3000 * math.Sin(2*math.Pi*float64(i)/period)) + } + + enc := NewEncoder() + enc.ltpCorr = 0.6 + voiced, pitchL, _, _, res, predGain := enc.findPitchLags(analysisBuf, fsKHz, nbSubfr, 200, 0) + require.True(t, voiced, "periodic signal should be voiced") + + shapeBuf := make([]float32, frameLength+2*laShape) + for i := range shapeBuf { + shapeBuf[i] = float32(3000 * math.Sin(2*math.Pi*float64(i)/period)) + } + pitchRes := res[ltpMemLength:] + qualityBands := [vadNBands]int32{20000, 20000, 15000, 15000} + + sr := enc.noiseShapeAnalysis( + shapeBuf, pitchRes, frameSignalTypeVoiced, pitchL, predGain, 20*128, 200, qualityBands, fsKHz, nbSubfr, subfrLength) + + require.NotNil(t, sr) + require.Len(t, sr.gains, nbSubfr) + for k, g := range sr.gains { + assert.Greaterf(t, g, float32(0), "gain %d", k) + assert.Falsef(t, math.IsNaN(float64(g)), "gain %d", k) + } + require.Len(t, sr.arQ13, nbSubfr*maxShapeLPCOrder) + require.Len(t, sr.tiltQ14, nbSubfr) + require.Len(t, sr.lfShpQ14, nbSubfr) + require.Len(t, sr.harmShapeQ14, nbSubfr) + // Voiced frames leave the sparseness measure untouched — quantOffset + // starts Low and processGains decides the real value below. + assert.Equal(t, frameQuantizationOffsetTypeLow, sr.quantOffset) + + resNrg := make([]float32, nbSubfr) + for k := range resNrg { + resNrg[k] = 500000 + } + + gainsQ16Int, gainIndices, lambdaQ10, quantOffset := enc.processGains( + sr, resNrg, frameSignalTypeVoiced, 0.6, 20*128, 200, 0, subfrLength, nbSubfr, false) + + require.Len(t, gainsQ16Int, nbSubfr) + require.Len(t, gainIndices, nbSubfr) + for k, g := range gainsQ16Int { + assert.Positivef(t, g, "gainsQ16Int %d", k) + } + assert.Positive(t, lambdaQ10) + assert.Contains(t, + []frameQuantizationOffsetType{frameQuantizationOffsetTypeLow, frameQuantizationOffsetTypeHigh}, quantOffset) +} + +// TestNoiseShapeAnalysisUnvoicedSparseness exercises the unvoiced sparseness +// measure directly: a bursty residual (high energy variation between 2ms +// segments) must pick the Low offset, a smooth one the High offset — this is +// exactly the branch that used to read the wrong buffer (shapeBuf instead of +// pitchRes). +func TestNoiseShapeAnalysisUnvoicedSparseness(t *testing.T) { + const ( + fsKHz = 16 + nbSubfr = 4 + subfrLength = 5 * fsKHz + laShape = 3 * fsKHz + frameLength = nbSubfr * subfrLength + ) + qualityBands := [vadNBands]int32{20000, 20000, 15000, 15000} + shapeBuf := make([]float32, frameLength+2*laShape) + pitchL := make([]int, nbSubfr) + + bursty := make([]float32, frameLength) + for i := range bursty { + if (i/(2*fsKHz))%2 == 0 { + bursty[i] = 5000 + } + } + enc := NewEncoder() + srBursty := enc.noiseShapeAnalysis( + shapeBuf, bursty, frameSignalTypeUnvoiced, pitchL, 1.0, 20*128, 200, qualityBands, fsKHz, nbSubfr, subfrLength) + assert.Equal(t, frameQuantizationOffsetTypeLow, srBursty.quantOffset, "bursty residual: expected Low") + + flat := make([]float32, frameLength) + for i := range flat { + flat[i] = float32(500 * math.Sin(2*math.Pi*float64(i)/23)) + } + enc2 := NewEncoder() + srFlat := enc2.noiseShapeAnalysis( + shapeBuf, flat, frameSignalTypeUnvoiced, pitchL, 1.0, 20*128, 200, qualityBands, fsKHz, nbSubfr, subfrLength) + assert.Equal(t, frameQuantizationOffsetTypeHigh, srFlat.quantOffset, "flat residual: expected High") +} diff --git a/internal/silk/pitch_ltp.go b/internal/silk/pitch_ltp.go new file mode 100644 index 0000000..96b48f8 --- /dev/null +++ b/internal/silk/pitch_ltp.go @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// encodePitchLags codes the primary pitch lag and the subframe pitch contour +// index, inverting Decoder.decodePitchLags. It is only called for voiced +// frames. lag is the primary pitch lag; contourIndex selects the contour VQ. +func (e *Encoder) encodePitchLags( + lag int, + contourIndex uint32, + bandwidth Bandwidth, + nanoseconds int, + isFirstSilkFrameInOpusFrame bool, +) { + lagAbsolute := isFirstSilkFrameInOpusFrame || !e.isPreviousFrameVoiced + lowPartICDF, lagScale, lagMin, _ := pitchLagCodebooks(bandwidth) + + switch { + case lagAbsolute: + e.encodeAbsolutePitchLag(lag, lowPartICDF, lagScale, lagMin) + default: + delta := lag - e.previousLag + 9 + if delta >= 1 && delta <= len(icdfPrimaryPitchLagChange)-2 { + e.rangeEncoder.EncodeSymbolWithICDF(icdfPrimaryPitchLagChange, uint32(delta)) //nolint:gosec // G115 + } else { + e.rangeEncoder.EncodeSymbolWithICDF(icdfPrimaryPitchLagChange, 0) + e.encodeAbsolutePitchLag(lag, lowPartICDF, lagScale, lagMin) + } + } + e.previousLag = lag + + _, lagIcdf := pitchContourCodebooks(bandwidth, nanoseconds) + e.rangeEncoder.EncodeSymbolWithICDF(lagIcdf, contourIndex) +} + +// encodeAbsolutePitchLag codes the primary lag as a high and low part. +func (e *Encoder) encodeAbsolutePitchLag(lag int, lowPartICDF []uint, lagScale, lagMin uint32) { + rel := uint32(lag) - lagMin //nolint:gosec // G115 + e.rangeEncoder.EncodeSymbolWithICDF(icdfPrimaryPitchLagHighPart, rel/lagScale) + e.rangeEncoder.EncodeSymbolWithICDF(lowPartICDF, rel%lagScale) +} + +// encodeLTPFilter codes the periodicity index and the per-subframe LTP filter +// index, inverting Decoder.decodeLTPFilterCoefficients. +func (e *Encoder) encodeLTPFilter(periodicityIndex uint32, filterIndices []uint32) { + e.rangeEncoder.EncodeSymbolWithICDF(icdfPeriodicityIndex, periodicityIndex) + filterICDF := ltpFilterIndexICDF(periodicityIndex) + for _, index := range filterIndices { + e.rangeEncoder.EncodeSymbolWithICDF(filterICDF, index) + } +} + +// ltpFilterIndexICDF returns the filter-index PDF for a periodicity index. +func ltpFilterIndexICDF(periodicityIndex uint32) []uint { + switch periodicityIndex { + case 1: + return icdfLTPFilterIndex1 + case 2: + return icdfLTPFilterIndex2 + default: + return icdfLTPFilterIndex0 + } +} + +// encodeLTPScaling codes the LTP scaling index. The caller emits it only for +// the first voiced SILK frame of an Opus frame. +func (e *Encoder) encodeLTPScaling(scaleIndex uint32) { + e.rangeEncoder.EncodeSymbolWithICDF(icdfLTPScalingParameter, scaleIndex) +} + +// encodeLCGSeed codes the excitation seed (RFC 6716 Section 4.2.7.7). +func (e *Encoder) encodeLCGSeed(seed uint32) { + e.rangeEncoder.EncodeSymbolWithICDF(icdfLinearCongruentialGeneratorSeed, seed) +} diff --git a/internal/silk/pitch_ltp_test.go b/internal/silk/pitch_ltp_test.go new file mode 100644 index 0000000..253b1e1 --- /dev/null +++ b/internal/silk/pitch_ltp_test.go @@ -0,0 +1,145 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "fmt" + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +func TestEncodePitchLagsRoundTrip(t *testing.T) { + configs := []struct { + bandwidth Bandwidth + nanoseconds int + }{ + {BandwidthNarrowband, nanoseconds10Ms}, + {BandwidthNarrowband, nanoseconds20Ms}, + {BandwidthMediumband, nanoseconds10Ms}, + {BandwidthMediumband, nanoseconds20Ms}, + {BandwidthWideband, nanoseconds10Ms}, + {BandwidthWideband, nanoseconds20Ms}, + } + + for _, cfg := range configs { + _, _, lagMin, lagMax := pitchLagCodebooks(cfg.bandwidth) + lagCb, lagIcdf := pitchContourCodebooks(cfg.bandwidth, cfg.nanoseconds) + contourMax := len(lagIcdf) - 2 + + scenarios := []struct { + name string + isFirst bool + prevVoiced bool + prevLag int + lag int + }{ + {"absolute", true, false, 100, int(lagMin) + 10}, + {"relative_near", false, true, int(lagMin) + 50, int(lagMin) + 52}, + {"relative_escape", false, true, int(lagMin) + 10, int(lagMax) - 5}, + } + + for _, sc := range scenarios { + for _, contourIndex := range []uint32{0, uint32(contourMax / 2), uint32(contourMax)} { //nolint:gosec + name := fmt.Sprintf("bw%d_ns%d_%s_c%d", cfg.bandwidth, cfg.nanoseconds, sc.name, contourIndex) + t.Run(name, func(t *testing.T) { + enc := NewEncoder() + enc.previousLag = sc.prevLag + enc.isPreviousFrameVoiced = sc.prevVoiced + enc.rangeEncoder.Init() + enc.encodePitchLags(sc.lag, contourIndex, cfg.bandwidth, cfg.nanoseconds, sc.isFirst) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.previousLag = sc.prevLag + dec.isPreviousFrameVoiced = sc.prevVoiced + dec.rangeDecoder.Init(data) + _, pitchLags := dec.decodePitchLags(frameSignalTypeVoiced, cfg.bandwidth, cfg.nanoseconds, sc.isFirst) + + require.Len(t, pitchLags, subframeCount(cfg.nanoseconds)) + for k := range pitchLags { + want := clamp(int32(lagMin), int32(sc.lag+int(lagCb[contourIndex][k])), int32(lagMax)) //nolint:gosec + require.Equalf(t, int(want), pitchLags[k], "subframe %d", k) + } + require.Equal(t, sc.lag, dec.previousLag) + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange(), "range coder desync") + }) + } + } + } +} + +func TestEncodeLTPFilterRoundTrip(t *testing.T) { + codebooks := [][][]int8{ + codebookLTPFilterPeriodicityIndex0, + codebookLTPFilterPeriodicityIndex1, + codebookLTPFilterPeriodicityIndex2, + } + + for periodicity := range uint32(3) { + codebook := codebooks[periodicity] + for _, subframes := range []int{2, 4} { + name := fmt.Sprintf("p%d_sf%d", periodicity, subframes) + t.Run(name, func(t *testing.T) { + filterIndices := make([]uint32, subframes) + for i := range filterIndices { + filterIndices[i] = uint32((i*3 + 1) % len(codebook)) //nolint:gosec + } + + enc := NewEncoder() + enc.rangeEncoder.Init() + enc.encodeLTPFilter(periodicity, filterIndices) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.rangeDecoder.Init(data) + bQ7 := dec.decodeLTPFilterCoefficients(frameSignalTypeVoiced, subframes) + + require.Len(t, bQ7, subframes) + for i := range filterIndices { + require.Equalf(t, codebook[filterIndices[i]], bQ7[i], "subframe %d", i) + } + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange(), "range coder desync") + }) + } + } +} + +func TestEncodeLTPScalingRoundTrip(t *testing.T) { + want := []float32{15565, 12288, 8192} + for scaleIndex := range uint32(3) { + enc := NewEncoder() + enc.rangeEncoder.Init() + enc.encodeLTPScaling(scaleIndex) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.rangeDecoder.Init(data) + got := dec.decodeLTPScalingParameter(frameSignalTypeVoiced, true) + + require.Equal(t, want[scaleIndex], got) + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange()) + } +} + +func TestEncodeLCGSeedRoundTrip(t *testing.T) { + for seed := range uint32(4) { + enc := NewEncoder() + enc.rangeEncoder.Init() + enc.encodeLCGSeed(seed) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.rangeDecoder.Init(data) + got := dec.decodeLinearCongruentialGeneratorSeed() + + require.Equal(t, seed, got) + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange()) + } +} diff --git a/internal/silk/pulses.go b/internal/silk/pulses.go new file mode 100644 index 0000000..37849f7 --- /dev/null +++ b/internal/silk/pulses.go @@ -0,0 +1,251 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +const ( + shellCodecFrameLength = 16 + log2ShellCodecFrameLength = 4 + silkMaxPulses = 16 + nRateLevels = 10 +) + +// silkMaxPulsesTable bounds the pulse sum per partition size (2, 4, 8, 16) +// before the block is downscaled by an extra LSB (gain_quant.c thresholds). +var silkMaxPulsesTable = [4]int32{8, 10, 12, 16} //nolint:gochecknoglobals // per-partition pulse-sum bounds. + +// silkRateLevelsBITSQ5[signalType>>1] and silkPulsesPerBlockBITSQ5 give the Q5 +// bit cost of each rate level, used to pick the cheapest one. +var silkRateLevelsBITSQ5 = [2][9]int32{ //nolint:gochecknoglobals // Q5 rate-level bit costs. + {131, 74, 141, 79, 80, 138, 95, 104, 134}, + {95, 99, 91, 125, 93, 76, 123, 115, 123}, +} + +var silkPulsesPerBlockBITSQ5 = [9][18]int32{ //nolint:gochecknoglobals // Q5 pulses-per-block bit costs. + {31, 57, 107, 160, 205, 205, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255}, + {69, 47, 67, 111, 166, 205, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255}, + {82, 74, 79, 95, 109, 128, 145, 160, 173, 205, 205, 205, 224, 255, 255, 224, 255, 224}, + {125, 74, 59, 69, 97, 141, 182, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255}, + {173, 115, 85, 73, 76, 92, 115, 145, 173, 205, 224, 224, 255, 255, 255, 255, 255, 255}, + {166, 134, 113, 102, 101, 102, 107, 118, 125, 138, 145, 155, 166, 182, 192, 192, 205, 150}, + {224, 182, 134, 101, 83, 79, 85, 97, 120, 145, 173, 205, 224, 255, 255, 255, 255, 255}, + {255, 224, 192, 150, 120, 101, 92, 89, 93, 102, 118, 134, 160, 182, 192, 224, 224, 224}, + {255, 224, 224, 182, 155, 134, 118, 109, 104, 102, 106, 111, 118, 131, 145, 160, 173, 131}, +} + +// silkSignICDF holds the excitation sign thresholds, indexed by +// 7*(quantOffsetType + 2*signalType) + min(pulseCount, 6). The decoder's named +// sign tables are the same values in {256, 256-threshold, 256} form. +var silkSignICDF = [42]int32{ //nolint:gochecknoglobals // excitation sign thresholds. + 254, 49, 67, 77, 82, 93, 99, + 198, 11, 18, 24, 31, 36, 45, + 255, 46, 66, 78, 87, 94, 104, + 208, 14, 21, 32, 42, 51, 66, + 255, 94, 104, 109, 112, 115, 118, + 248, 53, 69, 80, 88, 95, 102, +} + +// encodePulses range-encodes the excitation quantization indices, inverting the +// decoder's excitation path (RFC 6716 Section 4.2.7.8). pulses holds one signed +// index per sample; frameLength is padded up to a whole number of 16-sample +// shell blocks. +// +//nolint:gocognit,gocyclo,cyclop,maintidx // faithful port of silk_encode_pulses. +func (e *Encoder) encodePulses( + signalType frameSignalType, + quantOffsetType frameQuantizationOffsetType, + pulses []int8, + frameLength int, +) { + iter := frameLength >> log2ShellCodecFrameLength + if iter*shellCodecFrameLength < frameLength { + iter++ + } + + absPulses := make([]int32, iter*shellCodecFrameLength) + for i := range absPulses { + if i < len(pulses) { + absPulses[i] = absInt32(int32(pulses[i])) + } + } + + // Sum the pulses per block, downscaling by an LSB whenever any partition + // exceeds its silkMaxPulsesTable limit. + sumPulses := make([]int32, iter) + nRshifts := make([]int32, iter) + var comb [8]int32 + for i := range iter { + block := absPulses[i*shellCodecFrameLength : (i+1)*shellCodecFrameLength] + for { + scaleDown := combineAndCheck(comb[:8], block, silkMaxPulsesTable[0], 8) + scaleDown += combineAndCheck(comb[:4], comb[:8], silkMaxPulsesTable[1], 4) + scaleDown += combineAndCheck(comb[:2], comb[:4], silkMaxPulsesTable[2], 2) + scaleDown += combineAndCheck(sumPulses[i:i+1], comb[:2], silkMaxPulsesTable[3], 1) + if scaleDown == 0 { + break + } + nRshifts[i]++ + for k := range block { + block[k] >>= 1 + } + } + } + + // Pick the rate level with the fewest bits for the pulse-count symbols. + sigTypeIndex := 0 + if signalType == frameSignalTypeVoiced { + sigTypeIndex = 1 + } + minSumBits := int32(math.MaxInt32) + rateLevel := 0 + for k := range nRateLevels - 1 { + sumBits := silkRateLevelsBITSQ5[sigTypeIndex][k] + for i := range iter { + if nRshifts[i] > 0 { + sumBits += silkPulsesPerBlockBITSQ5[k][silkMaxPulses+1] + } else { + sumBits += silkPulsesPerBlockBITSQ5[k][sumPulses[i]] + } + } + if sumBits < minSumBits { + minSumBits = sumBits + rateLevel = k + } + } + if signalType == frameSignalTypeVoiced { + e.rangeEncoder.EncodeSymbolWithICDF(icdfRateLevelVoiced, uint32(rateLevel)) //nolint:gosec // G115 + } else { + e.rangeEncoder.EncodeSymbolWithICDF(icdfRateLevelUnvoiced, uint32(rateLevel)) //nolint:gosec // G115 + } + + // Pulse counts, with the value 17 escaping to extra LSBs. + for i := range iter { + if nRshifts[i] == 0 { + e.rangeEncoder.EncodeSymbolWithICDF(icdfPulseCount[rateLevel], uint32(sumPulses[i])) //nolint:gosec // G115 + } else { + e.rangeEncoder.EncodeSymbolWithICDF(icdfPulseCount[rateLevel], silkMaxPulses+1) + for range nRshifts[i] - 1 { + e.rangeEncoder.EncodeSymbolWithICDF(icdfPulseCount[nRateLevels-1], silkMaxPulses+1) + } + e.rangeEncoder.EncodeSymbolWithICDF(icdfPulseCount[nRateLevels-1], uint32(sumPulses[i])) //nolint:gosec // G115 + } + } + + // Pulse locations (shell code) for blocks that have pulses. + for i := range iter { + if sumPulses[i] > 0 { + e.encodePulseLocation(absPulses[i*shellCodecFrameLength:(i+1)*shellCodecFrameLength], sumPulses[i]) + } + } + + // LSBs shifted out during downscaling, most significant first. + for i := range iter { + if nRshifts[i] == 0 { + continue + } + nLS := nRshifts[i] - 1 + for k := range shellCodecFrameLength { + // LSBs come from the original magnitude, not the downscaled block. + absQ := int32(0) + if idx := i*shellCodecFrameLength + k; idx < len(pulses) { + absQ = absInt32(int32(pulses[idx])) + } + for j := nLS; j > 0; j-- { + e.rangeEncoder.EncodeSymbolWithICDF(icdfExcitationLSB, uint32((absQ>>j)&1)) //nolint:gosec // G115 + } + e.rangeEncoder.EncodeSymbolWithICDF(icdfExcitationLSB, uint32(absQ&1)) //nolint:gosec // G115 + } + } + + e.encodeExcitationSigns(pulses, frameLength, signalType, quantOffsetType, sumPulses) +} + +// combineAndCheck sums adjacent pairs of in into out and reports whether any +// sum exceeds maxPulses (in which case the block must be downscaled). +func combineAndCheck(out, in []int32, maxPulses int32, length int) int32 { + scaleDown := int32(0) + for k := range length { + sum := in[2*k] + in[2*k+1] + if sum > maxPulses { + scaleDown = 1 + } + out[k] = sum + } + + return scaleDown +} + +// encodePulseLocation is the shell encoder: it recursively splits a 16-sample +// block and codes the left-child pulse count at each split, inverting +// Decoder.decodePulseLocation. +func (e *Encoder) encodePulseLocation(block []int32, total int32) { + var s2 [8]int32 + for t := range 8 { + s2[t] = block[2*t] + block[2*t+1] + } + var s4 [4]int32 + for u := range 4 { + s4[u] = s2[2*u] + s2[2*u+1] + } + var s8 [2]int32 + for v := range 2 { + s8[v] = s4[2*v] + s4[2*v+1] + } + + e.encodeSplit(icdfPulseCountSplit16SamplePartitions, total, s8[0]) + for j := range 2 { + e.encodeSplit(icdfPulseCountSplit8SamplePartitions, s8[j], s4[2*j]) + for k := range 2 { + idx4 := 2*j + k + e.encodeSplit(icdfPulseCountSplit4SamplePartitions, s4[idx4], s2[2*idx4]) + for l := range 2 { + idx2 := 2*idx4 + l + e.encodeSplit(icdfPulseCountSplit2SamplePartitions, s2[idx2], block[2*idx2]) + } + } + } +} + +// encodeSplit codes how many of a partition's pulses fall in its first half. +func (e *Encoder) encodeSplit(icdf [][]uint, block, leftChild int32) { + if block > 0 { + e.rangeEncoder.EncodeSymbolWithICDF(icdf[block-1], uint32(leftChild)) //nolint:gosec // G115 + } +} + +// encodeExcitationSigns codes a sign for every non-zero pulse, using the PDF +// selected by signal type, quantization offset type and block pulse count. +func (e *Encoder) encodeExcitationSigns( + pulses []int8, + frameLength int, + signalType frameSignalType, + quantOffsetType frameQuantizationOffsetType, + sumPulses []int32, +) { + blocks := (frameLength + shellCodecFrameLength/2) >> log2ShellCodecFrameLength + offsetIndex := int(quantOffsetType) - int(frameQuantizationOffsetTypeLow) + signalIndex := int(signalType) - int(frameSignalTypeInactive) + base := 7 * (offsetIndex + 2*signalIndex) + for i := range blocks { + p := sumPulses[i] + if p <= 0 { + continue + } + count := int(p) + count = min(count, 6) + icdf := []uint{256, uint(256 - silkSignICDF[base+count]), 256} //nolint:gosec // G115 + for j := range shellCodecFrameLength { + idx := i*shellCodecFrameLength + j + if idx >= len(pulses) || pulses[idx] == 0 { + continue + } + sym := uint32(0) + if pulses[idx] > 0 { + sym = 1 + } + e.rangeEncoder.EncodeSymbolWithICDF(icdf, sym) + } + } +} diff --git a/internal/silk/pulses_test.go b/internal/silk/pulses_test.go new file mode 100644 index 0000000..41c3c59 --- /dev/null +++ b/internal/silk/pulses_test.go @@ -0,0 +1,70 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "fmt" + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// lcgPulses generates deterministic signed pulses in [-maxMag, maxMag]. +func lcgPulses(seed uint32, n int, maxMag int32) []int8 { + pulses := make([]int8, n) + state := seed + for i := range pulses { + state = 1664525*state + 1013904223 + pulses[i] = int8(int32(state>>24)%(2*maxMag+1) - maxMag) //nolint:gosec // G115: bounded to [-maxMag,maxMag]. + } + + return pulses +} + +// TestEncodePulsesRoundTrip encodes quantized excitation indices and decodes +// them back through the real decoder path, asserting the signed pulses and the +// range coder state match. Higher magnitudes force LSB downscaling, exercising +// the escape and LSB paths. +func TestEncodePulsesRoundTrip(t *testing.T) { + signalTypes := []frameSignalType{frameSignalTypeInactive, frameSignalTypeUnvoiced, frameSignalTypeVoiced} + offsets := []frameQuantizationOffsetType{frameQuantizationOffsetTypeLow, frameQuantizationOffsetTypeHigh} + lengths := []int{80, 160, 320} + mags := []int32{1, 3, 8} + + for _, signalType := range signalTypes { + for _, offset := range offsets { + for _, length := range lengths { + for _, mag := range mags { + name := fmt.Sprintf("t%d_o%d_len%d_mag%d", signalType, offset, length, mag) + t.Run(name, func(t *testing.T) { + seed := uint32(int(signalType)*131 + int(offset)*17 + length + int(mag)) //nolint:gosec // G115 + pulses := lcgPulses(seed, length, mag) + + enc := NewEncoder() + enc.rangeEncoder.Init() + enc.encodePulses(signalType, offset, pulses, length) + encRange := enc.rangeEncoder.FinalRange() + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + dec.rangeDecoder.Init(data) + shellblocks := length / shellCodecFrameLength + rateLevel := dec.decodeRatelevel(signalType == frameSignalTypeVoiced) + pulsecounts, lsbcounts := dec.decodePulseAndLSBCounts(shellblocks, rateLevel) + eRaw := dec.decodePulseLocation(pulsecounts) + dec.decodeExcitationLSB(eRaw, lsbcounts) + dec.decodeExcitationSign(eRaw, signalType, offset, pulsecounts) + + require.Len(t, eRaw, len(pulses)) + for i := range pulses { + require.Equalf(t, int32(pulses[i]), eRaw[i], "pulse %d", i) + } + assert.Equal(t, encRange, dec.rangeDecoder.FinalRange(), "range coder desync") + }) + } + } + } + } +}