diff --git a/encoder.go b/encoder.go index 2517efb..342d27c 100644 --- a/encoder.go +++ b/encoder.go @@ -8,6 +8,7 @@ import ( "fmt" "github.com/pion/opus/internal/celt" + "github.com/pion/opus/internal/silk" ) const ( @@ -22,9 +23,18 @@ const ( // is separate (bit 2 of the TOC) and not part of this constant. const celtOnlyFullband20msConfig = 31 +// silkOnlyWideband20msConfig is the TOC config number for SILK-only, wideband, +// 20 ms frames (RFC 6716 Table 2). +const silkOnlyWideband20msConfig = 9 + +// silkWidebandSampleCount is the number of 16 kHz samples in a 20 ms SILK +// wideband frame. +const silkWidebandSampleCount = 320 + // Encoder encodes PCM into Opus packets. type Encoder struct { celtEncoder celt.Encoder + silkEncoder silk.Encoder sampleRate int channels int bitrate int @@ -100,6 +110,7 @@ func WithComplexity(complexity int) EncoderOption { func NewEncoder(opts ...EncoderOption) (*Encoder, error) { encoder := &Encoder{ celtEncoder: celt.NewEncoder(), + silkEncoder: silk.NewEncoder(), sampleRate: celtSampleRate, channels: 1, bitrate: defaultBitrate, @@ -209,6 +220,41 @@ func splitChannels(in []float32, numChannels, frameSamples int) [][]float32 { return ch } +// EncodeSILK encodes one 20 ms mono SILK frame into a SILK-only Opus packet. +// pcm must hold exactly one 20 ms frame of mono s16 samples at the bandwidth's +// internal rate: 160 (NB/8 kHz), 240 (MB/12 kHz), or 320 (WB/16 kHz). The +// encoder covers voiced/LTP prediction, noise shaping and NLSF interpolation +// (see internal/silk); the delayed-decision NSQ, stereo and hybrid mode are +// not yet implemented. +func (e *Encoder) EncodeSILK(pcm []int16, bandwidth Bandwidth, out []byte) (int, error) { + var config int + switch bandwidth { + case BandwidthNarrowband: + config = 1 + case BandwidthMediumband: + config = 5 + case BandwidthWideband: + config = silkOnlyWideband20msConfig + default: + return 0, fmt.Errorf("%w: bandwidth %d", errInvalidSampleRate, bandwidth) + } + + want := bandwidth.SampleRate() / 50 // 20 ms + if len(pcm) != want { + return 0, fmt.Errorf("%w: got %d samples, want %d", errInvalidFrameSize, len(pcm), want) + } + + payload := e.silkEncoder.Encode(pcm, silk.Bandwidth(bandwidth), e.bitrate) + if len(out) < len(payload)+1 { + return 0, errOutBufferTooSmall + } + + out[0] = byte(config<<3) | byte(frameCodeOneFrame) // mono, one frame + n := copy(out[1:], payload) + + return n + 1, nil +} + func (e *Encoder) frameBytes() int { return int(int64(e.bitrate) * frame20msNS / 1000000000 / 8) } diff --git a/internal/silk/a2nlsf.go b/internal/silk/a2nlsf.go new file mode 100644 index 0000000..9512c59 --- /dev/null +++ b/internal/silk/a2nlsf.go @@ -0,0 +1,196 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// LPC-to-NLSF conversion (silk_A2NLSF, A2NLSF.c). The whitening filter's even +// and odd polynomials are root-found on a cosine grid to recover the NLSFs. + +const ( + binDivStepsA2NLSF = 3 + maxIterA2NLSF = 16 + lsfCosTabSz = 128 +) + +// silkLSFCosTabFIXQ12 is the Q12 cosine table used for the root search. +// +//nolint:gochecknoglobals +var silkLSFCosTabFIXQ12 = [lsfCosTabSz + 1]int16{ + 8192, 8190, 8182, 8170, 8152, 8130, 8104, 8072, 8034, 7994, 7946, 7896, + 7840, 7778, 7714, 7644, 7568, 7490, 7406, 7318, 7226, 7128, 7026, 6922, + 6812, 6698, 6580, 6458, 6332, 6204, 6070, 5934, 5792, 5648, 5502, 5352, + 5198, 5040, 4880, 4718, 4552, 4382, 4212, 4038, 3862, 3684, 3502, 3320, + 3136, 2948, 2760, 2570, 2378, 2186, 1990, 1794, 1598, 1400, 1202, 1002, + 802, 602, 402, 202, 0, -202, -402, -602, -802, -1002, -1202, -1400, + -1598, -1794, -1990, -2186, -2378, -2570, -2760, -2948, -3136, -3320, + -3502, -3684, -3862, -4038, -4212, -4382, -4552, -4718, -4880, -5040, + -5198, -5352, -5502, -5648, -5792, -5934, -6070, -6204, -6332, -6458, + -6580, -6698, -6812, -6922, -7026, -7128, -7226, -7318, -7406, -7490, + -7568, -7644, -7714, -7778, -7840, -7896, -7946, -7994, -8034, -8072, + -8104, -8130, -8152, -8170, -8182, -8190, -8192, +} + +// smlaww returns a + ((b * c) >> 16). +func smlaww(a, b, c int32) int32 { + return a + smulww(b, c) +} + +// a2nlsfTransPoly transforms a polynomial from cos(n*f) to cos(f)^n. +func a2nlsfTransPoly(p []int32, dd int) { + for k := 2; k <= dd; k++ { + for n := dd; n > k; n-- { + p[n-2] -= p[n] + } + p[k-2] -= p[k] << 1 + } +} + +// a2nlsfEvalPoly evaluates a Q16 polynomial at the Q12 point x. +func a2nlsfEvalPoly(p []int32, x int32, dd int) int32 { + y32 := p[dd] + xQ16 := x << 4 + for n := dd - 1; n >= 0; n-- { + y32 = smlaww(p[n], y32, xQ16) + } + + return y32 +} + +// a2nlsfInit builds the even (P) and odd (Q) polynomials from the LPC coefs. +func a2nlsfInit(aQ16, p, q []int32, dd int) { + p[dd] = 1 << 16 + q[dd] = 1 << 16 + for k := range dd { + p[k] = -aQ16[dd-k-1] - aQ16[dd+k] + q[k] = -aQ16[dd-k-1] + aQ16[dd+k] + } + for k := dd; k > 0; k-- { + p[k-1] -= p[k] + q[k-1] += q[k] + } + a2nlsfTransPoly(p, dd) + a2nlsfTransPoly(q, dd) +} + +// bwexpander32 applies a Q16 chirp to int32 AR coefficients. +func bwexpander32(ar []int32, d int, chirpQ16 int32) { + chirpMinusOneQ16 := chirpQ16 - 65536 + for i := range d - 1 { + ar[i] = smulww(chirpQ16, ar[i]) + chirpQ16 += rshiftRound32(chirpQ16*chirpMinusOneQ16, 16) + } + ar[d-1] = smulww(chirpQ16, ar[d-1]) +} + +// a2nlsf converts monic whitening filter coefficients (Q16, modified in place) +// to Q15 NLSFs, bandwidth-expanding until the roots converge. +// +//nolint:gocognit,gocyclo,cyclop,maintidx,varnamelen // faithful port of silk_A2NLSF. +func a2nlsf(nlsf []int16, aQ16 []int32, d int) { + dd := d >> 1 + pPoly := make([]int32, dd+1) + qPoly := make([]int32, dd+1) + a2nlsfInit(aQ16, pPoly, qPoly, dd) + polys := [2][]int32{pPoly, qPoly} + + p := pPoly //nolint:varnamelen // p tracks the active P/Q polynomial, as in the C reference. + xlo := int32(silkLSFCosTabFIXQ12[0]) + ylo := a2nlsfEvalPoly(p, xlo, dd) + + var rootIx int + if ylo < 0 { + nlsf[0] = 0 + p = qPoly + ylo = a2nlsfEvalPoly(p, xlo, dd) + rootIx = 1 + } + + k := 1 //nolint:varnamelen // k indexes the cosine grid, as in the C reference. + iterations := 0 + thr := int32(0) + for { + xhi := int32(silkLSFCosTabFIXQ12[k]) + yhi := a2nlsfEvalPoly(p, xhi, dd) + + if (ylo <= 0 && yhi >= thr) || (ylo >= 0 && yhi <= -thr) { //nolint:nestif // faithful port of silk_A2NLSF. + if yhi == 0 { + thr = 1 + } else { + thr = 0 + } + ffrac := int32(-256) + for m := range binDivStepsA2NLSF { + xmid := rshiftRound32(xlo+xhi, 1) + ymid := a2nlsfEvalPoly(p, xmid, dd) + if (ylo <= 0 && ymid >= 0) || (ylo >= 0 && ymid <= 0) { + xhi = xmid + yhi = ymid + } else { + xlo = xmid + ylo = ymid + ffrac += 128 >> m + } + } + if absInt32(ylo) < 65536 { + den := ylo - yhi + nom := (ylo << (8 - binDivStepsA2NLSF)) + (den >> 1) + if den != 0 { + ffrac += nom / den + } + } else { + ffrac += ylo / ((ylo - yhi) >> (8 - binDivStepsA2NLSF)) + } + nlsf[rootIx] = int16(min((int32(k)<<8)+ffrac, math.MaxInt16)) + + rootIx++ + if rootIx >= d { + return + } + p = polys[rootIx&1] + xlo = int32(silkLSFCosTabFIXQ12[k-1]) + ylo = int32(1-(rootIx&2)) << 12 + } else { + k++ + xlo = xhi + ylo = yhi + thr = 0 + + if k > lsfCosTabSz { + iterations++ + if iterations > maxIterA2NLSF { + nlsf[0] = int16((1 << 15) / int32(d+1)) //nolint:gosec // G115: quotient fits int16. + for k = 1; k < d; k++ { + nlsf[k] = nlsf[k-1] + nlsf[0] + } + + return + } + bwexpander32(aQ16, d, 65536-(1< +// SPDX-License-Identifier: MIT + +package silk + +import ( + "fmt" + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// TestA2NLSFRoundTrip converts a valid NLSF vector to LPC via the decoder and +// back with A2NLSF, expecting to recover the input within the fixed-point +// tolerance. +func TestA2NLSFRoundTrip(t *testing.T) { + cases := []struct { + bandwidth Bandwidth + order int + }{ + {BandwidthNarrowband, 10}, + {BandwidthWideband, 16}, + } + + for _, tc := range cases { + for seed := range 8 { + t.Run(fmt.Sprintf("bw%d_seed%d", tc.bandwidth, seed), func(t *testing.T) { + nlsf := genNLSF(tc.order, uint32(seed*131+tc.order)) //nolint:gosec // G115: non-negative test seed. + stabilizeNLSF(nlsf, tc.order, tc.bandwidth) + + dec := NewDecoder() + a32Q17 := dec.convertNormalizedLSFsToLPCCoefficients(nlsf, tc.bandwidth) + aFloat := make([]float32, tc.order) + for i := range aFloat { + aFloat[i] = float32(a32Q17[i]) / (1 << 17) + } + + recovered := make([]int16, tc.order) + a2nlsfFLP(recovered, aFloat, tc.order) + + // Round-trip through two different fixed-point approximations + // (decoder NLSF->LPC and A2NLSF LPC->NLSF); band-edge roots are + // the least precise. ~1% of full scale is expected. + for k := range nlsf { + assert.InDeltaf(t, nlsf[k], recovered[k], 300, "coefficient %d", k) + } + // Result must be a valid increasing NLSF vector. + for k := 1; k < tc.order; k++ { + require.Greaterf(t, recovered[k], recovered[k-1], "not increasing at %d", k) + } + }) + } + } +} diff --git a/internal/silk/control_snr.go b/internal/silk/control_snr.go new file mode 100644 index 0000000..5487687 --- /dev/null +++ b/internal/silk/control_snr.go @@ -0,0 +1,77 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// Target residual-quantizer SNR as a function of bitrate (silk_control_SNR, +// control_SNR.c). + +//nolint:gochecknoglobals // bitrate->SNR lookup tables (control_SNR.c). +var ( + silkTargetRateNB21 = []uint8{ + 0, 15, 39, 52, 61, 68, + 74, 79, 84, 88, 92, 95, 99, 102, 105, 108, 111, 114, 117, 119, 122, 124, + 126, 129, 131, 133, 135, 137, 139, 142, 143, 145, 147, 149, 151, 153, 155, 157, + 158, 160, 162, 163, 165, 167, 168, 170, 171, 173, 174, 176, 177, 179, 180, 182, + 183, 185, 186, 187, 189, 190, 192, 193, 194, 196, 197, 199, 200, 201, 203, 204, + 205, 207, 208, 209, 211, 212, 213, 215, 216, 217, 219, 220, 221, 223, 224, 225, + 227, 228, 230, 231, 232, 234, 235, 236, 238, 239, 241, 242, 243, 245, 246, 248, + 249, 250, 252, 253, 255, + } + + silkTargetRateMB21 = []uint8{ + 0, 0, 28, 43, 52, 59, + 65, 70, 74, 78, 81, 85, 87, 90, 93, 95, 98, 100, 102, 105, 107, 109, + 111, 113, 115, 116, 118, 120, 122, 123, 125, 127, 128, 130, 131, 133, 134, 136, + 137, 138, 140, 141, 143, 144, 145, 147, 148, 149, 151, 152, 153, 154, 156, 157, + 158, 159, 160, 162, 163, 164, 165, 166, 167, 168, 169, 171, 172, 173, 174, 175, + 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 188, 189, 190, + 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 203, 204, 205, + 206, 207, 208, 209, 210, 211, 212, 213, 214, 214, 215, 216, 217, 218, 219, 220, + 221, 222, 223, 224, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, + 236, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, + 251, 252, 253, 254, 255, + } + + silkTargetRateWB21 = []uint8{ + 0, 0, 0, 8, 29, 41, + 49, 56, 62, 66, 70, 74, 77, 80, 83, 86, 88, 91, 93, 95, 97, 99, + 101, 103, 105, 107, 108, 110, 112, 113, 115, 116, 118, 119, 121, 122, 123, 125, + 126, 127, 129, 130, 131, 132, 134, 135, 136, 137, 138, 140, 141, 142, 143, 144, + 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 156, 157, 158, 159, 159, 160, + 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 171, 172, 173, 174, 175, + 176, 177, 177, 178, 179, 180, 181, 181, 182, 183, 184, 185, 185, 186, 187, 188, + 189, 189, 190, 191, 192, 192, 193, 194, 195, 195, 196, 197, 198, 198, 199, 200, + 200, 201, 202, 203, 203, 204, 205, 206, 206, 207, 208, 209, 209, 210, 211, 211, + 212, 213, 214, 214, 215, 216, 216, 217, 218, 219, 219, 220, 221, 221, 222, 223, + 224, 224, 225, 226, 226, 227, 228, 229, 229, 230, 231, 232, 232, 233, 234, 234, + 235, 236, 237, 237, 238, 239, 240, 240, 241, 242, 243, 243, 244, 245, 246, 246, + 247, 248, 249, 249, 250, 251, 252, 253, 255, + } +) + +// controlSNR returns the target quantizer SNR in Q7 dB for a bitrate and +// bandwidth (silk_control_SNR). +func controlSNR(fsKHz, nbSubfr, targetRateBps int) int32 { + if nbSubfr == 2 { + targetRateBps -= 2000 + fsKHz/16 + } + + var table []uint8 + switch fsKHz { + case 8: + table = silkTargetRateNB21 + case 12: + table = silkTargetRateMB21 + default: + table = silkTargetRateWB21 + } + + id := (targetRateBps + 200) / 400 + id = min(id-10, len(table)-1) + if id <= 0 { + return 0 + } + + return int32(table[id]) * 21 +} diff --git a/internal/silk/decoder.go b/internal/silk/decoder.go index 8cd0b40..71de104 100644 --- a/internal/silk/decoder.go +++ b/internal/silk/decoder.go @@ -616,147 +616,8 @@ func (d *Decoder) normalizeLineSpectralFrequencyCoefficients( // percentile of a large training set). // // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 -// -//nolint:cyclop func (d *Decoder) normalizeLSFStabilization(nlsfQ15 []int16, dLPC int, bandwidth Bandwidth) { - // Let NDeltaMin_Q15[k] be the minimum required spacing for the current - // audio bandwidth from Table 25. - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - NDeltaMinQ15 := codebookMinimumSpacingForNormalizedLSCoefficientsNarrowbandAndMediumband - if bandwidth == BandwidthWideband { - NDeltaMinQ15 = codebookMinimumSpacingForNormalizedLSCoefficientsWideband - } - - // The procedure starts off by trying to make small adjustments that - // attempt to minimize the amount of distortion introduced. After 20 - // such adjustments, it falls back to a more direct method that - // guarantees the constraints are enforced but may require large - // adjustments. - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - for adjustment := 0; adjustment <= 19; adjustment++ { - // First, the procedure finds the index - // i where NLSF_Q15[i] - NLSF_Q15[i-1] - NDeltaMin_Q15[i] is the - // smallest, breaking ties by using the lower value of i. - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - i := 0 - iValue := int(math.MaxInt) - - for nlsfIndex := 0; nlsfIndex <= len(nlsfQ15); nlsfIndex++ { - // For the purposes of computing this spacing for the first and last coefficient, - // NLSF_Q15[-1] is taken to be 0 and NLSF_Q15[d_LPC] is taken to be 32768 - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - previousNLSF := 0 - currentNLSF := 32768 - if nlsfIndex != 0 { - previousNLSF = int(nlsfQ15[nlsfIndex-1]) // #nosec G602 - } - if nlsfIndex != len(nlsfQ15) { - currentNLSF = int(nlsfQ15[nlsfIndex]) - } - - spacingValue := currentNLSF - previousNLSF - NDeltaMinQ15[nlsfIndex] - if spacingValue < iValue { - i = nlsfIndex - iValue = spacingValue - } - } - - switch { - // If this value is non-negative, then the stabilization stops; the coefficients - // satisfy all the constraints. - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - case iValue >= 0: - return - // if i == 0, it sets NLSF_Q15[0] to NDeltaMin_Q15[0] - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - case i == 0: - nlsfQ15[0] = int16(NDeltaMinQ15[0]) //nolint:gosec // G115 - - continue - // if i == d_LPC, it sets - // NLSF_Q15[d_LPC-1] to (32768 - NDeltaMin_Q15[d_LPC]) - // - // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.4 - case i == dLPC: - nlsfQ15[dLPC-1] = int16(32768 - NDeltaMinQ15[dLPC]) //nolint:gosec // G115 - - continue - } - - // For all other values of i, both NLSF_Q15[i-1] and NLSF_Q15[i] are updated as - // follows: - // i-1 - // __ - // min_center_Q15 = (NDeltaMin_Q15[i]>>1) + \ NDeltaMin_Q15[k] - // /_ - // k=0 - // - minCenterQ15 := NDeltaMinQ15[i] >> 1 - for k := 0; k <= i-1; k++ { - minCenterQ15 += NDeltaMinQ15[k] - } - - // d_LPC - // __ - // max_center_Q15 = 32768 - (NDeltaMin_Q15[i]>>1) - \ NDeltaMin_Q15[k] - // /_ - // k=i+1 - maxCenterQ15 := 32768 - (NDeltaMinQ15[i] >> 1) - for k := i + 1; k <= dLPC; k++ { - maxCenterQ15 -= NDeltaMinQ15[k] - } - - // center_freq_Q15 = clamp(min_center_Q15[i], - // (NLSF_Q15[i-1] + NLSF_Q15[i] + 1)>>1 - // max_center_Q15[i]) - centerFreqQ15 := int(clamp( - int32(minCenterQ15), //nolint:gosec // G115 - int32((int(nlsfQ15[i-1])+int(nlsfQ15[i])+1)>>1), //nolint:gosec // G115 - int32(maxCenterQ15)), //nolint:gosec // G115 - ) - - // NLSF_Q15[i-1] = center_freq_Q15 - (NDeltaMin_Q15[i]>>1) - // NLSF_Q15[i] = NLSF_Q15[i-1] + NDeltaMin_Q15[i] - nlsfQ15[i-1] = int16(centerFreqQ15 - NDeltaMinQ15[i]>>1) //nolint:gosec // G115 - nlsfQ15[i] = nlsfQ15[i-1] + int16(NDeltaMinQ15[i]) //nolint:gosec // G115 - } - - // After the 20th repetition of the above procedure, the following - // fallback procedure executes once. First, the values of NLSF_Q15[k] - // for 0 <= k < d_LPC are sorted in ascending order. Then, for each - // value of k from 0 to d_LPC-1, NLSF_Q15[k] is set to - slices.Sort(nlsfQ15) - - // Then, for each value of k from 0 to d_LPC-1, NLSF_Q15[k] is set to - // - // max(NLSF_Q15[k], NLSF_Q15[k-1] + NDeltaMin_Q15[k]) - for k := 0; k <= dLPC-1; k++ { - prevNLSF := int16(0) - if k != 0 { - prevNLSF = nlsfQ15[k-1] // #nosec G602 - } - - nlsfQ15[k] = maxInt16(nlsfQ15[k], saturatingAddInt16(prevNLSF, int16(NDeltaMinQ15[k]))) //nolint:gosec // G115 - } - - // Next, for each value of k from d_LPC-1 down to 0, NLSF_Q15[k] is set - // to - // - // min(NLSF_Q15[k], NLSF_Q15[k+1] - NDeltaMin_Q15[k+1]) - for k := dLPC - 1; k >= 0; k-- { - nextNLSF := 32768 - if k != dLPC-1 { - nextNLSF = int(nlsfQ15[k+1]) - } - - nlsfQ15[k] = minInt16(nlsfQ15[k], int16(nextNLSF-NDeltaMinQ15[k+1])) //nolint:gosec // G115 - } + stabilizeNLSF(nlsfQ15, dLPC, bandwidth) } // https://datatracker.ietf.org/doc/html/rfc6716#section-4.2.7.5.5 diff --git a/internal/silk/enc_frame.go b/internal/silk/enc_frame.go new file mode 100644 index 0000000..32ea012 --- /dev/null +++ b/internal/silk/enc_frame.go @@ -0,0 +1,284 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// This file assembles a full SILK frame: analysis, quantization, the NSQ, and +// range coding of every field in decode order. It encodes mono, 20 ms, +// SILK-only frames with voiced/LTP prediction, faithful noise shaping and NLSF +// interpolation. The delayed-decision NSQ, stereo, and the rate-control loop +// are follow-up refinements. + +const ( + silkGainFloor = 1.0 + silkGainCeil = 32767.0 + silkVADThreshold = 100 // speech_activity_Q8 above which a frame is treated as active + silkLambdaQ10 = 1024 + silkLTPScaleQ14 = 15565 + burgMinInvGain = 1e-4 // 1 / MAX_PREDICTION_POWER_GAIN +) + +// silkInternalRate returns the SILK internal sample rate in kHz. +func silkInternalRate(bandwidth Bandwidth) int { + switch bandwidth { + case BandwidthNarrowband: + return 8 + case BandwidthMediumband: + return 12 + default: + return 16 + } +} + +// silkLPCOrder returns the prediction LPC order for the bandwidth. +func silkLPCOrder(bandwidth Bandwidth) int { + if bandwidth == BandwidthWideband { + return 16 + } + + return 10 +} + +// nlsfToLPCQ12 reconstructs Q12 LPC coefficients from quantized NLSFs, reusing +// the decoder's NLSF->LPC path so the encoder predicts with the exact +// coefficients the decoder will use. +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 +} + +// Encode encodes one 20 ms mono SILK frame from internal-rate PCM and returns +// the range-coded SILK payload (the SILK header plus frame, without the Opus +// TOC byte). +func (e *Encoder) Encode(input []int16, bandwidth Bandwidth, targetBitrate int) []byte { + if targetBitrate > 0 { + e.targetBitrate = targetBitrate + } + e.rangeEncoder.Init() + e.encodeSILKFrame(input, bandwidth) + + return e.rangeEncoder.Done() +} + +// encodeSILKFrame encodes one 20 ms mono SILK frame to the range encoder. +// input is PCM at the internal rate for the bandwidth. +// +//nolint:gocyclo,cyclop // the frame encoder threads many stages in decode order. +func (e *Encoder) encodeSILKFrame(input []int16, bandwidth Bandwidth) { + fsKHz := silkInternalRate(bandwidth) + order := silkLPCOrder(bandwidth) + subfrCount := subframeCount(nanoseconds20Ms) + subfrLength := 5 * fsKHz + frameLength := subfrCount * subfrLength + ltpMemLength := 20 * fsKHz + + // Voice activity. + saQ8, tiltQ15, quality := e.vad.getSpeechActivityQ8(input, frameLength, fsKHz) + active := saQ8 > silkVADThreshold + + // Pitch analysis on the whitening residual (with LTP-memory history). + if len(e.xBuf) != ltpMemLength { + e.xBuf = make([]float32, ltpMemLength) + } + analysis := make([]float32, ltpMemLength+frameLength+ltpOrder) + copy(analysis, e.xBuf) + for i := range frameLength { + analysis[ltpMemLength+i] = float32(input[i]) + } + voiced, pitchL, lagIndex, contourIndex, res, predGain := e.findPitchLags( + analysis[:ltpMemLength+frameLength], fsKHz, subfrCount, saQ8, tiltQ15) + // Keep a few zero samples of headroom after the residual for find_LTP. + res = append(res, make([]float32, ltpOrder)...) + + signalType := frameSignalTypeInactive + switch { + case voiced: + signalType = frameSignalTypeVoiced + active = true + case active: + signalType = frameSignalTypeUnvoiced + } + + // Noise-shaping analysis: AR shaping filters, initial gains, spectral tilt, + // low-frequency and harmonic shaping. + snrDBQ7 := controlSNR(fsKHz, subfrCount, e.targetBitrate) + laShape := laShapeMSLowComplex * fsKHz + shapeBuf := make([]float32, laShape+frameLength+laShape) + copy(shapeBuf, e.xBuf[ltpMemLength-laShape:ltpMemLength]) + for i := range frameLength { + shapeBuf[laShape+i] = float32(input[i]) + } + sr := e.noiseShapeAnalysis( + shapeBuf, signalType, pitchL, predGain, snrDBQ7, saQ8, quality, fsKHz, subfrCount, subfrLength) + + // Prediction coefficients (find_pred_coefs). Build LPC_in_pre: the LTP + // residual for voiced, or the gain-normalized input for unvoiced. Both drive + // the short-term LPC and the residual energy. + invGains := make([]float32, subfrCount) + for k := range subfrCount { + invGains[k] = 1.0 / sr.gains[k] + } + ltpCoefQ14 := make([]int16, ltpOrder*subfrCount) + nsqPitchL := make([]int, subfrCount) + lpcInPre := make([]float32, subfrCount*(order+subfrLength)) + xBase := ltpMemLength - order + var periodicityIndex int + var filterIndices []int8 + var predGainDB float32 + ltpScaleIndex := 0 + ltpScaleQ14 := int32(silkLTPScaleQ14) + if voiced { + xxLTP := make([]float32, subfrCount*ltpMatrixSize) + xXLTP := make([]float32, subfrCount*ltpOrder) + findLTPFLP(xxLTP, xXLTP, res, ltpMemLength, pitchL, subfrLength, subfrCount) + ltpCoefQ14, filterIndices, periodicityIndex, predGainDB = e.quantLTPGains(xxLTP, xXLTP, subfrLength, subfrCount) + copy(nsqPitchL, pitchL) + ltpScaleIndex, ltpScaleQ14 = ltpScaleControl(predGainDB, snrDBQ7, e.packetLossPerc, 1, false) + + ltpCoefFloat := make([]float32, ltpOrder*subfrCount) + for i := range ltpCoefFloat { + ltpCoefFloat[i] = float32(ltpCoefQ14[i]) * (1.0 / 16384.0) + } + ltpAnalysisFilterFLP(lpcInPre, analysis, xBase, ltpCoefFloat, pitchL, invGains, subfrLength, subfrCount, order) + } else { + e.sumLogGainQ7 = 0 + for k := range subfrCount { + dst := k * (order + subfrLength) + src := xBase + k*subfrLength + for i := range order + subfrLength { + lpcInPre[dst+i] = analysis[src+i] * invGains[k] + } + } + } + + // Short-term prediction: Burg over LPC_in_pre, search the NLSF interpolation + // factor, then quantize and build both frame-half LPC sets. + minInvGain := predCoefsMinInvGain(e.firstFrameAfterReset, predGainDB, sr.codingQuality) + nlsfInterpQ2, nlsf := e.findLPCNLSF(lpcInPre, minInvGain, bandwidth, order, subfrCount, subfrLength) + stabilizeNLSF(nlsf, order, bandwidth) + index1, indices2, quantNLSF := quantizeNLSF(nlsf, bandwidth) + predCoefQ12 := nlsfToLPCQ12(quantNLSF, bandwidth) // second frame half + predCoefQ12Half0 := predCoefQ12 + if nlsfInterpQ2 < 4 { + nlsf0 := make([]int16, order) + interpolateNLSF(nlsf0, e.prevNLSFq, quantNLSF, nlsfInterpQ2, order) + predCoefQ12Half0 = nlsfToLPCQ12(nlsf0, bandwidth) // interpolated first half + } + predCoef2 := make([]int16, 2*maxLPCOrder) + copy(predCoef2, predCoefQ12Half0) + copy(predCoef2[maxLPCOrder:], predCoefQ12) + copy(e.prevNLSFq, quantNLSF) + + // Residual energy per subframe from the quantized LPC (gain soft-limit). + predCoefFloat0 := make([]float32, order) + predCoefFloat1 := make([]float32, order) + for j := range order { + predCoefFloat0[j] = float32(predCoefQ12Half0[j]) * (1.0 / 4096.0) + predCoefFloat1[j] = float32(predCoefQ12[j]) * (1.0 / 4096.0) + } + resNrg := make([]float32, subfrCount) + residualEnergyFLP(resNrg, lpcInPre, predCoefFloat0, predCoefFloat1, sr.gains, subfrLength, subfrCount, order) + + // Process gains: reduce for high LTP gain, soft-limit, quantize; Lambda + offset. + gainsQ16Int, gainIndices, lambdaQ10, quantOffsetType := e.processGains( + sr, resNrg, signalType, predGainDB, snrDBQ7, saQ8, tiltQ15, subfrLength, subfrCount, false) + + // Noise-shaping quantization. + pulses := make([]int8, frameLength) + seed := uint32(e.frameCounter & 3) //nolint:gosec // G115 + e.nsq.quantize(input, pulses, &nsqParams{ + predCoefQ12: predCoef2, + ltpCoefQ14: ltpCoefQ14, + arQ13: sr.arQ13, + harmShapeGainQ14: sr.harmShapeQ14, + tiltQ14: sr.tiltQ14, + lfShpQ14: sr.lfShpQ14, + gainsQ16: gainsQ16Int, + pitchL: nsqPitchL, + lambdaQ10: lambdaQ10, + ltpScaleQ14: ltpScaleQ14, + seed: int32(seed), //nolint:gosec // G115 + signalType: signalType, + quantOffsetType: quantOffsetType, + nlsfInterpCoefQ2: nlsfInterpQ2, + ltpMemLength: ltpMemLength, + frameLength: frameLength, + subfrLength: subfrLength, + nbSubfr: subfrCount, + predictLPCOrder: order, + shapingLPCOrder: shapeLPCOrderLowComplex, + }) + e.frameCounter++ + + // Emit every field in the order the decoder reads it. + vadBit := uint32(0) + if active { + vadBit = 1 + } + e.rangeEncoder.EncodeSymbolLogP(1, vadBit) + e.rangeEncoder.EncodeSymbolLogP(1, 0) // LBRR flag (no redundancy) + + e.emitFrameType(signalType, quantOffsetType, active) + e.emitGainIndices(gainIndices, signalType, false) + e.emitNLSFIndices(index1, indices2, bandwidth, voiced) + e.rangeEncoder.EncodeSymbolWithICDF(icdfNormalizedLSFInterpolationIndex, uint32(nlsfInterpQ2)) //nolint:gosec // G115 + if voiced { + primaryLag := int(lagIndex) + peMinLagMS*fsKHz + contour := uint32(contourIndex) //nolint:gosec // G115: contour index is non-negative. + period := uint32(periodicityIndex) //nolint:gosec // G115: periodicity index is non-negative. + scale := uint32(ltpScaleIndex) //nolint:gosec // G115: scale index is 0..2. + e.encodePitchLags(primaryLag, contour, bandwidth, nanoseconds20Ms, true) + e.encodeLTPFilter(period, toUint32(filterIndices)) + e.encodeLTPScaling(scale) + } + e.rangeEncoder.EncodeSymbolWithICDF(icdfLinearCongruentialGeneratorSeed, seed) + e.encodePulses(signalType, quantOffsetType, pulses, frameLength) + + // Carry state to the next frame. + copy(e.xBuf, analysis[frameLength:frameLength+ltpMemLength]) + e.isPreviousFrameVoiced = voiced + e.firstFrameAfterReset = false +} + +// toUint32 converts codebook indices to the type the emitters expect. +func toUint32(indices []int8) []uint32 { + out := make([]uint32, len(indices)) + for i, v := range indices { + out[i] = uint32(v) //nolint:gosec // G115: indices are small non-negative. + } + + return out +} + +// emitFrameType codes the signal type and quantization offset (RFC 6716 +// Section 4.2.7.3). +func (e *Encoder) emitFrameType(signalType frameSignalType, quantOffsetType frameQuantizationOffsetType, vad bool) { + high := quantOffsetType == frameQuantizationOffsetTypeHigh + if !vad { + sym := uint32(0) + if high { + sym = 1 + } + e.rangeEncoder.EncodeSymbolWithICDF(icdfFrameTypeVADInactive, sym) + + return + } + + var sym uint32 + switch { + case signalType == frameSignalTypeUnvoiced && high: + sym = 1 + case signalType == frameSignalTypeVoiced && !high: + sym = 2 + case signalType == frameSignalTypeVoiced && high: + sym = 3 + } + e.rangeEncoder.EncodeSymbolWithICDF(icdfFrameTypeVADActive, sym) +} diff --git a/internal/silk/enc_frame_test.go b/internal/silk/enc_frame_test.go new file mode 100644 index 0000000..36bf41a --- /dev/null +++ b/internal/silk/enc_frame_test.go @@ -0,0 +1,59 @@ +// 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" +) + +// TestEncodeSILKFrameDecodable encodes a SILK frame and decodes it with the +// real decoder, proving the bitstream is well-formed and self-consistent. This +// is the end-to-end gate short of opus_compare against libopus. +func TestEncodeSILKFrameDecodable(t *testing.T) { + for _, bandwidth := range []Bandwidth{BandwidthNarrowband, BandwidthMediumband, BandwidthWideband} { + fsKHz := silkInternalRate(bandwidth) + frameLength := 20 * fsKHz + + input := make([]int16, frameLength) + for i := range input { + input[i] = int16(4000*math.Sin(2*math.Pi*float64(i)/50) + 1500*math.Sin(2*math.Pi*float64(i)/13)) + } + + enc := NewEncoder() + enc.rangeEncoder.Init() + enc.encodeSILKFrame(input, bandwidth) + data := enc.rangeEncoder.Done() + require.NotEmpty(t, data) + + dec := NewDecoder() + out := make([]float32, frameLength) + err := dec.Decode(data, out, false, nanoseconds20Ms, bandwidth) + require.NoErrorf(t, err, "bandwidth %d", bandwidth) + + var energy float64 + for _, v := range out { + energy += float64(v) * float64(v) + } + assert.Positivef(t, energy, "bandwidth %d: decoded output is silent", bandwidth) + } +} + +// TestEncodeSILKFrameSilence checks a silent input encodes and decodes cleanly. +func TestEncodeSILKFrameSilence(t *testing.T) { + bandwidth := BandwidthWideband + frameLength := 20 * silkInternalRate(bandwidth) + + enc := NewEncoder() + enc.rangeEncoder.Init() + enc.encodeSILKFrame(make([]int16, frameLength), bandwidth) + data := enc.rangeEncoder.Done() + + dec := NewDecoder() + out := make([]float32, frameLength) + require.NoError(t, dec.Decode(data, out, false, nanoseconds20Ms, bandwidth)) +} diff --git a/internal/silk/encode_fidelity_test.go b/internal/silk/encode_fidelity_test.go new file mode 100644 index 0000000..faa8506 --- /dev/null +++ b/internal/silk/encode_fidelity_test.go @@ -0,0 +1,100 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "math" + "testing" + + "github.com/stretchr/testify/require" +) + +// voicedTestSignal builds a multi-frame speech-like signal: a pitch fundamental +// with a few harmonics (so the frame is classified voiced and the LTP path is +// exercised) plus a slow formant-like amplitude modulation. +func voicedTestSignal(fsKHz, frames int) []int16 { + fs := float64(fsKHz * 1000) + n := frames * 20 * fsKHz + out := make([]int16, n) + pitch := 140.0 + for i := range out { + t := float64(i) / fs + s := math.Sin(2*math.Pi*pitch*t) + + 0.5*math.Sin(2*math.Pi*2*pitch*t) + + 0.33*math.Sin(2*math.Pi*3*pitch*t) + + 0.2*math.Sin(2*math.Pi*5*pitch*t) + env := 0.6 + 0.4*math.Sin(2*math.Pi*3*t) + out[i] = int16(6000 * env * s / 2.03) + } + + return out +} + +// alignedSNR aligns b to a over a small lag range and returns the best +// energy-weighted SNR in dB over the aligned region. +func alignedSNR(a, b []int16, maxLag int) float64 { + best := math.Inf(-1) + lo := 20 * len(a) / 100 + hi := 80 * len(a) / 100 + for lag := range maxLag { + var sig, noise float64 + for i := lo; i < hi; i++ { + if i+lag >= len(b) { + break + } + d := float64(a[i]) - float64(b[i+lag]) + sig += float64(a[i]) * float64(a[i]) + noise += d * d + } + if noise == 0 { + return math.Inf(1) + } + if snr := 10 * math.Log10(sig/noise); snr > best { + best = snr + } + } + + return best +} + +// TestEncodeDecodeFidelity encodes a multi-frame voiced signal frame-by-frame, +// decodes it back and asserts the reconstruction quality per bandwidth. It +// locks in the analysis/quantization pipeline (voiced/LTP, noise shaping, NLSF +// interpolation) end to end. +func TestEncodeDecodeFidelity(t *testing.T) { + const frames = 25 + cases := []struct { + bandwidth Bandwidth + minSNR float64 + }{ + {BandwidthNarrowband, 12.0}, + {BandwidthMediumband, 12.0}, + {BandwidthWideband, 12.0}, + } + for _, tc := range cases { + fsKHz := silkInternalRate(tc.bandwidth) + frameLen := 20 * fsKHz + input := voicedTestSignal(fsKHz, frames) + + enc := NewEncoder() + dec := NewDecoder() + decoded := make([]int16, 0, len(input)) + for f := range frames { + frame := input[f*frameLen : (f+1)*frameLen] + payload := enc.Encode(frame, tc.bandwidth, 24000) + + out := make([]float32, frameLen) + require.NoErrorf(t, dec.Decode(payload, out, false, nanoseconds20Ms, tc.bandwidth), + "bandwidth %d frame %d", tc.bandwidth, f) + for _, v := range out { + s := math.Round(float64(v) * 32768) + decoded = append(decoded, int16(math.Max(-32768, math.Min(32767, s)))) + } + } + + snr := alignedSNR(input, decoded, 16) + t.Logf("bandwidth %d: SNR %.1f dB", tc.bandwidth, snr) + require.Greaterf(t, snr, tc.minSNR, "bandwidth %d reconstruction SNR too low", tc.bandwidth) + } +} diff --git a/internal/silk/encoder.go b/internal/silk/encoder.go new file mode 100644 index 0000000..ccbab70 --- /dev/null +++ b/internal/silk/encoder.go @@ -0,0 +1,71 @@ +// 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 + frameCounter int + targetBitrate int // target bitrate in bps (drives control_SNR) + packetLossPerc int // expected packet loss %, drives LTP state scaling + sumLogGainQ7 int32 // cumulative LTP gain limit (quant_LTP_gains) + xBuf []float32 // previous frame, as LTP-memory history for pitch analysis + 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/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..d3e7d21 --- /dev/null +++ b/internal/silk/find_pitch_lags_test.go @@ -0,0 +1,62 @@ +// 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, _, res, _ := 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") + } +} + +// 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..975d108 --- /dev/null +++ b/internal/silk/find_pred_coefs.go @@ -0,0 +1,121 @@ +// 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] + invGain := invGains[k] + 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)) + 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/gains.go b/internal/silk/gains.go new file mode 100644 index 0000000..a7e31a0 --- /dev/null +++ b/internal/silk/gains.go @@ -0,0 +1,109 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// Gain quantization constants from silk/define.h and the derived macros at the +// top of silk/gain_quant.c. +const ( + gainOffsetQ7 = 2090 // OFFSET = (MIN_QGAIN_DB*128)/6 + 16*128 + gainScaleQ16 = 2251 // SCALE_Q16 = (65536*(N_LEVELS_QGAIN-1)) / (((MAX-MIN)*128)/6) + gainInvScaleQ16 = 1907825 // INV_SCALE_Q16 = (65536*(((MAX-MIN)*128)/6)) / (N_LEVELS_QGAIN-1) == 0x1D1C71 + gainNLevels = 64 // N_LEVELS_QGAIN + gainMaxDelta = 36 // MAX_DELTA_GAIN_QUANT + gainMinDelta = -4 // MIN_DELTA_GAIN_QUANT + gainMaxLogQ7 = 3967 // 31 in Q7, the upper clamp used by silk_log2lin() +) + +// 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, _ := e.quantizeGains(gainsTargetQ16, subframeCount, conditional) + e.emitGainIndices(indices, signalType, conditional) + + return gainQ16 +} + +// quantizeGains is silk_gains_quant: it turns target gains into transmit +// indices and the dequantized gains, updating the running previousLogGain. For +// the independent first subframe the index is the full 6-bit gain index; for +// delta subframes it is the non-negative transmit index (0..40). +func (e *Encoder) quantizeGains( + gainsTargetQ16 []int32, + subframeCount int, + conditional bool, +) (indices []int8, gainQ16 []float32, gainQ16Int []int32) { + indices = make([]int8, subframeCount) + gainQ16 = make([]float32, subframeCount) + gainQ16Int = make([]int32, subframeCount) + + for subframeIndex := range subframeCount { + ind := smulwb(gainScaleQ16, lin2log(gainsTargetQ16[subframeIndex])-gainOffsetQ7) + if ind < e.previousLogGain { + ind++ + } + ind = clamp(0, ind, gainNLevels-1) + + if subframeIndex == 0 && !conditional { //nolint:nestif // faithful port of silk_gains_quant. + ind = clamp(e.previousLogGain+gainMinDelta, ind, gainNLevels-1) + e.previousLogGain = ind + indices[subframeIndex] = int8(ind) //nolint:gosec // G115: ind is in [0,63]. + } else { + delta := ind - e.previousLogGain + doubleStepThreshold := 2*gainMaxDelta - gainNLevels + e.previousLogGain + if delta > doubleStepThreshold { + delta = doubleStepThreshold + ((delta - doubleStepThreshold + 1) >> 1) + } + delta = clamp(gainMinDelta, delta, gainMaxDelta) + if delta > doubleStepThreshold { + e.previousLogGain += (delta << 1) - doubleStepThreshold + if e.previousLogGain > gainNLevels-1 { + e.previousLogGain = gainNLevels - 1 + } + } else { + e.previousLogGain += delta + } + indices[subframeIndex] = int8(delta - gainMinDelta) //nolint:gosec // G115: delta is in [gainMinDelta,gainMaxDelta]. + } + + inLogQ7 := (gainInvScaleQ16 * e.previousLogGain >> 16) + gainOffsetQ7 + inLogQ7 = min(inLogQ7, gainMaxLogQ7) + i := inLogQ7 >> 7 + f := inLogQ7 & 127 + gain := (1 << i) + ((-174*f*(128-f)>>16)+f)*((1<>7) + gainQ16Int[subframeIndex] = gain + gainQ16[subframeIndex] = float32(gain) + } + + e.haveEncoded = true + + return indices, gainQ16, gainQ16Int +} + +// 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 new file mode 100644 index 0000000..49b3e63 --- /dev/null +++ b/internal/silk/gains_test.go @@ -0,0 +1,129 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// TestEncodeSubframeGainsRoundTrip checks that encoded gains decode back to the +// same values through the decoder, with the range coder and gain state in sync. +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") + }) + } +} + +// TestLin2LogInverseOfLog2Lin 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). +func TestLin2LogRoundTripAgainstGainDequant(t *testing.T) { + // For every gain index 0..63, the log domain value fed to log2lin is + // exactly recoverable, so lin2log(log2lin(x)) must return x for the gain + // grid points used by the codec. + for logGain := range int32(gainNLevels) { + inLogQ7 := (gainInvScaleQ16 * logGain >> 16) + gainOffsetQ7 + inLogQ7 = min(inLogQ7, gainMaxLogQ7) + i := inLogQ7 >> 7 + f := inLogQ7 & 127 + gainQ16 := (1 << i) + ((-174*f*(128-f)>>16)+f)*((1<>7) + + got := lin2log(gainQ16) + // lin2log/log2lin are approximate inverses; the reference keeps the + // round-trip within a couple of Q7 units. + diff := got - inLogQ7 + if diff < 0 { + diff = -diff + } + assert.LessOrEqualf(t, diff, int32(3), "lin2log(log2lin(%d))=%d want ~%d", logGain, got, inLogQ7) + } +} diff --git a/internal/silk/lin2log.go b/internal/silk/lin2log.go new file mode 100644 index 0000000..128bc2d --- /dev/null +++ b/internal/silk/lin2log.go @@ -0,0 +1,47 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// Fixed-point log helpers needed by the encoder. silk_log2lin() is not here +// because it is spelled out inline wherever gains are dequantized. + +// ror32 rotates a 32-bit value right by rot bits; a negative rot rotates left. +func ror32(a32 int32, rot int) int32 { + x := uint32(a32) //nolint:gosec // G115: bit-level rotate, sign is irrelevant. + switch { + case rot == 0: + return a32 + case rot < 0: + m := uint(-rot) + + return int32((x << m) | (x >> (32 - m))) //nolint:gosec // G115 + default: + r := uint(rot) + + return int32((x << (32 - r)) | (x >> r)) //nolint:gosec // G115 + } +} + +// clzFrac returns the leading-zero count of in and the 7 bits after the +// leading one. +func clzFrac(in int32) (lz, fracQ7 int32) { + lz = int32(clz32(in)) //nolint:gosec // G115: leading-zero count is in [0,32]. + fracQ7 = ror32(in, 24-int(lz)) & 0x7f + + return lz, fracQ7 +} + +// smlawb returns a + smulwb(b, c). +func smlawb(a, b, c int32) int32 { + return a + smulwb(b, c) +} + +// lin2log approximates 128*log2(inLin) with a piecewise-parabolic fit; it is +// the near-inverse of the inline log2lin used for gains. +func lin2log(inLin int32) int32 { + lz, fracQ7 := clzFrac(inLin) + + // silk_ADD_LSHIFT32(silk_SMLAWB(frac_Q7, silk_MUL(frac_Q7, 128-frac_Q7), 179), 31-lz, 7) + return smlawb(fracQ7, fracQ7*(128-fracQ7), 179) + ((31 - lz) << 7) +} diff --git a/internal/silk/lpc_analysis.go b/internal/silk/lpc_analysis.go new file mode 100644 index 0000000..79143b1 --- /dev/null +++ b/internal/silk/lpc_analysis.go @@ -0,0 +1,142 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +const findLPCCondFac = 1e-5 + +// findLPC estimates LPC coefficients with Burg and converts them to NLSFs +// (silk_find_LPC_FLP). The optional NLSF-interpolation search is deferred; the +// interpolation index stays at 4 (no interpolation). x holds nb_subfr +// subframes of subfrLength samples each (including order preceding samples). +func findLPC(nlsfQ15 []int16, x []float32, minInvGain float32, subfrLength, nbSubfr, order int) { + a := make([]float32, order) + burgModifiedFLP(a, x, minInvGain, subfrLength, nbSubfr, order) + a2nlsfFLP(nlsfQ15, a, order) +} + +// burgModifiedFLP estimates LPC coefficients with Burg's method, stacked over +// nb_subfr subframes, limiting the prediction gain to 1/minInvGain. It writes +// the order coefficients into a and returns the residual energy +// (silk_burg_modified_FLP). +// +//nolint:gocognit,gocyclo,cyclop,maintidx,varnamelen // faithful port of silk_burg_modified_FLP. +func burgModifiedFLP(a, x []float32, minInvGain float32, subfrLength, nbSubfr, order int) float32 { + var cFirstRow, cLastRow [maxLPCOrder]float64 + var cAf, cAb [maxLPCOrder + 1]float64 + var af [maxLPCOrder]float64 + + c0 := energyFLP(x, nbSubfr*subfrLength) + for s := range nbSubfr { + xPtr := x[s*subfrLength:] + for n := 1; n < order+1; n++ { + cFirstRow[n-1] += innerProductFLP(xPtr, xPtr[n:], subfrLength-n) + } + } + copy(cLastRow[:], cFirstRow[:]) + + cAb[0] = c0 + findLPCCondFac*c0 + 1e-9 + cAf[0] = cAb[0] + invGain := 1.0 + reachedMaxGain := false + + for n := range order { + for s := range nbSubfr { + xPtr := x[s*subfrLength:] + tmp1 := float64(xPtr[n]) + tmp2 := float64(xPtr[subfrLength-n-1]) + for k := range n { + cFirstRow[k] -= float64(xPtr[n]) * float64(xPtr[n-k-1]) + cLastRow[k] -= float64(xPtr[subfrLength-n-1]) * float64(xPtr[subfrLength-n+k]) + atmp := af[k] + tmp1 += float64(xPtr[n-k-1]) * atmp + tmp2 += float64(xPtr[subfrLength-n+k]) * atmp + } + for k := 0; k <= n; k++ { + cAf[k] -= tmp1 * float64(xPtr[n-k]) + cAb[k] -= tmp2 * float64(xPtr[subfrLength-n+k-1]) + } + } + + tmp1 := cFirstRow[n] + tmp2 := cLastRow[n] + for k := range n { + atmp := af[k] + tmp1 += cLastRow[n-k-1] * atmp + tmp2 += cFirstRow[n-k-1] * atmp + } + cAf[n+1] = tmp1 //nolint:gosec // G602: n+1 <= order <= maxLPCOrder. + cAb[n+1] = tmp2 //nolint:gosec // G602: n+1 <= order <= maxLPCOrder. + + num := cAb[n+1] //nolint:gosec // G602: n+1 <= order <= maxLPCOrder. + nrgB := cAb[0] + nrgF := cAf[0] + for k := range n { + atmp := af[k] + num += cAb[n-k] * atmp + nrgB += cAb[k+1] * atmp + nrgF += cAf[k+1] * atmp + } + + rc := -2.0 * num / (nrgF + nrgB) + + gain := invGain * (1.0 - rc*rc) + if gain <= float64(minInvGain) { + rc = math.Sqrt(1.0 - float64(minInvGain)/invGain) + if num > 0 { + rc = -rc + } + invGain = float64(minInvGain) + reachedMaxGain = true + } else { + invGain = gain + } + + for k := range (n + 1) >> 1 { + t1 := af[k] + t2 := af[n-k-1] + af[k] = t1 + rc*t2 + af[n-k-1] = t2 + rc*t1 + } + af[n] = rc + + if reachedMaxGain { + for k := n + 1; k < order; k++ { + af[k] = 0.0 + } + + break + } + + for k := 0; k <= n+1; k++ { + t1 := cAf[k] + cAf[k] += rc * cAb[n-k+1] + cAb[n-k+1] += rc * t1 + } + } + + var residualEnergy float64 + if reachedMaxGain { + for k := range order { + a[k] = float32(-af[k]) + } + for s := range nbSubfr { + c0 -= energyFLP(x[s*subfrLength:], order) + } + residualEnergy = c0 * invGain + } else { + residualEnergy = cAf[0] + sumSq := 1.0 + for k := range order { + atmp := af[k] + residualEnergy += cAf[k+1] * atmp + sumSq += atmp * atmp + a[k] = float32(-atmp) + } + residualEnergy -= findLPCCondFac * c0 * sumSq + } + + return float32(residualEnergy) +} diff --git a/internal/silk/lpc_analysis_test.go b/internal/silk/lpc_analysis_test.go new file mode 100644 index 0000000..d0b801e --- /dev/null +++ b/internal/silk/lpc_analysis_test.go @@ -0,0 +1,82 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// arSignal generates a length-n AR(2) process x[i] = a1*x[i-1] + a2*x[i-2] + e. +func arSignal(n int, a1, a2 float32, seed uint32) []float32 { + x := make([]float32, n) + var prev1, prev2 float32 + state := seed + for i := range x { + state = 1664525*state + 1013904223 + e := float32(int32(state>>16)%2000-1000) / 1000 + cur := a1*prev1 + a2*prev2 + e + x[i] = cur + prev2, prev1 = prev1, cur + } + + return x +} + +func TestBurgModifiedFLPRecoversAR(t *testing.T) { + const ( + a1 = 0.5 + a2 = 0.2 + n = 512 + ) + x := arSignal(n, a1, a2, 42) + + a := make([]float32, 2) + resNrg := burgModifiedFLP(a, x, 1e-5, n, 1, 2) + + assert.InDelta(t, a1, a[0], 0.12) + assert.InDelta(t, a2, a[1], 0.12) + assert.Greater(t, resNrg, float32(0)) + assert.Less(t, resNrg, float32(energyFLP(x, n)), "AR signal should be predictable") +} + +func TestBurgModifiedFLPWhiteNoise(t *testing.T) { + const n = 512 + x := make([]float32, n) + state := uint32(777) + for i := range x { + state = 1664525*state + 1013904223 + x[i] = float32(int32(state>>16)%2000-1000) / 1000 + } + + a := make([]float32, 4) + burgModifiedFLP(a, x, 1e-5, n, 1, 4) + + // White noise has no predictable structure; coefficients should be small. + for k, c := range a { + require.InDeltaf(t, 0, c, 0.2, "coefficient %d", k) + } +} + +func TestFindLPC(t *testing.T) { + const ( + order = 16 + nbSubfr = 4 + subfrLength = 40 + order + ) + x := arSignal(nbSubfr*subfrLength, 0.6, 0.2, 99) + + nlsf := make([]int16, order) + findLPC(nlsf, x, 1e-4, subfrLength, nbSubfr, order) + + // A valid NLSF vector: non-decreasing and non-negative (Q15, <= int16 max). + for k := range nlsf { + require.GreaterOrEqual(t, nlsf[k], int16(0)) + } + for k := 1; k < order; k++ { + require.GreaterOrEqual(t, nlsf[k], nlsf[k-1]) + } +} 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) + } +} 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..e270f9b --- /dev/null +++ b/internal/silk/ltp_quant.go @@ -0,0 +1,208 @@ +// 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 ( + 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 ( + 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 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 { + 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_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..46b49ef --- /dev/null +++ b/internal/silk/ltp_scale_test.go @@ -0,0 +1,24 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" +) + +func TestLTPScaleControl(t *testing.T) { + const snrDBQ7 = 18 * 128 // ~18 dB + + // No packet loss: always minimum scaling (index 0, 15565). + idx, q14 := ltpScaleControl(20, snrDBQ7, 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, snrDBQ7, 25, 1, false) + assert.NotEqual(t, 0, idx, "high loss/gain: expected stronger scaling") + assert.Equal(t, ltpScalesTableQ14[idx], q14, "scale matches table entry") +} diff --git a/internal/silk/math_fix.go b/internal/silk/math_fix.go new file mode 100644 index 0000000..21b9262 --- /dev/null +++ b/internal/silk/math_fix.go @@ -0,0 +1,171 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// Additional fixed-point primitives from the RFC 6716 C macros, used by the +// analysis stages. + +// smulbb returns (int16)a * (int16)b. +func smulbb(a, b int32) int32 { + return int32(int16(a)) * int32(int16(b)) //nolint:gosec // G115: Q-format conversion. +} + +// smlabb returns a + (int16)b * (int16)c. +func smlabb(a, b, c int32) int32 { + return a + int32(int16(b))*int32(int16(c)) //nolint:gosec // G115: Q-format conversion. +} + +// smulww returns (a * b) >> 16 at 64-bit width. +func smulww(a, b int32) int32 { + return int32((int64(a) * int64(b)) >> 16) //nolint:gosec // G115: Q-format conversion. +} + +// smulwt returns (a * (b>>16)) >> 16. +func smulwt(a, b int32) int32 { + return int32((int64(a) * int64(b>>16)) >> 16) //nolint:gosec // G115: Q-format conversion. +} + +// smlawt returns a + smulwt(b, c). +func smlawt(a, b, c int32) int32 { + return a + smulwt(b, c) +} + +// sub32Ovflw / add32Ovflw are two's-complement wrapping arithmetic. +func sub32Ovflw(a, b int32) int32 { return int32(uint32(a) - uint32(b)) } //nolint:gosec // G115 +func add32Ovflw(a, b int32) int32 { return int32(uint32(a) + uint32(b)) } //nolint:gosec // G115 + +// addLShift32 returns a + (b << shift). +func addLShift32(a, b int32, shift uint) int32 { + return a + (b << shift) +} + +// addSat32 saturates the signed sum of a and b to int32. +func addSat32(a, b int32) int32 { + sum := int64(a) + int64(b) + if sum > math.MaxInt32 { + return math.MaxInt32 + } + if sum < math.MinInt32 { + return math.MinInt32 + } + + return int32(sum) +} + +// lshiftSat32 saturates a << shift to int32. +func lshiftSat32(a int32, shift uint) int32 { + return int32(clampI64(int64(a)< hi { + return hi + } + + return v +} + +// silkRand advances the LCG dither seed (silk_RAND). +func silkRand(seed int32) int32 { + return int32(uint32(907633515) + uint32(seed)*uint32(196314165)) //nolint:gosec // G115: wrapping LCG. +} + +// smlabbOvflw returns a + (int16)b * (int16)c with wrapping. +func smlabbOvflw(a, b, c int32) int32 { + return add32Ovflw(a, int32(int16(b))*int32(int16(c))) //nolint:gosec // G115: Q-format conversion. +} + +// div32VarQ approximates (a << qres) / b (silk_DIV32_varQ). +func div32VarQ(a, b int32, qres int) int32 { + aHeadrm := clz32(absInt32(a)) - 1 + a32Nrm := a << uint(aHeadrm) //nolint:gosec // G115: shift count is non-negative. + bHeadrm := clz32(absInt32(b)) - 1 + b32Nrm := b << uint(bHeadrm) //nolint:gosec // G115: shift count is non-negative. + b32Inv := (math.MaxInt32 >> 2) / (b32Nrm >> 16) + result := smulwb(a32Nrm, b32Inv) + a32Nrm = sub32Ovflw(a32Nrm, lshiftOvflw(smmul(b32Nrm, result), 3)) + result = smlawb(result, a32Nrm, b32Inv) + + lshift := 29 + aHeadrm - bHeadrm - qres + if lshift < 0 { + return lshiftSat32(result, uint(-lshift)) + } + if lshift < 32 { + return result >> uint(lshift) + } + + return 0 +} + +// sat16 saturates to the int16 range. +func sat16(a int32) int32 { + switch { + case a > math.MaxInt16: + return math.MaxInt16 + case a < math.MinInt16: + return math.MinInt16 + default: + return a + } +} + +// addPosSat32 saturates the sum of two non-negative values to int32. +func addPosSat32(a, b int32) int32 { + if (uint32(a)+uint32(b))&0x80000000 != 0 { //nolint:gosec // G115: bit test on the sum. + return math.MaxInt32 + } + + return a + b +} + +// sqrtApprox approximates the square root (silk_SQRT_APPROX). +func sqrtApprox(x int32) int32 { + if x <= 0 { + return 0 + } + lz, fracQ7 := clzFrac(x) + y := int32(46214) + if lz&1 != 0 { + y = 32768 + } + y >>= lz >> 1 + + return smlawb(y, y, smulbb(213, fracQ7)) +} + +//nolint:gochecknoglobals // fixed sigmoid lookup tables from sigm_Q15.c. +var ( + sigmLUTSlopeQ10 = [6]int32{237, 153, 73, 30, 12, 7} + sigmLUTPosQ15 = [6]int32{16384, 23955, 28861, 31213, 32178, 32548} + sigmLUTNegQ15 = [6]int32{16384, 8812, 3906, 1554, 589, 219} +) + +// sigmQ15 approximates the sigmoid in Q15 (silk_sigm_Q15). +func sigmQ15(inQ5 int32) int32 { + if inQ5 < 0 { + inQ5 = -inQ5 + if inQ5 >= 6*32 { + return 0 + } + ind := inQ5 >> 5 + + return sigmLUTNegQ15[ind] - smulbb(sigmLUTSlopeQ10[ind], inQ5&0x1F) + } + if inQ5 >= 6*32 { + return 32767 + } + ind := inQ5 >> 5 + + return sigmLUTPosQ15[ind] + smulbb(sigmLUTSlopeQ10[ind], inQ5&0x1F) +} diff --git a/internal/silk/math_fix_test.go b/internal/silk/math_fix_test.go new file mode 100644 index 0000000..c34d705 --- /dev/null +++ b/internal/silk/math_fix_test.go @@ -0,0 +1,45 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" +) + +func TestDiv32VarQ(t *testing.T) { + cases := []struct { + a, b int32 + q int + }{ + {1000, 7, 16}, + {-500000, 13, 10}, + {123456, 789, 10}, + {1 << 20, 3, 12}, + } + for _, c := range cases { + got := float64(div32VarQ(c.a, c.b, c.q)) + want := float64(c.a) * float64(int64(1)<>16)) >> 16 + assert.Equal(t, int32(int64(1000000)*int64(0x00030000>>16)>>16), smulwt(1000000, 0x00030000)) +} + +func TestAddSat32(t *testing.T) { + assert.Equal(t, int32(3), addSat32(1, 2)) + assert.Equal(t, int32(2147483647), addSat32(2147483647, 100)) + assert.Equal(t, int32(-2147483648), addSat32(-2147483648, -100)) +} + +func TestSilkRandDeterministic(t *testing.T) { + // Matches the decoder's LCG update: seed = 196314165*seed + 907633515. + seed := int32(12345) + want := int32(uint32(907633515) + uint32(seed)*uint32(196314165)) //nolint:gosec + assert.Equal(t, want, silkRand(seed)) +} diff --git a/internal/silk/nlsf_encode.go b/internal/silk/nlsf_encode.go new file mode 100644 index 0000000..e9ddb98 --- /dev/null +++ b/internal/silk/nlsf_encode.go @@ -0,0 +1,256 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "math" + "slices" +) + +const ( + nlsfQuantMaxAmplitude = 4 + nlsfQuantMaxAmplitudeExt = 10 + nlsfLevelAdjQ10 = 102 // round(NLSF_QUANT_LEVEL_ADJ * 1024) +) + +// nlsfStepSizes returns the Q16 quantization step and its Q6 inverse. +func nlsfStepSizes(bandwidth Bandwidth) (qstepQ16, invQstepQ6 int32) { + if bandwidth == BandwidthWideband { + return 9830, 427 + } + + return 11796, 356 +} + +// nlsfWeightQ9 computes the Q9 spectral weight for coefficient k from the +// stage-1 codebook vector, using the same formula as the decoder. +func nlsfWeightQ9(cb1 []uint, order, k int) int32 { + previous, next := uint(0), uint(256) + if k != 0 { + previous = cb1[k-1] + } + if k+1 != order { + next = cb1[k+1] + } + w2Q18 := (1024/(cb1[k]-previous) + 1024/(next-cb1[k])) << 16 + + i := ilog(int(w2Q18)) //nolint:gosec // G115 + f := int((w2Q18 >> (i - 8)) & 127) //nolint:gosec // G115 + y := 46214 + if i&1 != 0 { + y = 32768 + } + y >>= (32 - i) >> 1 + + return int32(int16(y + ((213 * f * y) >> 16))) //nolint:gosec // G115 +} + +// nlsfSecondOperand mirrors the decoder's residual reconstruction term for a +// stage-2 index. +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. +func quantizeNLSF(nlsfQ15 []int16, bandwidth Bandwidth) (int, []int8, []int16) { + order := len(nlsfQ15) + + cb1Set := codebookNormalizedLSFStageOneNarrowbandOrMediumband + predSelect := predictionWeightSelectionForNarrowbandAndMediumbandNormalizedLSF + predTable := predictionWeightForNarrowbandAndMediumbandNormalizedLSF + if bandwidth == BandwidthWideband { + cb1Set = codebookNormalizedLSFStageOneWideband + predSelect = predictionWeightSelectionForWidebandNormalizedLSF + predTable = predictionWeightForWidebandNormalizedLSF + } + qstepQ16, invQstepQ6 := nlsfStepSizes(bandwidth) + + bestIndex1 := 0 + bestIndices2 := make([]int8, order) + bestNLSF := make([]int16, order) + bestDistortion := int64(math.MaxInt64) + + indices2 := make([]int8, order) + weightsQ9 := make([]int32, order) + resReconQ10 := make([]int16, order) + candidate := make([]int16, order) + + for index1 := range cb1Set { + cb1 := cb1Set[index1] + for k := range order { + weightsQ9[k] = nlsfWeightQ9(cb1, order, k) + } + + // Greedily quantize the stage-2 residual backwards. + prevOut := int32(0) + for k := order - 1; k >= 0; k-- { //nolint:varnamelen // k indexes the coefficient, as in the C reference. + target := (int32(nlsfQ15[k]) - int32(cb1[k])<<7) * weightsQ9[k] >> 14 //nolint:gosec // G115 + predQ10 := int32(0) + if k+1 < order { + predQ10 = (int32(predTable[predSelect[index1][k]][k]) * prevOut) >> 8 //nolint:gosec // G115 + } + ind := clamp(-nlsfQuantMaxAmplitudeExt, (invQstepQ6*(target-predQ10))>>16, nlsfQuantMaxAmplitudeExt-1) + out0 := nlsfSecondOperand(ind, qstepQ16) + predQ10 + out1 := nlsfSecondOperand(ind+1, qstepQ16) + predQ10 + chosen, chosenOut := ind, out0 + if absInt32(target-out1) < absInt32(target-out0) { + chosen, chosenOut = ind+1, out1 + } + indices2[k] = int8(chosen) //nolint:gosec // G115 + resReconQ10[k] = int16(chosenOut) //nolint:gosec // G115 + prevOut = int32(resReconQ10[k]) + } + + // Reconstruct and stabilize as the decoder will, then score. + for k := range order { + candidate[k] = int16(clamp(0, //nolint:gosec // G115 + int32((int(cb1[k])<<7)+(int(resReconQ10[k])<<14)/int(weightsQ9[k])), 32767)) //nolint:gosec // G115 + } + stabilizeNLSF(candidate, order, bandwidth) + + var distortion int64 + for k := range order { + diff := int64(nlsfQ15[k]) - int64(candidate[k]) + distortion += int64(weightsQ9[k]) * diff * diff + } + if distortion < bestDistortion { + bestDistortion = distortion + bestIndex1 = index1 + copy(bestIndices2, indices2) + copy(bestNLSF, candidate) + } + } + + 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). +// +//nolint:cyclop // faithful port of silk_NLSF_stabilize. +func stabilizeNLSF(nlsfQ15 []int16, dLPC int, bandwidth Bandwidth) { + NDeltaMinQ15 := codebookMinimumSpacingForNormalizedLSCoefficientsNarrowbandAndMediumband + if bandwidth == BandwidthWideband { + NDeltaMinQ15 = codebookMinimumSpacingForNormalizedLSCoefficientsWideband + } + + for adjustment := 0; adjustment <= 19; adjustment++ { + i := 0 + iValue := int(math.MaxInt) + for nlsfIndex := 0; nlsfIndex <= len(nlsfQ15); nlsfIndex++ { + previousNLSF := 0 + currentNLSF := 32768 + if nlsfIndex != 0 { + previousNLSF = int(nlsfQ15[nlsfIndex-1]) + } + if nlsfIndex != len(nlsfQ15) { + currentNLSF = int(nlsfQ15[nlsfIndex]) + } + spacingValue := currentNLSF - previousNLSF - NDeltaMinQ15[nlsfIndex] + if spacingValue < iValue { + i = nlsfIndex + iValue = spacingValue + } + } + + switch { + case iValue >= 0: + return + case i == 0: + nlsfQ15[0] = int16(NDeltaMinQ15[0]) //nolint:gosec // G115 + + continue + case i == dLPC: + nlsfQ15[dLPC-1] = int16(32768 - NDeltaMinQ15[dLPC]) //nolint:gosec // G115 + + continue + } + + minCenterQ15 := NDeltaMinQ15[i] >> 1 + for k := 0; k <= i-1; k++ { + minCenterQ15 += NDeltaMinQ15[k] + } + maxCenterQ15 := 32768 - (NDeltaMinQ15[i] >> 1) + for k := i + 1; k <= dLPC; k++ { + maxCenterQ15 -= NDeltaMinQ15[k] + } + centerFreqQ15 := int(clamp( + int32(minCenterQ15), //nolint:gosec // G115 + int32((int(nlsfQ15[i-1])+int(nlsfQ15[i])+1)>>1), //nolint:gosec // G115 + int32(maxCenterQ15)), //nolint:gosec // G115 + ) + nlsfQ15[i-1] = int16(centerFreqQ15 - NDeltaMinQ15[i]>>1) //nolint:gosec // G115 + nlsfQ15[i] = nlsfQ15[i-1] + int16(NDeltaMinQ15[i]) //nolint:gosec // G115 + } + + slices.Sort(nlsfQ15) + for k := 0; k <= dLPC-1; k++ { + prevNLSF := int16(0) + if k != 0 { + prevNLSF = nlsfQ15[k-1] + } + nlsfQ15[k] = maxInt16(nlsfQ15[k], saturatingAddInt16(prevNLSF, int16(NDeltaMinQ15[k]))) //nolint:gosec // G115 + } + for k := dLPC - 1; k >= 0; k-- { + nextNLSF := 32768 + if k != dLPC-1 { + nextNLSF = int(nlsfQ15[k+1]) + } + nlsfQ15[k] = minInt16(nlsfQ15[k], int16(nextNLSF-NDeltaMinQ15[k+1])) //nolint:gosec // G115 + } +} diff --git a/internal/silk/nlsf_encode_test.go b/internal/silk/nlsf_encode_test.go new file mode 100644 index 0000000..1a03540 --- /dev/null +++ b/internal/silk/nlsf_encode_test.go @@ -0,0 +1,167 @@ +// 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" +) + +// genNLSF builds a deterministic, strictly increasing NLSF vector in Q15. +func genNLSF(order int, seed uint32) []int16 { + nlsf := make([]int16, order) + state := seed + cur := int32(400) + for k := range nlsf { + state = 1664525*state + 1013904223 + cur += int32(state>>20)%2000 + 700 + if cur > 32200 { + cur = 32200 + } + nlsf[k] = int16(cur) + } + + return nlsf +} + +// TestEncodeNLSFRoundTrip encodes an NLSF vector and decodes it back through +// the decoder, asserting the reconstructed vector and range coder match and +// the result is a valid (stabilized, increasing) NLSF vector. +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 +// that expose an out-of-range stage-2 index if the quantizer fails to clamp. +func genClusteredNLSF(order int, seed uint32) []int16 { + nlsf := make([]int16, order) + state := seed + state = 1664525*state + 1013904223 + center := int32(2000) + int32(state>>18)%28000 + cur := center + for k := range nlsf { + state = 1664525*state + 1013904223 + cur += int32(state>>28) % 3 // 0..2 Q15: far below the minimum spacing. + if cur > 32200 { + cur = 32200 + } + nlsf[k] = int16(cur) + } + + return nlsf +} + +// TestEncodeNLSFClusteredStaysInRange is a regression test for an inverted +// clamp in quantizeNLSF that let stage-2 indices exceed the extension range. +// The oversized index encodes an extension symbol the decoder cannot reproduce, +// so the range coder silently desynchronizes rather than failing loudly. Tightly +// clustered NLSF vectors trigger it most reliably. See pion/opus#147. +func TestEncodeNLSFClusteredStaysInRange(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 64 { + name := fmt.Sprintf("bw%d_voiced%t_seed%d", tc.bandwidth, voiced, seed) + t.Run(name, func(t *testing.T) { + input := genClusteredNLSF(tc.order, uint32(seed*131+tc.order+boolSeed(voiced))) //nolint:gosec // G115 + + // Stage-2 indices must stay within the extension range the + // decoder can reconstruct: [-Ext, Ext-1]. + stabilized := append([]int16(nil), input...) + stabilizeNLSF(stabilized, len(stabilized), tc.bandwidth) + _, indices2, _ := quantizeNLSF(stabilized, tc.bandwidth) + // The decoder reconstructs indices in [-Ext, Ext] inclusive + // (RFC 6716 4.2.7.5.2). quantizeNLSF clamps ind to [-Ext, Ext-1] + // and may then pick ind+1, so the emitted index reaches +Ext. + for k, v := range indices2 { + require.GreaterOrEqualf(t, int(v), -nlsfQuantMaxAmplitudeExt, "stage-2 index %d underflow", k) + require.LessOrEqualf(t, int(v), nlsfQuantMaxAmplitudeExt, "stage-2 index %d overflow", k) + } + + // And the full encode->decode round-trip must stay in sync. + 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") + }) + } + } + } +} diff --git a/internal/silk/noise_shape.go b/internal/silk/noise_shape.go new file mode 100644 index 0000000..d505fc2 --- /dev/null +++ b/internal/silk/noise_shape.go @@ -0,0 +1,41 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// Noise-shaping analysis primitives. The main silk_noise_shape_analysis_FLP +// orchestration (which populates the encoder control struct) is wired up with +// the rest of the frame encoder. + +const maxShapeLPCOrder = 24 + +// sigmoid returns 1/(1+exp(-x)) (silk_sigmoid). +func sigmoid(x float32) float32 { + return float32(1.0 / (1.0 + math.Exp(float64(-x)))) +} + +// warpedAutocorrelationFLP correlates the input after a chain of first-order +// all-pass warping sections (silk_warped_autocorrelation_FLP). order must be +// even; corr receives order+1 taps. warping 0 gives the regular autocorrelation. +func warpedAutocorrelationFLP(corr, input []float32, warping float32, length, order int) { + var state, c [maxShapeLPCOrder + 1]float64 //nolint:varnamelen // c is the correlation accumulator. + w := float64(warping) + for n := range length { + tmp1 := float64(input[n]) + for i := 0; i < order; i += 2 { + tmp2 := state[i] + w*state[i+1] - w*tmp1 //nolint:gosec // G602: order is even and < maxShapeLPCOrder. + state[i] = tmp1 + c[i] += state[0] * tmp1 + tmp1 = state[i+1] + w*state[i+2] - w*tmp2 //nolint:gosec // G602: order is even and < maxShapeLPCOrder. + state[i+1] = tmp2 //nolint:gosec // G602: order is even and < maxShapeLPCOrder. + c[i+1] += state[0] * tmp2 //nolint:gosec // G602: order is even and < maxShapeLPCOrder. + } + state[order] = tmp1 + c[order] += state[0] * tmp1 + } + for i := 0; i < order+1; i++ { + corr[i] = float32(c[i]) //nolint:gosec // G602: i <= order. + } +} diff --git a/internal/silk/noise_shape_analysis.go b/internal/silk/noise_shape_analysis.go new file mode 100644 index 0000000..efa27a4 --- /dev/null +++ b/internal/silk/noise_shape_analysis.go @@ -0,0 +1,308 @@ +// 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. +// +//nolint:cyclop // faithful port of the noise-shape analysis stage. +func (e *Encoder) noiseShapeAnalysis( + shapeBuf []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(shapeBuf[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 := e.quantizeGains(gainsTargetQ16, 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_test.go b/internal/silk/noise_shape_test.go new file mode 100644 index 0000000..19ae264 --- /dev/null +++ b/internal/silk/noise_shape_test.go @@ -0,0 +1,43 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import ( + "testing" + + "github.com/stretchr/testify/assert" +) + +func TestSigmoid(t *testing.T) { + assert.InDelta(t, 0.5, sigmoid(0), 1e-6) + assert.Greater(t, sigmoid(2), sigmoid(0)) + assert.Less(t, sigmoid(-2), sigmoid(0)) + assert.InDelta(t, 1, sigmoid(20), 1e-6) + assert.InDelta(t, 0, sigmoid(-20), 1e-6) +} + +// TestWarpedAutocorrelationZeroWarping checks that with no warping the result +// matches the regular autocorrelation. +func TestWarpedAutocorrelationZeroWarping(t *testing.T) { + const ( + length = 64 + order = 8 + ) + x := make([]float32, length) + state := uint32(7) + for i := range x { + state = 1664525*state + 1013904223 + x[i] = float32(int32(state>>16)%2000-1000) / 1000 + } + + warped := make([]float32, order+1) + warpedAutocorrelationFLP(warped, x, 0, length, order) + + regular := make([]float32, order+1) + autocorrelationFLP(regular, x, length, order+1) + + for i := range warped { + assert.InDeltaf(t, regular[i], warped[i], 1e-3, "tap %d", i) + } +} diff --git a/internal/silk/nsq.go b/internal/silk/nsq.go new file mode 100644 index 0000000..556ae39 --- /dev/null +++ b/internal/silk/nsq.go @@ -0,0 +1,359 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// Noise-shaping quantizer (silk_NSQ_c, NSQ.c). It applies short- and long-term +// prediction plus noise shaping to the scaled input and quantizes the +// excitation to integer pulses. This is the non-delayed-decision variant; +// NSQ_del_dec is a later quality refinement. + +const ( + maxFSKHz = 16 + maxSubFrameLength = 5 * maxFSKHz // 80 + maxFrameLength = 20 * maxFSKHz // 320 + nsqLPCBufLength = maxLPCOrder // 16 + harmShapeFIRTaps = 3 + quantLevelAdjustQ10 = 80 +) + +// quantizationOffsetsQ10[signalType>>1][quantOffsetType] (silk_Quantization_Offsets_Q10). +// +//nolint:gochecknoglobals +var quantizationOffsetsQ10 = [2][2]int32{{100, 240}, {32, 100}} + +// nsqState is the persistent noise-shaping quantizer state (silk_nsq_state). +type nsqState struct { + xq []int16 // quantized output with history, 2*maxFrameLength + sLTPShpQ14 []int32 // 2*maxFrameLength + sLPCQ14 []int32 // maxSubFrameLength + nsqLPCBufLength + sAR2Q14 [maxShapeLPCOrder]int32 + sLFARShpQ14 int32 + sDiffShpQ14 int32 + lagPrev int + sLTPBufIdx int + sLTPShpBufIdx int + randSeed int32 + prevGainQ16 int32 + rewhiteFlag int +} + +func newNSQState() *nsqState { + return &nsqState{ + xq: make([]int16, 2*maxFrameLength), + sLTPShpQ14: make([]int32, 2*maxFrameLength), + sLPCQ14: make([]int32, maxSubFrameLength+nsqLPCBufLength), + prevGainQ16: 65536, + } +} + +// nsqParams bundles the per-frame quantizer inputs from the encoder control. +type nsqParams struct { + predCoefQ12 []int16 // 2*maxLPCOrder + ltpCoefQ14 []int16 // ltpOrder*nbSubfr + arQ13 []int16 // nbSubfr*maxShapeLPCOrder + harmShapeGainQ14 []int32 + tiltQ14 []int32 + lfShpQ14 []int32 + gainsQ16 []int32 + pitchL []int + lambdaQ10 int32 + ltpScaleQ14 int32 + seed int32 + signalType frameSignalType + quantOffsetType frameQuantizationOffsetType + nlsfInterpCoefQ2 int + ltpMemLength int + frameLength int + subfrLength int + nbSubfr int + predictLPCOrder int + shapingLPCOrder int +} + +// quantize runs the NSQ over one frame, filling pulses with the quantized +// excitation indices. +// +//nolint:varnamelen // p is the NSQ parameter block, as in the C reference. +func (nsq *nsqState) quantize(x16 []int16, pulses []int8, p *nsqParams) { + nsq.randSeed = p.seed + lag := nsq.lagPrev + + sigIndex := 0 + if p.signalType == frameSignalTypeVoiced { + sigIndex = 1 + } + offsetIndex := 0 + if p.quantOffsetType == frameQuantizationOffsetTypeHigh { + offsetIndex = 1 + } + offsetQ10 := quantizationOffsetsQ10[sigIndex][offsetIndex] + + lsfInterpFlag := 1 + if p.nlsfInterpCoefQ2 == 4 { + lsfInterpFlag = 0 + } + + sLTPQ15 := make([]int32, p.ltpMemLength+p.frameLength) + sLTP := make([]int16, p.ltpMemLength+p.frameLength) + xScQ10 := make([]int32, p.subfrLength) + + nsq.sLTPShpBufIdx = p.ltpMemLength + nsq.sLTPBufIdx = p.ltpMemLength + pxqIndex := p.ltpMemLength + + for k := range p.nbSubfr { //nolint:varnamelen // k indexes the subframe. + aQ12 := p.predCoefQ12[((k>>1)|(1-lsfInterpFlag))*maxLPCOrder:] + bQ14 := p.ltpCoefQ14[k*ltpOrder:] + arShpQ13 := p.arQ13[k*maxShapeLPCOrder:] + + harmShapeFIRPackedQ14 := p.harmShapeGainQ14[k] >> 2 + harmShapeFIRPackedQ14 |= (p.harmShapeGainQ14[k] >> 1) << 16 + + nsq.rewhiteFlag = 0 + if p.signalType == frameSignalTypeVoiced { + lag = p.pitchL[k] + if k&(3-(lsfInterpFlag<<1)) == 0 { + startIdx := p.ltpMemLength - lag - p.predictLPCOrder - ltpOrder/2 + lpcAnalysisFilterFixed(sLTP[startIdx:], nsq.xq[startIdx+k*p.subfrLength:], + aQ12, p.ltpMemLength-startIdx, p.predictLPCOrder) + nsq.rewhiteFlag = 1 + nsq.sLTPBufIdx = p.ltpMemLength + } + } + + nsq.scaleStates(x16[k*p.subfrLength:], xScQ10, sLTP, sLTPQ15, k, lag, p) + nsq.noiseShapeQuantizer(xScQ10, pulses[k*p.subfrLength:], pxqIndex, sLTPQ15, + aQ12, bQ14, arShpQ13, lag, harmShapeFIRPackedQ14, p.tiltQ14[k], p.lfShpQ14[k], + p.gainsQ16[k], p.lambdaQ10, offsetQ10, p) + + pxqIndex += p.subfrLength + } + + nsq.lagPrev = p.pitchL[p.nbSubfr-1] + + // Shift the quantized-speech and shaping buffers, keeping the LTP memory. + copy(nsq.xq, nsq.xq[p.frameLength:p.frameLength+p.ltpMemLength]) + copy(nsq.sLTPShpQ14, nsq.sLTPShpQ14[p.frameLength:p.frameLength+p.ltpMemLength]) +} + +//nolint:gocyclo,cyclop,varnamelen // faithful port of the dense NSQ inner loop. +func (nsq *nsqState) noiseShapeQuantizer( + xScQ10 []int32, pulses []int8, xqIndex int, sLTPQ15 []int32, + aQ12, bQ14, arShpQ13 []int16, lag int, + harmShapeFIRPackedQ14, tiltQ14, lfShpQ14, gainQ16, lambdaQ10, offsetQ10 int32, + p *nsqParams, +) { + shpLagPtr := nsq.sLTPShpBufIdx - lag + harmShapeFIRTaps/2 + predLagPtr := nsq.sLTPBufIdx - lag + ltpOrder/2 + gainQ10 := gainQ16 >> 6 + psLPCIndex := nsqLPCBufLength - 1 + + for i := range p.subfrLength { + nsq.randSeed = silkRand(nsq.randSeed) + lpcPredQ10 := nsqShortPrediction(nsq.sLPCQ14, psLPCIndex, aQ12, p.predictLPCOrder) + + var ltpPredQ13 int32 + if p.signalType == frameSignalTypeVoiced { + ltpPredQ13 = 2 + ltpPredQ13 = smlawb(ltpPredQ13, sLTPQ15[predLagPtr], int32(bQ14[0])) + ltpPredQ13 = smlawb(ltpPredQ13, sLTPQ15[predLagPtr-1], int32(bQ14[1])) + ltpPredQ13 = smlawb(ltpPredQ13, sLTPQ15[predLagPtr-2], int32(bQ14[2])) + ltpPredQ13 = smlawb(ltpPredQ13, sLTPQ15[predLagPtr-3], int32(bQ14[3])) + ltpPredQ13 = smlawb(ltpPredQ13, sLTPQ15[predLagPtr-4], int32(bQ14[4])) + predLagPtr++ + } + + nARQ12 := nsqNoiseShapeFeedbackLoop(nsq.sDiffShpQ14, nsq.sAR2Q14[:], arShpQ13, p.shapingLPCOrder) + nARQ12 = smlawb(nARQ12, nsq.sLFARShpQ14, tiltQ14) + + nLFQ12 := smulwb(nsq.sLTPShpQ14[nsq.sLTPShpBufIdx-1], lfShpQ14) + nLFQ12 = smlawt(nLFQ12, nsq.sLFARShpQ14, lfShpQ14) + + tmp1 := sub32Ovflw(lpcPredQ10<<2, nARQ12) + tmp1 = sub32Ovflw(tmp1, nLFQ12) + if lag > 0 { + nLTPQ13 := smulwb(addSat32(nsq.sLTPShpQ14[shpLagPtr], nsq.sLTPShpQ14[shpLagPtr-2]), harmShapeFIRPackedQ14) + nLTPQ13 = smlawt(nLTPQ13, nsq.sLTPShpQ14[shpLagPtr-1], harmShapeFIRPackedQ14) + nLTPQ13 <<= 1 + shpLagPtr++ + tmp2 := ltpPredQ13 - nLTPQ13 + tmp1 = add32Ovflw(tmp2, lshiftOvflw(tmp1, 1)) + tmp1 = rshiftRound32(tmp1, 3) + } else { + tmp1 = rshiftRound32(tmp1, 2) + } + + rQ10 := xScQ10[i] - tmp1 + if nsq.randSeed < 0 { + rQ10 = -rQ10 + } + rQ10 = clamp(-(31 << 10), rQ10, 30<<10) + + q1Q10 := rQ10 - offsetQ10 + q1Q0 := q1Q10 >> 10 + if lambdaQ10 > 2048 { + rdoOffset := lambdaQ10/2 - 512 + switch { + case q1Q10 > rdoOffset: + q1Q0 = (q1Q10 - rdoOffset) >> 10 + case q1Q10 < -rdoOffset: + q1Q0 = (q1Q10 + rdoOffset) >> 10 + case q1Q10 < 0: + q1Q0 = -1 + default: + q1Q0 = 0 + } + } + + var q2Q10, rd1Q20, rd2Q20 int32 + switch { + case q1Q0 > 0: + q1Q10 = (q1Q0 << 10) - quantLevelAdjustQ10 + offsetQ10 + q2Q10 = q1Q10 + 1024 + rd1Q20 = smulbb(q1Q10, lambdaQ10) + rd2Q20 = smulbb(q2Q10, lambdaQ10) + case q1Q0 == 0: + q1Q10 = offsetQ10 + q2Q10 = q1Q10 + 1024 - quantLevelAdjustQ10 + rd1Q20 = smulbb(q1Q10, lambdaQ10) + rd2Q20 = smulbb(q2Q10, lambdaQ10) + case q1Q0 == -1: + q2Q10 = offsetQ10 + q1Q10 = q2Q10 - (1024 - quantLevelAdjustQ10) + rd1Q20 = smulbb(-q1Q10, lambdaQ10) + rd2Q20 = smulbb(q2Q10, lambdaQ10) + default: + q1Q10 = (q1Q0 << 10) + quantLevelAdjustQ10 + offsetQ10 + q2Q10 = q1Q10 + 1024 + rd1Q20 = smulbb(-q1Q10, lambdaQ10) + rd2Q20 = smulbb(-q2Q10, lambdaQ10) + } + rrQ10 := rQ10 - q1Q10 + rd1Q20 = smlabb(rd1Q20, rrQ10, rrQ10) + rrQ10 = rQ10 - q2Q10 + rd2Q20 = smlabb(rd2Q20, rrQ10, rrQ10) + if rd2Q20 < rd1Q20 { + q1Q10 = q2Q10 + } + + pulses[i] = int8(rshiftRound32(q1Q10, 10)) //nolint:gosec // G115 + + excQ14 := q1Q10 << 4 + if nsq.randSeed < 0 { + excQ14 = -excQ14 + } + + lpcExcQ14 := addLShift32(excQ14, ltpPredQ13, 1) + xqQ14 := add32Ovflw(lpcExcQ14, lpcPredQ10<<4) + + nsq.xq[xqIndex+i] = int16(sat16(rshiftRound32(smulww(xqQ14, gainQ10), 8))) //nolint:gosec // G115 + + psLPCIndex++ + nsq.sLPCQ14[psLPCIndex] = xqQ14 + nsq.sDiffShpQ14 = sub32Ovflw(xqQ14, lshiftOvflw(xScQ10[i], 4)) + sLFARShpQ14 := sub32Ovflw(nsq.sDiffShpQ14, lshiftOvflw(nARQ12, 2)) + nsq.sLFARShpQ14 = sLFARShpQ14 + nsq.sLTPShpQ14[nsq.sLTPShpBufIdx] = sub32Ovflw(sLFARShpQ14, lshiftOvflw(nLFQ12, 2)) + sLTPQ15[nsq.sLTPBufIdx] = lpcExcQ14 << 1 + nsq.sLTPShpBufIdx++ + nsq.sLTPBufIdx++ + nsq.randSeed = add32Ovflw(nsq.randSeed, int32(pulses[i])) + } + + copy(nsq.sLPCQ14, nsq.sLPCQ14[p.subfrLength:p.subfrLength+nsqLPCBufLength]) +} + +// scaleStates scales the input and LTP state by 1/gain and adjusts for gain +// changes (silk_nsq_scale_states). +// +//nolint:cyclop,varnamelen // faithful port of silk_nsq_scale_states. +func (nsq *nsqState) scaleStates( + x16 []int16, xScQ10 []int32, sLTP []int16, sLTPQ15 []int32, subfr, lag int, p *nsqParams, +) { + invGainQ31 := inverse32VarQ(max(p.gainsQ16[subfr], 1), 47) + invGainQ26 := rshiftRound32(invGainQ31, 5) + for i := range p.subfrLength { + xScQ10[i] = smulww(int32(x16[i]), invGainQ26) + } + + if nsq.rewhiteFlag != 0 { + if subfr == 0 { + invGainQ31 = smulwb(invGainQ31, p.ltpScaleQ14) << 2 + } + for i := nsq.sLTPBufIdx - lag - ltpOrder/2; i < nsq.sLTPBufIdx; i++ { + sLTPQ15[i] = smulwb(invGainQ31, int32(sLTP[i])) + } + } + + if p.gainsQ16[subfr] != nsq.prevGainQ16 { + gainAdjQ16 := div32VarQ(nsq.prevGainQ16, p.gainsQ16[subfr], 16) + for i := nsq.sLTPShpBufIdx - p.ltpMemLength; i < nsq.sLTPShpBufIdx; i++ { + nsq.sLTPShpQ14[i] = smulww(gainAdjQ16, nsq.sLTPShpQ14[i]) + } + if p.signalType == frameSignalTypeVoiced && nsq.rewhiteFlag == 0 { + for i := nsq.sLTPBufIdx - lag - ltpOrder/2; i < nsq.sLTPBufIdx; i++ { + sLTPQ15[i] = smulww(gainAdjQ16, sLTPQ15[i]) + } + } + nsq.sLFARShpQ14 = smulww(gainAdjQ16, nsq.sLFARShpQ14) + nsq.sDiffShpQ14 = smulww(gainAdjQ16, nsq.sDiffShpQ14) + for i := range nsqLPCBufLength { + nsq.sLPCQ14[i] = smulww(gainAdjQ16, nsq.sLPCQ14[i]) + } + for i := range maxShapeLPCOrder { + nsq.sAR2Q14[i] = smulww(gainAdjQ16, nsq.sAR2Q14[i]) + } + nsq.prevGainQ16 = p.gainsQ16[subfr] + } +} + +// nsqShortPrediction returns the short-term LPC prediction +// (silk_noise_shape_quantizer_short_prediction). +func nsqShortPrediction(buf []int32, bufIdx int, coef []int16, order int) int32 { + out := int32(order >> 1) //nolint:gosec // G115 + for i := range order { + out = smlawb(out, buf[bufIdx-i], int32(coef[i])) + } + + return out +} + +// nsqNoiseShapeFeedbackLoop applies the noise-shaping AR filter, updating its +// state (silk_NSQ_noise_shape_feedback_loop). +func nsqNoiseShapeFeedbackLoop(data0 int32, data1 []int32, coef []int16, order int) int32 { + tmp2 := data0 + tmp1 := data1[0] + data1[0] = tmp2 + + out := int32(order >> 1) //nolint:gosec // G115 + out = smlawb(out, tmp2, int32(coef[0])) + for j := 2; j < order; j += 2 { + tmp2 = data1[j-1] + data1[j-1] = tmp1 + out = smlawb(out, tmp1, int32(coef[j-1])) + tmp1 = data1[j] + data1[j] = tmp2 + out = smlawb(out, tmp2, int32(coef[j])) + } + data1[order-1] = tmp1 + out = smlawb(out, tmp1, int32(coef[order-1])) + + return out << 1 +} + +// lpcAnalysisFilterFixed is the fixed-point LPC residual filter used for +// re-whitening (silk_LPC_analysis_filter). b holds Q12 coefficients. +func lpcAnalysisFilterFixed(out, in []int16, b []int16, length, order int) { + for ix := order; ix < length; ix++ { + out32Q12 := smulbb(int32(in[ix-1]), int32(b[0])) + for j := 1; j < order; j++ { + out32Q12 = smlabbOvflw(out32Q12, int32(in[ix-1-j]), int32(b[j])) + } + out32Q12 = sub32Ovflw(int32(in[ix])<<12, out32Q12) + out[ix] = int16(sat16(rshiftRound32(out32Q12, 12))) //nolint:gosec // G115 + } + for j := range order { + out[j] = 0 + } +} diff --git a/internal/silk/nsq_test.go b/internal/silk/nsq_test.go new file mode 100644 index 0000000..ae9323d --- /dev/null +++ b/internal/silk/nsq_test.go @@ -0,0 +1,74 @@ +// 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" +) + +// TestNSQUnvoicedSmoke exercises the NSQ end to end for an unvoiced frame: +// it must run without panicking and produce a bounded, non-trivial pulse train. +// Bit-exact correctness is validated later end-to-end via opus_compare. +func TestNSQUnvoicedSmoke(t *testing.T) { + const ( + fsKHz = 16 + nbSubfr = 4 + subfrLength = 5 * fsKHz + frameLength = nbSubfr * subfrLength + ltpMemLength = 20 * fsKHz + predictLPCOrder = 16 + shapingLPCOrder = 16 + ) + + x16 := make([]int16, frameLength) + for i := range x16 { + x16[i] = int16(3000 * math.Sin(2*math.Pi*float64(i)/37)) + } + pulses := make([]int8, frameLength) + + gains := make([]int32, nbSubfr) + for k := range gains { + gains[k] = 100 * 65536 // gain ~ signal scale so the excitation is non-degenerate + } + + p := &nsqParams{ //nolint:varnamelen // p is the NSQ parameter block. + predCoefQ12: make([]int16, 2*maxLPCOrder), + ltpCoefQ14: make([]int16, ltpOrder*nbSubfr), + arQ13: make([]int16, nbSubfr*maxShapeLPCOrder), + harmShapeGainQ14: make([]int32, nbSubfr), + tiltQ14: make([]int32, nbSubfr), + lfShpQ14: make([]int32, nbSubfr), + gainsQ16: gains, + pitchL: make([]int, nbSubfr), + lambdaQ10: 1024, + ltpScaleQ14: 15565, + seed: 1, + signalType: frameSignalTypeUnvoiced, + quantOffsetType: frameQuantizationOffsetTypeLow, + nlsfInterpCoefQ2: 4, + ltpMemLength: ltpMemLength, + frameLength: frameLength, + subfrLength: subfrLength, + nbSubfr: nbSubfr, + predictLPCOrder: predictLPCOrder, + shapingLPCOrder: shapingLPCOrder, + } + + nsq := newNSQState() + require.NotPanics(t, func() { + nsq.quantize(x16, pulses, p) + }) + + nonZero := 0 + for _, v := range pulses { + if v != 0 { + nonZero++ + } + } + assert.Positive(t, nonZero, "expected a non-trivial pulse train") +} diff --git a/internal/silk/pitch_analysis.go b/internal/silk/pitch_analysis.go new file mode 100644 index 0000000..25f89c7 --- /dev/null +++ b/internal/silk/pitch_analysis.go @@ -0,0 +1,399 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// stage3Array holds the per-subframe, per-codebook, per-lag correlations or +// energies computed for the stage-3 refinement. +type stage3Array = [peMaxNBSubfr][peNBCbksStage3Max][peNBStage3Lags]float32 + +// silkLog2 approximates log2 (silk_log2). +func silkLog2(x float64) float32 { + return float32(3.32192809488736 * math.Log10(x)) +} + +// stage3Params selects the codebook/range tables for the stage-3 search. +func stage3Params(nbSubfr, complexity int) (lagRange [][2]int8, lagCB [][]int8, nbCbkSearch int) { + if nbSubfr == peMaxNBSubfr { + rng := silkLagRangeStage3[complexity] + lagRange = rng[:] + lagCB = silkCBLagsStage3[:] + + return lagRange, lagCB, int(silkNbCbkSearchsStage3[complexity]) + } + rng := silkLagRangeStage3_10ms + lagRange = rng[:] + lagCB = silkCBLagsStage3_10ms[:] + + return lagRange, lagCB, peNBCbksStage3_10ms +} + +// calcCorrST3 fills the stage-3 cross-correlation array +// (silk_P_Ana_calc_corr_st3). +func calcCorrST3(crossCorr *stage3Array, frame []float32, startLag, sfLength, nbSubfr, complexity int) { + lagRange, lagCB, nbCbkSearch := stage3Params(nbSubfr, complexity) + + var scratch [scratchSizePitch]float32 + var xcorr [scratchSizePitch]float32 + targetOffset := sfLength << 2 + for k := range nbSubfr { //nolint:varnamelen // k indexes the subframe. + lagLow := int(lagRange[k][0]) + lagHigh := int(lagRange[k][1]) + pitchXcorr(frame[targetOffset:], frame[targetOffset-startLag-lagHigh:], xcorr[:], sfLength, lagHigh-lagLow+1) + + lagCounter := 0 + for j := lagLow; j <= lagHigh; j++ { + scratch[lagCounter] = xcorr[lagHigh-j] + lagCounter++ + } + + delta := lagLow + for i := range nbCbkSearch { + idx := int(lagCB[k][i]) - delta + for j := range peNBStage3Lags { + crossCorr[k][i][j] = scratch[idx+j] + } + } + targetOffset += sfLength + } +} + +// calcEnergyST3 fills the stage-3 energy array (silk_P_Ana_calc_energy_st3). +func calcEnergyST3(energies *stage3Array, frame []float32, startLag, sfLength, nbSubfr, complexity int) { + lagRange, lagCB, nbCbkSearch := stage3Params(nbSubfr, complexity) + + var scratch [scratchSizePitch]float32 + targetOffset := sfLength << 2 + for k := range nbSubfr { //nolint:varnamelen // k indexes the subframe. + lagCounter := 0 + basisOffset := targetOffset - (startLag + int(lagRange[k][0])) + energy := energyFLP(frame[basisOffset:], sfLength) + 1e-3 + scratch[lagCounter] = float32(energy) + lagCounter++ + + lagDiff := int(lagRange[k][1]) - int(lagRange[k][0]) + 1 + for i := 1; i < lagDiff; i++ { + energy -= float64(frame[basisOffset+sfLength-i]) * float64(frame[basisOffset+sfLength-i]) + energy += float64(frame[basisOffset-i]) * float64(frame[basisOffset-i]) + scratch[lagCounter] = float32(energy) + lagCounter++ + } + + delta := int(lagRange[k][0]) + for i := range nbCbkSearch { + idx := int(lagCB[k][i]) - delta + for j := range peNBStage3Lags { + energies[k][i][j] = scratch[idx+j] + } + } + targetOffset += sfLength + } +} + +// pitchAnalysisCore estimates the pitch lag (silk_pitch_analysis_core_FLP). +// It returns the lag index, contour index, per-subframe lags in pitchOut, and +// whether the frame is voiced. ltpCorr is updated in place. Only 8 and 16 kHz +// are handled; 12 kHz (medium-band) is pending down2_3. +// +//nolint:gocognit,gocyclo,cyclop,maintidx,unparam // faithful port; complexity is fixed at 2 today. +func pitchAnalysisCore( + frame []float32, + pitchOut []int, + ltpCorr *float32, + prevLag int, + searchThres1, searchThres2 float32, + fsKHz, complexity, nbSubfr int, +) (lagIndex int16, contourIndex int8, voiced bool) { + frameLength := (peLTPMemLengthMS + nbSubfr*peSubfrLengthMS) * fsKHz + frameLength4kHz := (peLTPMemLengthMS + nbSubfr*peSubfrLengthMS) * 4 + frameLength8kHz := (peLTPMemLengthMS + nbSubfr*peSubfrLengthMS) * 8 + sfLength := peSubfrLengthMS * fsKHz + sfLength4kHz := peSubfrLengthMS * 4 + sfLength8kHz := peSubfrLengthMS * 8 + minLag := peMinLagMS * fsKHz + minLag4kHz := peMinLagMS * 4 + minLag8kHz := peMinLagMS * 8 + maxLag := peMaxLagMS*fsKHz - 1 + maxLag4kHz := peMaxLagMS * 4 + maxLag8kHz := peMaxLagMS*8 - 1 + + // Resample to 8 kHz, then decimate to 4 kHz. + frame8FIX := make([]int16, frameLength8kHz) + if fsKHz == 16 { + frame16FIX := make([]int16, frameLength) + float2ShortArray(frame16FIX, frame[:frameLength]) + var state [2]int32 + resamplerDown2(&state, frame8FIX, frame16FIX) + } else { + float2ShortArray(frame8FIX, frame[:frameLength8kHz]) + } + + frame4FIX := make([]int16, frameLength4kHz) + var state [2]int32 + resamplerDown2(&state, frame4FIX, frame8FIX) + + frame8kHz := make([]float32, frameLength8kHz) + frame4kHz := make([]float32, frameLength4kHz) + short2FloatArray(frame8kHz, frame8FIX) + short2FloatArray(frame4kHz, frame4FIX) + + // Low-pass filter (differentiator's inverse), int16-saturating. + for i := frameLength4kHz - 1; i > 0; i-- { + frame4kHz[i] = float32(sat16(int32(frame4kHz[i]) + int32(frame4kHz[i-1]))) + } + + var c [peMaxNBSubfr][(peMaxLag >> 1) + 5]float32 //nolint:varnamelen // c is the correlation grid. + xcorr := make([]float32, maxLag4kHz-minLag4kHz+1) + + // First stage at 4 kHz. + targetOffset := sfLength4kHz << 2 + for k := 0; k < nbSubfr>>1; k++ { + basisOffset := targetOffset - minLag4kHz + pitchXcorr( + frame4kHz[targetOffset:], frame4kHz[targetOffset-maxLag4kHz:], xcorr, sfLength8kHz, maxLag4kHz-minLag4kHz+1, + ) + + crossCorr := float64(xcorr[maxLag4kHz-minLag4kHz]) + normalizer := energyFLP(frame4kHz[targetOffset:], sfLength8kHz) + + energyFLP(frame4kHz[basisOffset:], sfLength8kHz) + + float64(sfLength8kHz)*4000.0 + c[0][minLag4kHz] += float32(2 * crossCorr / normalizer) + + for d := minLag4kHz + 1; d <= maxLag4kHz; d++ { + basisOffset-- + crossCorr = float64(xcorr[maxLag4kHz-d]) + normalizer += float64(frame4kHz[basisOffset])*float64(frame4kHz[basisOffset]) - + float64(frame4kHz[basisOffset+sfLength8kHz])*float64(frame4kHz[basisOffset+sfLength8kHz]) + c[0][d] += float32(2 * crossCorr / normalizer) + } + targetOffset += sfLength8kHz + } + + // Short-lag bias. + for i := maxLag4kHz; i >= minLag4kHz; i-- { + c[0][i] -= c[0][i] * float32(i) / 4096.0 + } + + lengthDSrch := 4 + 2*complexity + dSrch := make([]int, peDSrchLength) + insertionSortDecreasingFLP(c[0][minLag4kHz:], dSrch, maxLag4kHz-minLag4kHz+1, lengthDSrch) + + cmax := c[0][minLag4kHz] + if cmax < 0.2 { + clearPitch(pitchOut, nbSubfr) + *ltpCorr = 0 + + return 0, 0, false + } + + threshold := searchThres1 * cmax + for i := range lengthDSrch { + if c[0][minLag4kHz+i] > threshold { + dSrch[i] = (dSrch[i] + minLag4kHz) << 1 //nolint:gosec // G602: i < lengthDSrch. + } else { + lengthDSrch = i + + break + } + } + + dComp := make([]int16, (peMaxLag>>1)+5) + for i := 0; i < lengthDSrch; i++ { + dComp[dSrch[i]] = 1 + } + for i := maxLag8kHz + 3; i >= minLag8kHz; i-- { + dComp[i] += dComp[i-1] + dComp[i-2] + } + lengthDSrch = 0 + for i := minLag8kHz; i < maxLag8kHz+1; i++ { + if dComp[i+1] > 0 { //nolint:gosec // G602: i+1 <= maxLag8kHz+1 < len(dComp). + dSrch[lengthDSrch] = i + lengthDSrch++ + } + } + for i := maxLag8kHz + 3; i >= minLag8kHz; i-- { + dComp[i] += dComp[i-1] + dComp[i-2] + dComp[i-3] + } + lengthDComp := 0 + for i := minLag8kHz; i < maxLag8kHz+4; i++ { + if dComp[i] > 0 { //nolint:gosec // G602: i < maxLag8kHz+4 < len(dComp). + dComp[lengthDComp] = int16(i - 2) + lengthDComp++ + } + } + + // Second stage at 8 kHz. + for k := range c { + for d := range c[k] { //nolint:gosec // G602: k < len(c). + c[k][d] = 0 //nolint:gosec // G602: indices within c. + } + } + src8kHz := frame8kHz + if fsKHz == 8 { + src8kHz = frame + } + targetOffset = peLTPMemLengthMS * 8 + for k := range nbSubfr { + energyTmp := energyFLP(src8kHz[targetOffset:], sfLength8kHz) + 1.0 + for j := 0; j < lengthDComp; j++ { + d := int(dComp[j]) + basisOffset := targetOffset - d + crossCorr := innerProductFLP(src8kHz[basisOffset:], src8kHz[targetOffset:], sfLength8kHz) + if crossCorr > 0 { + energy := energyFLP(src8kHz[basisOffset:], sfLength8kHz) + c[k][d] = float32(2 * crossCorr / (energy + energyTmp)) + } else { + c[k][d] = 0 + } + } + targetOffset += sfLength8kHz + } + + // Search the lag range and codebook. + ccmax := float32(0) + ccmaxB := float32(-1000) + cbimax := 0 + lag := -1 + + var prevLagLog2 float32 + if prevLag > 0 { + if fsKHz == 16 { + prevLag >>= 1 + } + prevLagLog2 = silkLog2(float64(prevLag)) + } + + var lagCB [][]int8 + var nbCbkSearch int + if nbSubfr == peMaxNBSubfr { + lagCB = silkCBLagsStage2[:] + if fsKHz == 8 && complexity > 0 { + nbCbkSearch = peNBCbksStage2Ext + } else { + nbCbkSearch = peNBCbksStage2 + } + } else { + lagCB = silkCBLagsStage2_10ms[:] + nbCbkSearch = peNBCbksStage2_10ms + } + + cc := make([]float32, peNBCbksStage2Ext) + for k := 0; k < lengthDSrch; k++ { + d := dSrch[k] //nolint:varnamelen // d is the candidate lag, as in the C reference. + for j := range nbCbkSearch { + cc[j] = 0 + for i := range nbSubfr { + cc[j] += c[i][d+int(lagCB[i][j])] + } + } + ccmaxNew := float32(-1000) + cbimaxNew := 0 + for i := range nbCbkSearch { + if cc[i] > ccmaxNew { + ccmaxNew = cc[i] + cbimaxNew = i + } + } + + lagLog2 := silkLog2(float64(d)) + ccmaxNewB := ccmaxNew - peShortlagBias*float32(nbSubfr)*lagLog2 + if prevLag > 0 { + deltaLagLog2Sqr := lagLog2 - prevLagLog2 + deltaLagLog2Sqr *= deltaLagLog2Sqr + ccmaxNewB -= pePrevlagBias * float32(nbSubfr) * (*ltpCorr) * deltaLagLog2Sqr / (deltaLagLog2Sqr + 0.5) + } + + if ccmaxNewB > ccmaxB && ccmaxNew > float32(nbSubfr)*searchThres2 { + ccmaxB = ccmaxNewB + ccmax = ccmaxNew + lag = d + cbimax = cbimaxNew + } + } + + if lag == -1 { + clearPitch(pitchOut, nbSubfr) + *ltpCorr = 0 + + return 0, 0, false + } + + *ltpCorr = ccmax / float32(nbSubfr) + + if fsKHz > 8 { + lag <<= 1 // fsKHz == 16 + lag = int(clamp(int32(minLag), int32(lag), int32(maxLag))) //nolint:gosec // G115 + startLag := max(lag-2, minLag) + endLag := min(lag+2, maxLag) + lagNew := lag + cbimax = 0 + ccmax = -1000 + + var crossCorrST3, energiesST3 stage3Array + calcCorrST3(&crossCorrST3, frame, startLag, sfLength, nbSubfr, complexity) + calcEnergyST3(&energiesST3, frame, startLag, sfLength, nbSubfr, complexity) + + contourBias := float32(peFlatcontourBias) / float32(lag) + lagCB3, nbCbkSearch3, cbkSize3 := stage3Codebook(nbSubfr, complexity) + _ = cbkSize3 + + targetOffset := peLTPMemLengthMS * fsKHz + energyTmp := energyFLP(frame[targetOffset:], nbSubfr*sfLength) + 1.0 + lagCounter := 0 + for d := startLag; d <= endLag; d++ { //nolint:varnamelen // d is the candidate lag, as in the C reference. + for j := range nbCbkSearch3 { //nolint:varnamelen // j indexes the codebook. + crossCorr := 0.0 + energy := energyTmp + for k := range nbSubfr { + crossCorr += float64(crossCorrST3[k][j][lagCounter]) + energy += float64(energiesST3[k][j][lagCounter]) + } + var ccmaxNew float32 + if crossCorr > 0 { + ccmaxNew = float32(2 * crossCorr / energy) + ccmaxNew *= 1.0 - contourBias*float32(j) + } + if ccmaxNew > ccmax && d+int(silkCBLagsStage3[0][j]) <= maxLag { + ccmax = ccmaxNew + lagNew = d + cbimax = j + } + } + lagCounter++ + } + + for k := range nbSubfr { + pitchOut[k] = lagNew + int(lagCB3[k][cbimax]) + pitchOut[k] = int(clamp(int32(minLag), int32(pitchOut[k]), int32(peMaxLagMS*fsKHz))) //nolint:gosec // G115 + } + lagIndex = int16(lagNew - minLag) //nolint:gosec // G115 + contourIndex = int8(cbimax) + } else { + for k := range nbSubfr { + pitchOut[k] = lag + int(lagCB[k][cbimax]) + pitchOut[k] = int(clamp(int32(minLag8kHz), int32(pitchOut[k]), int32(peMaxLagMS*8))) + } + lagIndex = int16(lag - minLag8kHz) + contourIndex = int8(cbimax) + } + + return lagIndex, contourIndex, true +} + +// stage3Codebook returns the stage-3 lag codebook and search count. +func stage3Codebook(nbSubfr, complexity int) (lagCB [][]int8, nbCbkSearch, cbkSize int) { + if nbSubfr == peMaxNBSubfr { + return silkCBLagsStage3[:], int(silkNbCbkSearchsStage3[complexity]), peNBCbksStage3Max + } + + return silkCBLagsStage3_10ms[:], peNBCbksStage3_10ms, peNBCbksStage3_10ms +} + +// clearPitch zeroes the first nbSubfr pitch lags. +func clearPitch(pitchOut []int, nbSubfr int) { + for k := range nbSubfr { + pitchOut[k] = 0 + } +} diff --git a/internal/silk/pitch_analysis_test.go b/internal/silk/pitch_analysis_test.go new file mode 100644 index 0000000..564bda4 --- /dev/null +++ b/internal/silk/pitch_analysis_test.go @@ -0,0 +1,93 @@ +// 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 TestSilkLog2(t *testing.T) { + assert.InDelta(t, 3, silkLog2(8), 1e-4) + assert.InDelta(t, 0, silkLog2(1), 1e-4) +} + +// TestPitchAnalysisCoreDetectsPeriod feeds a periodic signal and checks the +// detected lag matches the period and the frame reads as voiced. +func TestPitchAnalysisCoreDetectsPeriod(t *testing.T) { + const ( + fsKHz = 16 + nbSubfr = 4 + period = 80 // 16 kHz / 80 = 200 Hz + ) + frameLength := (peLTPMemLengthMS + nbSubfr*peSubfrLengthMS) * fsKHz + frame := make([]float32, frameLength) + for i := range frame { + frame[i] = float32(5000 * math.Sin(2*math.Pi*float64(i)/period)) + } + + pitchOut := make([]int, nbSubfr) + ltpCorr := float32(0) + lagIndex, contourIndex, voiced := pitchAnalysisCore(frame, pitchOut, <pCorr, 0, 0.4, 0.3, fsKHz, 2, nbSubfr) + + require.True(t, voiced, "a clean periodic tone should be voiced") + for k, p := range pitchOut { + assert.InDeltaf(t, period, p, 8, "subframe %d lag", k) + } + assert.GreaterOrEqual(t, lagIndex, int16(0)) + assert.GreaterOrEqual(t, contourIndex, int8(0)) + assert.Greater(t, ltpCorr, float32(0)) +} + +// TestPitchAnalysisCoreEscapesOnNoise checks the low-correlation escape path +// runs without error and yields zeroed lags when unvoiced. +func TestPitchAnalysisCoreEscapesOnNoise(t *testing.T) { + const ( + fsKHz = 16 + nbSubfr = 4 + ) + frameLength := (peLTPMemLengthMS + nbSubfr*peSubfrLengthMS) * fsKHz + frame := make([]float32, frameLength) + state := uint32(987654321) + for i := range frame { + state = 1664525*state + 1013904223 + frame[i] = float32(int32(state>>16)%400 - 200) + } + + pitchOut := make([]int, nbSubfr) + ltpCorr := float32(0) + _, _, voiced := pitchAnalysisCore(frame, pitchOut, <pCorr, 0, 0.7, 0.6, fsKHz, 2, nbSubfr) + + if !voiced { + for _, p := range pitchOut { + assert.Equal(t, 0, p) + } + assert.Equal(t, float32(0), ltpCorr) + } +} + +func TestPitchAnalysisCore8kHz(t *testing.T) { + const ( + fsKHz = 8 + nbSubfr = 4 + period = 48 + ) + frameLength := (peLTPMemLengthMS + nbSubfr*peSubfrLengthMS) * fsKHz + frame := make([]float32, frameLength) + for i := range frame { + frame[i] = float32(5000 * math.Sin(2*math.Pi*float64(i)/period)) + } + + pitchOut := make([]int, nbSubfr) + ltpCorr := float32(0) + _, _, voiced := pitchAnalysisCore(frame, pitchOut, <pCorr, 0, 0.4, 0.3, fsKHz, 2, nbSubfr) + + require.True(t, voiced) + for k, p := range pitchOut { + assert.InDeltaf(t, period, p, 6, "subframe %d lag", k) + } +} 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/pitch_tables.go b/internal/silk/pitch_tables.go new file mode 100644 index 0000000..970d36b --- /dev/null +++ b/internal/silk/pitch_tables.go @@ -0,0 +1,69 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +// Pitch estimator constants and codebooks (pitch_est_defines.h, +// pitch_est_tables.c). +const ( + peMaxNBSubfr = 4 + peSubfrLengthMS = 5 + peLTPMemLengthMS = 20 + peMaxFSKHz = 16 + peMaxLagMS = 18 + peMinLagMS = 2 + peMaxLag = peMaxLagMS * peMaxFSKHz + peDSrchLength = 24 + peNBStage3Lags = 5 + peNBCbksStage2 = 3 + peNBCbksStage2Ext = 11 + peNBCbksStage2_10ms = 3 + peNBCbksStage3Max = 34 + peNBCbksStage3_10ms = 12 + silkPEMaxComplex = 2 + scratchSizePitch = 22 + + peShortlagBias = 0.2 + pePrevlagBias = 0.2 + peFlatcontourBias = 0.05 +) + +//nolint:gochecknoglobals // pitch codebook tables from pitch_est_tables.c. +var ( + silkCBLagsStage2_10ms = [peMaxNBSubfr >> 1][]int8{ + {0, 1, 0}, + {0, 0, 1}, + } + + silkCBLagsStage3_10ms = [peMaxNBSubfr >> 1][]int8{ + {0, 0, 1, -1, 1, -1, 2, -2, 2, -2, 3, -3}, + {0, 1, 0, 1, -1, 2, -1, 2, -2, 3, -2, 3}, + } + + silkLagRangeStage3_10ms = [peMaxNBSubfr >> 1][2]int8{ + {-3, 7}, + {-2, 7}, + } + + silkCBLagsStage2 = [peMaxNBSubfr][]int8{ + {0, 2, -1, -1, -1, 0, 0, 1, 1, 0, 1}, + {0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0}, + {0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0}, + {0, -1, 2, 1, 0, 1, 1, 0, 0, -1, -1}, + } + + silkCBLagsStage3 = [peMaxNBSubfr][]int8{ + {0, 0, 1, -1, 0, 1, -1, 0, -1, 1, -2, 2, -2, -2, 2, -3, 2, 3, -3, -4, 3, -4, 4, 4, -5, 5, -6, -5, 6, -7, 6, 5, 8, -9}, + {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 1, -1, 0, 1, -1, -1, 1, -1, 2, 1, -1, 2, -2, -2, 2, -2, 2, 2, 3, -3}, + {0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, -1, 1, 0, 0, 2, 1, -1, 2, -1, -1, 2, -1, 2, 2, -1, 3, -2, -2, -2, 3}, + {0, 1, 0, 0, 1, 0, 1, -1, 2, -1, 2, -1, 2, 3, -2, 3, -2, -2, 4, 4, -3, 5, -3, -4, 6, -4, 6, 5, -5, 8, -6, -5, -7, 9}, + } + + silkLagRangeStage3 = [silkPEMaxComplex + 1][peMaxNBSubfr][2]int8{ + {{-5, 8}, {-1, 6}, {-1, 6}, {-4, 10}}, + {{-6, 10}, {-2, 6}, {-1, 6}, {-5, 10}}, + {{-9, 12}, {-3, 7}, {-2, 7}, {-7, 13}}, + } + + silkNbCbkSearchsStage3 = [silkPEMaxComplex + 1]int8{16, 24, 34} +) diff --git a/internal/silk/pitch_util.go b/internal/silk/pitch_util.go new file mode 100644 index 0000000..74f6be1 --- /dev/null +++ b/internal/silk/pitch_util.go @@ -0,0 +1,130 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +// Helpers used by the pitch analysis core: 2:1 downsampler, float/short +// conversions, energy, cross-correlation, and index-tracking insertion sorts. + +const ( + resamplerDown2Coef0 = 9872 + resamplerDown2Coef1 = 39809 - 65536 // -25727 +) + +// resamplerDown2 halves the sample rate with a second-order all-pass filter +// (silk_resampler_down2). out must have len(in)/2 samples; state persists +// across calls. +func resamplerDown2(state *[2]int32, out, in []int16) { + for k := range len(in) >> 1 { //nolint:varnamelen // k indexes the sample pair. + in32 := int32(in[2*k]) << 10 + y := in32 - state[0] + x := smlawb(y, y, resamplerDown2Coef1) + out32 := state[0] + x + state[0] = in32 + x + + in32 = int32(in[2*k+1]) << 10 + y = in32 - state[1] + x = smulwb(y, resamplerDown2Coef0) + out32 = out32 + state[1] + x + state[1] = in32 + x + + out[k] = int16(sat16(rshiftRound32(out32, 11))) //nolint:gosec // G115 + } +} + +// float2ShortArray rounds and saturates float samples to int16. +func float2ShortArray(out []int16, in []float32) { + for i := range out { + out[i] = int16(sat16(int32(math.RoundToEven(float64(in[i]))))) //nolint:gosec // G115 + } +} + +// short2FloatArray widens int16 samples to float. +func short2FloatArray(out []float32, in []int16) { + for i := range out { + out[i] = float32(in[i]) + } +} + +// energyFLP returns the signal energy in double precision. +func energyFLP(data []float32, n int) float64 { + var result float64 + for i := range n { + result += float64(data[i]) * float64(data[i]) + } + + return result +} + +// pitchXcorr computes xcorr[j] = for j in [0, maxPitch). +func pitchXcorr(x, y, xcorr []float32, length, maxPitch int) { + for j := range maxPitch { + xcorr[j] = float32(innerProductFLP(x, y[j:], length)) + } +} + +// insertionSortDecreasingFLP sorts the first K of L values into decreasing +// order, tracking their original indices (silk_insertion_sort_decreasing_FLP). +// +//nolint:dupl,varnamelen // twin of the increasing-order sort below; l/k are lengths, as in the C reference. +func insertionSortDecreasingFLP(a []float32, idx []int, l, k int) { + for i := range k { + idx[i] = i + } + for i := 1; i < k; i++ { + value := a[i] + j := i - 1 + for ; j >= 0 && value > a[j]; j-- { + a[j+1] = a[j] + idx[j+1] = idx[j] + } + a[j+1] = value + idx[j+1] = i + } + for i := k; i < l; i++ { + value := a[i] + if value > a[k-1] { //nolint:gosec // G602: k-1 >= 0 for k >= 1. + j := k - 2 + for ; j >= 0 && value > a[j]; j-- { + a[j+1] = a[j] + idx[j+1] = idx[j] + } + a[j+1] = value + idx[j+1] = i + } + } +} + +// insertionSortIncreasing sorts the first K of L values into increasing order, +// tracking their original indices (silk_insertion_sort_increasing). +// +//nolint:dupl,varnamelen // twin of the decreasing-order sort above; l/k are lengths, as in the C reference. +func insertionSortIncreasing(a []int32, idx []int, l, k int) { + for i := range k { + idx[i] = i + } + for i := 1; i < k; i++ { + value := a[i] + j := i - 1 + for ; j >= 0 && value < a[j]; j-- { + a[j+1] = a[j] + idx[j+1] = idx[j] + } + a[j+1] = value + idx[j+1] = i + } + for i := k; i < l; i++ { + value := a[i] + if value < a[k-1] { //nolint:gosec // G602: k-1 >= 0 for k >= 1. + j := k - 2 + for ; j >= 0 && value < a[j]; j-- { + a[j+1] = a[j] + idx[j+1] = idx[j] + } + a[j+1] = value + idx[j+1] = i + } + } +} diff --git a/internal/silk/pitch_util_test.go b/internal/silk/pitch_util_test.go new file mode 100644 index 0000000..91da2e3 --- /dev/null +++ b/internal/silk/pitch_util_test.go @@ -0,0 +1,82 @@ +// 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 TestResamplerDown2(t *testing.T) { + in := make([]int16, 64) + for i := range in { + in[i] = 1000 + } + out := make([]int16, len(in)/2) + var state [2]int32 + resamplerDown2(&state, out, in) + + // A DC input should settle near the same DC level after the transient. + assert.InDelta(t, 1000, out[len(out)-1], 30) + + // Deterministic for identical state and input. + out2 := make([]int16, len(in)/2) + var state2 [2]int32 + resamplerDown2(&state2, out2, in) + assert.Equal(t, out, out2) +} + +func TestFloatShortConversion(t *testing.T) { + out16 := make([]int16, 4) + float2ShortArray(out16, []float32{0.4, 1.6, -2.5, 100000}) + assert.Equal(t, int16(0), out16[0]) + assert.Equal(t, int16(2), out16[1]) + assert.Equal(t, int16(-2), out16[2]) // round half to even + assert.Equal(t, int16(math.MaxInt16), out16[3]) + + outF := make([]float32, 2) + short2FloatArray(outF, []int16{-5, 7}) + assert.Equal(t, []float32{-5, 7}, outF) +} + +func TestEnergyFLP(t *testing.T) { + assert.InDelta(t, 1+4+9, energyFLP([]float32{1, 2, 3}, 3), 1e-9) +} + +func TestPitchXcorr(t *testing.T) { + // A period-4 signal correlates most strongly at lag 4. + length := 32 + sig := make([]float32, length+8) + for i := range sig { + sig[i] = float32(math.Sin(2 * math.Pi * float64(i) / 4)) + } + xcorr := make([]float32, 8) + pitchXcorr(sig, sig, xcorr, length, 8) + assert.Greater(t, xcorr[4], xcorr[1]) + assert.Greater(t, xcorr[4], xcorr[2]) + assert.Greater(t, xcorr[4], xcorr[3]) +} + +func TestInsertionSortDecreasingFLP(t *testing.T) { + a := []float32{3, 1, 4, 1, 5, 9, 2, 6} + idx := make([]int, 4) + insertionSortDecreasingFLP(a, idx, len(a), 4) + assert.Equal(t, []float32{9, 6, 5, 4}, a[:4]) + assert.Equal(t, 5, idx[0]) // 9 was at index 5 + assert.Equal(t, 7, idx[1]) // 6 was at index 7 + assert.Equal(t, 4, idx[2]) // 5 was at index 4 + assert.Equal(t, 2, idx[3]) // 4 was at index 2 +} + +func TestInsertionSortIncreasing(t *testing.T) { + a := []int32{3, 1, 4, 1, 5, 9, 2, 6} + idx := make([]int, 4) + insertionSortIncreasing(a, idx, len(a), 4) + require.Equal(t, []int32{1, 1, 2, 3}, a[:4]) + assert.Equal(t, 6, idx[2]) // value 2 was at index 6 + assert.Equal(t, 0, idx[3]) // value 3 was at index 0 +} 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") + }) + } + } + } + } +} diff --git a/internal/silk/vad.go b/internal/silk/vad.go new file mode 100644 index 0000000..bcef45e --- /dev/null +++ b/internal/silk/vad.go @@ -0,0 +1,266 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package silk + +import "math" + +const ( + vadNBands = 4 + vadInternalSubframesLog2 = 2 + vadInternalSubframes = 1 << vadInternalSubframesLog2 + vadNoiseLevelSmoothCoefQ16 = 1024 + vadNoiseLevelsBias = 50 + vadNegativeOffsetQ5 = 128 + vadSNRFactorQ16 = 45000 + vadSNRSmoothCoefQ18 = 4096 + anaFilterBankA20 = 5394 << 1 + anaFilterBankA21 = -24290 + variableHPMinCutoffHz = 60 + variableHPMaxCutoffHz = 100 + variableHPMinCutoffQ16 = variableHPMinCutoffHz << 16 + variableHPMaxDeltaFreqQ7 = 51 // round(0.4 * 128) + variableHPSmoothCoef1Q16 = 6554 // round(0.1 * 65536) + variableHPSmoothInitDefault = 0 +) + +// tiltWeights weights each band's SNR in the frequency-tilt measure. +// +//nolint:gochecknoglobals +var tiltWeights = [vadNBands]int32{30000, 6000, -12000, -12000} + +// vadState holds the per-channel voice-activity detector state +// (silk_VAD_state). +type vadState struct { + anaState [2]int32 + anaState1 [2]int32 + anaState2 [2]int32 + xnrgSubfr [vadNBands]int32 + nrgRatioSmthQ8 [vadNBands]int32 + hpState int16 + nl [vadNBands]int32 + invNL [vadNBands]int32 + noiseLevelBias [vadNBands]int32 + counter int32 +} + +// newVADState initializes the VAD with approximate pink-noise levels. +func newVADState() vadState { + var v vadState //nolint:varnamelen // v is the VAD state under construction. + for b := range vadNBands { + v.noiseLevelBias[b] = max(vadNoiseLevelsBias/int32(b+1), 1) + } + for b := range vadNBands { + v.nl[b] = 100 * v.noiseLevelBias[b] + v.invNL[b] = math.MaxInt32 / v.nl[b] + } + v.counter = 15 + for b := range vadNBands { + v.nrgRatioSmthQ8[b] = 100 * 256 + } + + return v +} + +// getSpeechActivityQ8 returns the Q8 speech activity for one frame and updates +// the VAD state (silk_VAD_GetSA_Q8). It also returns the frequency tilt and +// per-band input-quality measures used by later analysis stages. +// +//nolint:cyclop // faithful port of silk_VAD_GetSA_Q8. +func (v *vadState) getSpeechActivityQ8( + pcm []int16, frameLength, fsKHz int, +) (saQ8 int, tiltQ15 int32, qualityQ15 [vadNBands]int32) { + decimatedFramelength1 := frameLength >> 1 + decimatedFramelength2 := frameLength >> 2 + decimatedFramelength := frameLength >> 3 + + var xOffset [vadNBands]int + xOffset[1] = decimatedFramelength + decimatedFramelength2 + xOffset[2] = xOffset[1] + decimatedFramelength + xOffset[3] = xOffset[2] + decimatedFramelength2 + x := make([]int16, xOffset[3]+decimatedFramelength1) //nolint:varnamelen // x holds the analysis subbands. + + anaFiltBank1(pcm, &v.anaState, x, x[xOffset[3]:], frameLength) + anaFiltBank1(x, &v.anaState1, x, x[xOffset[2]:], decimatedFramelength1) + anaFiltBank1(x, &v.anaState2, x, x[xOffset[1]:], decimatedFramelength2) + + // High-pass differentiator on the lowest band. + x[decimatedFramelength-1] >>= 1 + hpStateTmp := x[decimatedFramelength-1] + for i := decimatedFramelength - 1; i > 0; i-- { + x[i-1] >>= 1 + x[i] -= x[i-1] + } + x[0] -= v.hpState + v.hpState = hpStateTmp + + // Energy per band. + var xnrg [vadNBands]int32 + for b := range vadNBands { //nolint:varnamelen // b indexes the frequency band. + bandLength := frameLength >> min(vadNBands-b, vadNBands-1) + subframeLength := bandLength >> vadInternalSubframesLog2 + offset := 0 + xnrg[b] = v.xnrgSubfr[b] + var sumSquared int32 + for s := range vadInternalSubframes { + sumSquared = 0 + for i := range subframeLength { + xTmp := int32(x[xOffset[b]+i+offset]) >> 3 + sumSquared = smlabb(sumSquared, xTmp, xTmp) + } + if s < vadInternalSubframes-1 { + xnrg[b] = addPosSat32(xnrg[b], sumSquared) + } else { + xnrg[b] = addPosSat32(xnrg[b], sumSquared>>1) + } + offset += subframeLength + } + v.xnrgSubfr[b] = sumSquared + } + + v.getNoiseLevels(&xnrg) + + // Signal-plus-noise to noise ratio. + var sumSquared, inputTilt int32 + var nrgToNoiseRatioQ8 [vadNBands]int32 + for b := range vadNBands { //nolint:varnamelen // b indexes the frequency band. + speechNrg := xnrg[b] - v.nl[b] + if speechNrg <= 0 { + nrgToNoiseRatioQ8[b] = 256 + + continue + } + if uint32(xnrg[b])&0xFF800000 == 0 { //nolint:gosec // G115: xnrg is non-negative here. + nrgToNoiseRatioQ8[b] = (xnrg[b] << 8) / (v.nl[b] + 1) + } else { + nrgToNoiseRatioQ8[b] = xnrg[b] / ((v.nl[b] >> 8) + 1) + } + snrQ7 := lin2log(nrgToNoiseRatioQ8[b]) - 8*128 + sumSquared = smlabb(sumSquared, snrQ7, snrQ7) + if speechNrg < 1<<20 { + snrQ7 = smulwb(sqrtApprox(speechNrg)<<6, snrQ7) + } + inputTilt = smlawb(inputTilt, tiltWeights[b], snrQ7) + } + + sumSquared /= vadNBands + pSNRdBQ7 := int32(int16(3 * sqrtApprox(sumSquared))) //nolint:gosec // G115 + saQ15 := sigmQ15(smulwb(vadSNRFactorQ16, pSNRdBQ7) - vadNegativeOffsetQ5) + tiltQ15 = (sigmQ15(inputTilt) - 16384) << 1 + + // Scale by power level. + var powerNrg int32 + for b := range vadNBands { + powerNrg += int32(b+1) * ((xnrg[b] - v.nl[b]) >> 4) + } + if frameLength == 20*fsKHz { + powerNrg >>= 1 + } + switch { + case powerNrg <= 0: + saQ15 >>= 1 + case powerNrg < 16384: + powerNrg = sqrtApprox(powerNrg << 16) + saQ15 = smulwb(32768+powerNrg, saQ15) + } + saQ8 = min(int(saQ15>>7), math.MaxUint8) + + // Smoothed energy-to-noise ratio and per-band quality. + smoothCoefQ16 := smulwb(vadSNRSmoothCoefQ18, smulwb(saQ15, saQ15)) + if frameLength == 10*fsKHz { + smoothCoefQ16 >>= 1 + } + for b := range vadNBands { + v.nrgRatioSmthQ8[b] = smlawb(v.nrgRatioSmthQ8[b], nrgToNoiseRatioQ8[b]-v.nrgRatioSmthQ8[b], smoothCoefQ16) + snrQ7 := 3 * (lin2log(v.nrgRatioSmthQ8[b]) - 8*128) + qualityQ15[b] = sigmQ15((snrQ7 - 16*128) >> 4) + } + + return saQ8, tiltQ15, qualityQ15 +} + +// getNoiseLevels updates the smoothed per-band noise estimate +// (silk_VAD_GetNoiseLevels). +func (v *vadState) getNoiseLevels(pX *[vadNBands]int32) { + minCoef := int32(0) + if v.counter < 1000 { + minCoef = math.MaxInt16 / ((v.counter >> 4) + 1) + v.counter++ + } + for k := range vadNBands { //nolint:varnamelen // k indexes the frequency band. + nl := v.nl[k] + nrg := addPosSat32(pX[k], v.noiseLevelBias[k]) + invNrg := math.MaxInt32 / nrg + + var coef int32 + switch { + case nrg > nl<<3: + coef = vadNoiseLevelSmoothCoefQ16 >> 3 + case nrg < nl: + coef = vadNoiseLevelSmoothCoefQ16 + default: + coef = smulwb(smulww(invNrg, nl), vadNoiseLevelSmoothCoefQ16<<1) + } + coef = max(coef, minCoef) + + v.invNL[k] = smlawb(v.invNL[k], invNrg-v.invNL[k], coef) + nl = min(math.MaxInt32/v.invNL[k], 0x00FFFFFF) + v.nl[k] = nl + } +} + +// anaFiltBank1 splits a signal into low and high bands (silk_ana_filt_bank_1). +func anaFiltBank1(in []int16, state *[2]int32, outL, outH []int16, n int) { + for k := range n >> 1 { //nolint:varnamelen // k indexes the sample pair. + in32 := int32(in[2*k]) << 10 + y := in32 - state[0] + xVal := smlawb(y, y, anaFilterBankA21) + out1 := state[0] + xVal + state[0] = in32 + xVal + + in32 = int32(in[2*k+1]) << 10 + y = in32 - state[1] + xVal = smulwb(y, anaFilterBankA20) + out2 := state[1] + xVal + state[1] = in32 + xVal + + outL[k] = int16(sat16(rshiftRound32(out2+out1, 11))) //nolint:gosec // G115 + outH[k] = int16(sat16(rshiftRound32(out2-out1, 11))) //nolint:gosec // G115 + } +} + +// hpVariableCutoff returns the updated high-pass smoother state, adapting the +// cutoff to the previous frame's pitch (silk_HP_variable_cutoff). It only +// changes when the previous frame was voiced. +func hpVariableCutoff( + smth1Q15 int32, + prevSignalType frameSignalType, + prevLag, fsKHz int, + qualityBand0Q15 int32, + speechActivityQ8 int, +) int32 { + if prevSignalType != frameSignalTypeVoiced { + return smth1Q15 + } + + pitchFreqHzQ16 := (int32(fsKHz*1000) << 16) / int32(prevLag) //nolint:gosec // G115 + pitchFreqLogQ7 := lin2log(pitchFreqHzQ16) - (16 << 7) + + pitchFreqLogQ7 = smlawb(pitchFreqLogQ7, + smulwb((-qualityBand0Q15)<<2, qualityBand0Q15), + pitchFreqLogQ7-(lin2log(variableHPMinCutoffQ16)-(16<<7))) + + deltaFreqQ7 := pitchFreqLogQ7 - (smth1Q15 >> 8) + if deltaFreqQ7 < 0 { + deltaFreqQ7 *= 3 + } + deltaFreqQ7 = clamp(-variableHPMaxDeltaFreqQ7, deltaFreqQ7, variableHPMaxDeltaFreqQ7) + + smth1Q15 = smlawb( + smth1Q15, + smulbb(int32(speechActivityQ8), deltaFreqQ7), //nolint:gosec // G115 + variableHPSmoothCoef1Q16, + ) + + return clamp(lin2log(variableHPMinCutoffHz)<<8, smth1Q15, lin2log(variableHPMaxCutoffHz)<<8) +} diff --git a/internal/silk/vad_test.go b/internal/silk/vad_test.go new file mode 100644 index 0000000..2c7b6b1 --- /dev/null +++ b/internal/silk/vad_test.go @@ -0,0 +1,102 @@ +// 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 TestSigmQ15(t *testing.T) { + assert.Equal(t, int32(16384), sigmQ15(0)) + assert.Equal(t, int32(32767), sigmQ15(1000)) + assert.Equal(t, int32(0), sigmQ15(-1000)) + assert.Greater(t, sigmQ15(32), sigmQ15(0)) + assert.Less(t, sigmQ15(-32), sigmQ15(0)) +} + +func TestSqrtApprox(t *testing.T) { + assert.Equal(t, int32(0), sqrtApprox(0)) + assert.Equal(t, int32(0), sqrtApprox(-5)) + // The approximation is documented as within ~10% for small outputs and + // ~2.5% for outputs above 120. + for _, x := range []int32{1000, 50000, 1 << 20, 1 << 28} { + got := float64(sqrtApprox(x)) + want := math.Sqrt(float64(x)) + assert.InEpsilonf(t, want, got, 0.11, "sqrtApprox(%d)", x) + } +} + +// tone builds a full-scale-ish sine at freq Hz for the given length. +func tone(length, freqHz, fsHz int, amp float64) []int16 { + pcm := make([]int16, length) + for i := range pcm { + pcm[i] = int16(amp * math.Sin(2*math.Pi*float64(freqHz)*float64(i)/float64(fsHz))) + } + + return pcm +} + +func TestVADSpeechActivityBounds(t *testing.T) { + const fsKHz = 16 + frameLength := 20 * fsKHz + + silence := make([]int16, frameLength) + loud := tone(frameLength, 300, fsKHz*1000, 8000) + + silentVAD := newVADState() + loudVAD := newVADState() + var silentSA, loudSA int + for range 10 { + var quality [vadNBands]int32 + silentSA, _, quality = silentVAD.getSpeechActivityQ8(silence, frameLength, fsKHz) + require.GreaterOrEqual(t, silentSA, 0) + require.LessOrEqual(t, silentSA, 255) + for _, q := range quality { + require.GreaterOrEqual(t, q, int32(0)) + require.LessOrEqual(t, q, int32(32767)) + } + + loudSA, _, _ = loudVAD.getSpeechActivityQ8(loud, frameLength, fsKHz) + require.GreaterOrEqual(t, loudSA, 0) + require.LessOrEqual(t, loudSA, 255) + } + + assert.Greater(t, loudSA, silentSA, "a loud tone should read as more active than silence") + assert.Less(t, silentSA, 64, "silence should read as low activity") +} + +func TestVADDeterministic(t *testing.T) { + const fsKHz = 16 + frameLength := 10 * fsKHz + pcm := tone(frameLength, 500, fsKHz*1000, 5000) + + a := newVADState() + b := newVADState() + for range 5 { + saA, tiltA, qualA := a.getSpeechActivityQ8(pcm, frameLength, fsKHz) + saB, tiltB, qualB := b.getSpeechActivityQ8(pcm, frameLength, fsKHz) + require.Equal(t, saA, saB) + require.Equal(t, tiltA, tiltB) + require.Equal(t, qualA, qualB) + } +} + +func TestHPVariableCutoff(t *testing.T) { + lo := lin2log(variableHPMinCutoffHz) << 8 + hi := lin2log(variableHPMaxCutoffHz) << 8 + initial := (lo + hi) / 2 + + // Unvoiced previous frame leaves the smoother untouched. + got := hpVariableCutoff(initial, frameSignalTypeUnvoiced, 120, 16, 20000, 200) + assert.Equal(t, initial, got) + + // Voiced previous frame keeps the smoother within the cutoff range. + got = hpVariableCutoff(initial, frameSignalTypeVoiced, 120, 16, 20000, 200) + assert.GreaterOrEqual(t, got, lo) + assert.LessOrEqual(t, got, hi) +} diff --git a/silk_encode_test.go b/silk_encode_test.go new file mode 100644 index 0000000..b273467 --- /dev/null +++ b/silk_encode_test.go @@ -0,0 +1,53 @@ +// SPDX-FileCopyrightText: 2026 The Pion community +// SPDX-License-Identifier: MIT + +package opus + +import ( + "math" + "testing" + + "github.com/stretchr/testify/assert" + "github.com/stretchr/testify/require" +) + +// TestEncodeSILKRoundTrip encodes a SILK wideband frame with the public +// encoder and decodes the resulting Opus packet with the public decoder, +// checking the packet is well-formed and decodes to non-silent audio. +func TestEncodeSILKRoundTrip(t *testing.T) { + enc, err := NewEncoder() + require.NoError(t, err) + + pcm := make([]int16, silkWidebandSampleCount) + for i := range pcm { + pcm[i] = int16(5000*math.Sin(2*math.Pi*float64(i)/48) + 1200*math.Sin(2*math.Pi*float64(i)/11)) + } + + packet := make([]byte, 1275) + n, err := enc.EncodeSILK(pcm, BandwidthWideband, packet) + require.NoError(t, err) + require.Greater(t, n, 1) + packet = packet[:n] + + dec, err := NewDecoderWithOutput(16000, 1) + require.NoError(t, err) + + out := make([]float32, silkWidebandSampleCount) + got, err := dec.DecodeToFloat32(packet, out) + require.NoError(t, err) + require.Positive(t, got) + + var energy float64 + for _, v := range out[:got] { + energy += float64(v) * float64(v) + } + assert.Positive(t, energy, "decoded output is silent") +} + +func TestEncodeSILKInvalidLength(t *testing.T) { + enc, err := NewEncoder() + require.NoError(t, err) + + _, err = enc.EncodeSILK(make([]int16, 100), BandwidthWideband, make([]byte, 1275)) + require.Error(t, err) +}