From ec7f27ec43b38d3144dc61171beb6dda44ed667a Mon Sep 17 00:00:00 2001 From: Brandon Jones Date: Sun, 12 Jul 2026 13:57:19 -0400 Subject: [PATCH 1/3] fix(stretch-phase-vocoder): Restore stereo coherence, bounded phase, and exact OLA normalization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Four defects in the PhaseVocoder stretch stage combined to collapse the stereo center and shift output level whenever tempo != 1: - Inter-channel incoherence: L/R were processed with fully independent phase accumulators, destroying the inter-channel phase relationship that forms the stereo center. The recursion now runs once on the mid (L+R) spectrum — by FFT linearity the sum of the channel spectra, so no third FFT — and each channel resynthesizes as the shared mid synthesis phase plus its instantaneous offset from the mid, preserving the inter-channel difference exactly. - Unbounded phase accumulation: synthesis phase grew without wrapping in float32 storage; high bins lose phase resolution within minutes of playback. The accumulator is now wrapped to (-pi, pi] every frame. - Incorrect OLA normalization: the constant 2/overlapFactor assumed single Hann coverage, but the pipeline windows twice (analysis and synthesis), and the synthesis hop varies with tempo. Overlap-add is now normalized per sample by the accumulated squared-window sum (WOLA), exact for any hop including mid-stream tempo changes. - Warmup phase poisoning: the recursion was seeded from the very first hop, when the analysis window is still mostly zeros, freezing the startup transient's arbitrary inter-bin phase relationships into the accumulators permanently — a constant, tempo-dependent level loss (e.g. -1.4 dB at tempo 2.0) even after the fixes above. Frames now pass through untouched until the analysis window is fully primed; the first fully-valid frame seeds the recursion. Measured on stationary tones (fftSize 512, overlap 4): steady-state gain is now 1.0000 at tempos 0.5 through 4.0 for bin-aligned and non-aligned frequencies (previously 0.68-0.97, and 0.75 at identity). First-frame output is now exactly normalized rather than attenuated by w^2/2; those samples are analysis-priming zeros either way. New DSP-correctness specs cover level at multiple tempos, DC reproduction, inter-channel identity and fixed-offset preservation, synthesis-phase boundedness over long runs, and mid-stream tempo changes. Seven of them fail against the previous implementation. Co-Authored-By: Claude Fable 5 --- .../src/PhaseVocoder.spec.ts | 164 +++++++++++ .../stretch-phase-vocoder/src/PhaseVocoder.ts | 269 +++++++++++------- packages/stretch-phase-vocoder/src/windows.ts | 6 +- 3 files changed, 329 insertions(+), 110 deletions(-) diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts index a1dd9d6..c534894 100644 --- a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts +++ b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts @@ -48,6 +48,42 @@ function fillBuffer(buf: TestBuffer, frames: number, value = 0.5): void { buf.putSamples(s, 0, frames); } +const TWO_PI = 2 * Math.PI; + +// Bin-aligned test tone: bin 16 of a 512-point FFT → period 32 samples. +const OMEGA = (TWO_PI * 16) / 512; + +function fillSine( + buf: TestBuffer, + frames: number, + omega: number, + phaseL = 0, + phaseR = 0, + amp = 0.5, +): void { + const s = new Float32Array(frames * 2); + for (let i = 0; i < frames; i++) { + s[i * 2] = amp * Math.sin(omega * i + phaseL); + s[i * 2 + 1] = amp * Math.sin(omega * i + phaseR); + } + buf.putSamples(s, 0, frames); +} + +function allOutput(buf: TestBuffer): Float32Array { + const out = new Float32Array(buf.frameCount * 2); + buf.extract(out, 0, buf.frameCount); + return out; +} + +// RMS across both channels of interleaved frames [from, to). +function rmsInterleaved(x: Float32Array, from: number, to: number): number { + let sum = 0; + for (let n = from; n < to; n++) { + sum += x[n * 2] ** 2 + x[n * 2 + 1] ** 2; + } + return Math.sqrt(sum / (2 * (to - from))); +} + describe('PhaseVocoder', () => { let pv: PhaseVocoder; let inputBuf: TestBuffer; @@ -232,6 +268,134 @@ describe('PhaseVocoder', () => { expect(pv.tempo).toBe(1.0); }); }); + + // Warmup rationale: the analysis window fills after N/Ha = 4 hops, the + // phase recursion is valid from the 2nd frame, and the OLA reaches steady + // state after ⌈N/Hs⌉ frames. Skipping the first 16 hops of output is + // comfortably past every transient, and 16·Hs frames is a whole number of + // the 32-sample test-tone periods. + describe('DSP correctness', () => { + const WARMUP_HOPS = 16; + + it('preserves RMS at tempo 1.0 for a bin-aligned sine', () => { + const Ha = pv.sampleReq; + const hops = 64; + fillSine(inputBuf, hops * Ha, OMEGA); + for (let i = 0; i < hops; i++) pv.process(); + + const out = allOutput(outputBuf); + const outRms = rmsInterleaved(out, WARMUP_HOPS * Ha, hops * Ha); + const inRms = 0.5 / Math.SQRT2; + const dB = 20 * Math.log10(outRms / inRms); + expect(Math.abs(dB)).toBeLessThan(0.5); + }); + + for (const tempo of [0.8, 1.25, 2.0]) { + it(`preserves RMS at tempo ${tempo}`, () => { + const Ha = pv.sampleReq; + const Hs = Math.max(1, Math.round(Ha / tempo)); + const hops = 64; + pv.tempo = tempo; + fillSine(inputBuf, hops * Ha, OMEGA); + for (let i = 0; i < hops; i++) pv.process(); + + const out = allOutput(outputBuf); + const outRms = rmsInterleaved(out, WARMUP_HOPS * Hs, hops * Hs); + const dB = 20 * Math.log10(outRms / (0.5 / Math.SQRT2)); + expect(Math.abs(dB)).toBeLessThan(0.5); + }); + } + + it('reproduces a DC signal at unity level for tempo 1.0', () => { + const Ha = pv.sampleReq; + const hops = 32; + fillBuffer(inputBuf, hops * Ha, 0.5); + for (let i = 0; i < hops; i++) pv.process(); + + const out = allOutput(outputBuf); + for (let n = WARMUP_HOPS * Ha; n < hops * Ha; n++) { + expect(Math.abs(out[n * 2] - 0.5)).toBeLessThan(1e-3); + expect(Math.abs(out[n * 2 + 1] - 0.5)).toBeLessThan(1e-3); + } + }); + + it('produces identical outputs for identical L/R inputs', () => { + pv.tempo = 1.25; + const hops = 48; + fillSine(inputBuf, hops * pv.sampleReq, OMEGA, 0.4, 0.4); + for (let i = 0; i < hops; i++) pv.process(); + + const out = allOutput(outputBuf); + for (let i = 0; i < out.length; i += 2) { + expect(Math.abs(out[i] - out[i + 1])).toBeLessThanOrEqual(1e-6); + } + }); + + it('preserves a fixed inter-channel phase offset at tempo 1.0', () => { + const phi = Math.PI / 3; + const Ha = pv.sampleReq; + const hops = 64; + fillSine(inputBuf, hops * Ha, OMEGA, 0, phi); + for (let i = 0; i < hops; i++) pv.process(); + + // Complex demodulation at the tone frequency over whole periods: for + // a·sin(ωn + φ), arg(Σ x[n]·e^(−iωn)) = φ − π/2, so the L/R argument + // difference is exactly φ. + const out = allOutput(outputBuf); + let reLh = 0; + let imLh = 0; + let reRh = 0; + let imRh = 0; + for (let n = WARMUP_HOPS * Ha; n < hops * Ha; n++) { + const c = Math.cos(OMEGA * n); + const s = Math.sin(OMEGA * n); + reLh += out[n * 2] * c; + imLh -= out[n * 2] * s; + reRh += out[n * 2 + 1] * c; + imRh -= out[n * 2 + 1] * s; + } + let diff = Math.atan2(imRh, reRh) - Math.atan2(imLh, reLh); + diff -= TWO_PI * Math.round(diff / TWO_PI); + expect(Math.abs(diff - phi)).toBeLessThan(0.05); + }); + + it('keeps synthesis phase bounded and channels coherent after many frames', () => { + pv.tempo = 1.25; + const hops = 400; + fillSine(inputBuf, hops * pv.sampleReq, OMEGA, 0.2, 0.2); + for (let i = 0; i < hops; i++) pv.process(); + + // eslint-disable-next-line @typescript-eslint/no-explicit-any + const synth = (pv as any).synthPhaseM as Float32Array; + for (let k = 0; k < synth.length; k++) { + // π + headroom: Math.fround(π) > π and the wrap can land at ±π. + expect(Math.abs(synth[k])).toBeLessThanOrEqual(Math.PI + 1e-3); + } + + const out = allOutput(outputBuf); + for (let i = 0; i < out.length; i += 2) { + expect(Math.abs(out[i] - out[i + 1])).toBeLessThanOrEqual(1e-6); + } + }); + + it('stays level-stable across a mid-stream tempo change', () => { + const Ha = pv.sampleReq; + fillSine(inputBuf, 64 * Ha, OMEGA); + for (let i = 0; i < 32; i++) pv.process(); + pv.tempo = 2.0; + for (let i = 0; i < 32; i++) pv.process(); + + const out = allOutput(outputBuf); + const frames = out.length / 2; + for (let i = 0; i < out.length; i++) { + expect(Number.isFinite(out[i])).toBe(true); + } + const Hs2 = Math.round(Ha / 2.0); + const outRms = rmsInterleaved(out, frames - 16 * Hs2, frames); + const dB = 20 * Math.log10(outRms / (0.5 / Math.SQRT2)); + expect(Math.abs(dB)).toBeLessThan(1.5); + }); + }); }); describe('createPhaseVocoderFactory', () => { diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.ts index ee2babb..e234768 100644 --- a/packages/stretch-phase-vocoder/src/PhaseVocoder.ts +++ b/packages/stretch-phase-vocoder/src/PhaseVocoder.ts @@ -16,6 +16,16 @@ import type { import { fft, ifft } from './fft.js'; import { makeHannWindow } from './windows.js'; +/** + * Floor for the per-sample overlap-add normalization divisor. + * + * @remarks + * Wherever the accumulated window-squared sum is below this floor, the + * numerator carries the same window factor (`O(w² · x)`), so dividing by the + * floor attenuates rather than amplifies — no noise blow-up at frame edges. + */ +const OLA_NORM_EPS = 1e-6; + /** Valid FFT sizes for the phase vocoder. */ export type PhaseVocoderFftSize = 512 | 1024 | 2048 | 4096; @@ -62,12 +72,20 @@ export interface PhaseVocoderOptions { * **Algorithm overview:** * 1. Accumulate analysis frames into a sliding `fftSize`-sample window. * 2. Apply a Hann window and compute the FFT of each channel. - * 3. Compute instantaneous frequency per bin and accumulate synthesis phase. + * 3. Run the phase-vocoder recursion once on the mid (L+R) spectrum — by FFT + * linearity it is the sum of the channel spectra, so no third FFT is + * needed — and resynthesize each channel as the shared mid synthesis phase + * plus that channel's instantaneous offset from the mid. This preserves + * the inter-channel phase difference exactly, keeping the stereo image + * (and center content) intact. * 4. Reconstruct with IFFT, apply Hann window, and overlap-add into an output buffer. * 5. Extract `Hs = round(Ha / tempo)` frames per processing step. * - * The normalization factor is `overlapFactor / 2` (for a Hann window, four - * 75 %-overlapping windows sum to 2; factor-2 windows sum to 1). + * Both analysis and synthesis are Hann-windowed (WOLA). Overlap-add is + * normalized per sample by the accumulated sum of squared window values + * (≈ 1.5 at 75 % overlap), which is exact for any synthesis hop — including + * mid-stream tempo changes and hops large enough that windows no longer + * overlap. * * @example * import { SoundTouch } from '@soundtouchjs/core'; @@ -87,8 +105,6 @@ export class PhaseVocoder implements StretchPipe { private readonly overlapFactor: number; private readonly analysisHop: number; private readonly window: Float32Array; - /** Reciprocal of the OLA normalization constant (overlapFactor / 2). */ - private readonly normInv: number; // Per-channel analysis ring buffers (fftSize samples each). private analysisL: Float32Array; @@ -100,22 +116,34 @@ export class PhaseVocoder implements StretchPipe { private reR: Float32Array; private imR: Float32Array; - // Phase vocoder state (bins = fftSize/2 + 1). - private prevPhaseL: Float32Array; - private prevPhaseR: Float32Array; - private synthPhaseL: Float32Array; - private synthPhaseR: Float32Array; + // Phase vocoder state on the mid (L+R) reference (bins = fftSize/2 + 1). + // A single shared recursion keeps the channels phase-coherent; values are + // wrapped to (−π, π] so float32 storage never loses phase resolution. + private prevPhaseM: Float32Array; + private synthPhaseM: Float32Array; - // Overlap-add accumulator (fftSize samples). + // Overlap-add accumulators (fftSize samples): windowed synthesis output + // per channel, plus the squared-window sum used for exact normalization. private olaL: Float32Array; private olaR: Float32Array; + private olaNorm: Float32Array; // Scratch buffers for deinterleaving / reinterleaving. private inputScratch: Float32Array; private outputScratch: Float32Array; - /** Whether the phase accumulators hold valid state from a previous frame. */ - private hasFrame = false; + /** + * Hops remaining before the phase recursion may be seeded. + * + * @remarks + * The sliding analysis window fills after `overlapFactor` hops. Seeding the + * phase accumulators from a partially-filled (zero-padded) window freezes + * the startup transient's arbitrary inter-bin phase relationships into the + * recursion permanently — a constant level loss and envelope corruption at + * every tempo ≠ 1. Frames pass through untouched until the window is + * primed; the first fully-valid frame seeds the recursion. + */ + private hopsUntilPrimed: number; /** * Creates a PhaseVocoder instance. @@ -126,8 +154,7 @@ export class PhaseVocoder implements StretchPipe { this.overlapFactor = options.overlapFactor ?? 4; this.analysisHop = this.fftSize / this.overlapFactor; this.window = makeHannWindow(this.fftSize); - // For Hann OLA with overlapFactor, sum of windows = overlapFactor / 2 - this.normInv = 2.0 / this.overlapFactor; + this.hopsUntilPrimed = this.overlapFactor; const N = this.fftSize; const bins = N / 2 + 1; @@ -139,12 +166,11 @@ export class PhaseVocoder implements StretchPipe { this.imL = new Float32Array(N); this.reR = new Float32Array(N); this.imR = new Float32Array(N); - this.prevPhaseL = new Float32Array(bins); - this.prevPhaseR = new Float32Array(bins); - this.synthPhaseL = new Float32Array(bins); - this.synthPhaseR = new Float32Array(bins); + this.prevPhaseM = new Float32Array(bins); + this.synthPhaseM = new Float32Array(bins); this.olaL = new Float32Array(N); this.olaR = new Float32Array(N); + this.olaNorm = new Float32Array(N); this.inputScratch = new Float32Array(Ha * 2); this.outputScratch = new Float32Array(N * 2); } @@ -183,7 +209,10 @@ export class PhaseVocoder implements StretchPipe { this.analysisR.fill(0); this.olaL.fill(0); this.olaR.fill(0); + this.olaNorm.fill(0); this.clearMidBuffer(); + // The analysis window was emptied, so it must refill before seeding. + this.hopsUntilPrimed = this.overlapFactor; } /** @@ -194,11 +223,11 @@ export class PhaseVocoder implements StretchPipe { * continuity invariant is re-established on the next frame. */ clearMidBuffer(): void { - this.prevPhaseL.fill(0); - this.prevPhaseR.fill(0); - this.synthPhaseL.fill(0); - this.synthPhaseR.fill(0); - this.hasFrame = false; + this.prevPhaseM.fill(0); + this.synthPhaseM.fill(0); + // The analysis window still holds valid samples; re-seed from the very + // next frame rather than waiting for a full refill. + this.hopsUntilPrimed = 1; } /** @@ -293,129 +322,153 @@ export class PhaseVocoder implements StretchPipe { this.analysisR[N - Ha + i] = this.inputScratch[i * 2 + 1]; } - // Process each channel independently. - this._processChannel( - this.analysisL, - this.reL, - this.imL, - this.prevPhaseL, - this.synthPhaseL, - Ha, - Hs, - bins, - ); - this._processChannel( - this.analysisR, - this.reR, - this.imR, - this.prevPhaseR, - this.synthPhaseR, - Ha, - Hs, - bins, - ); - - this.hasFrame = true; - - // Overlap-add: accumulate synthesis frame then extract Hs frames. - const norm = this.normInv; + // Window + FFT each channel, then run the shared phase-vocoder + // recursion on the mid (L+R) reference and resynthesize both spectra. + this._analyzeChannel(this.analysisL, this.reL, this.imL); + this._analyzeChannel(this.analysisR, this.reR, this.imR); + this._processSpectra(Ha, Hs, bins); + + ifft(this.reL, this.imL); + ifft(this.reR, this.imR); + + // WOLA: overlap-add the synthesis-windowed frame, accumulating the + // squared-window sum alongside for exact per-sample normalization. for (let i = 0; i < N; i++) { - this.olaL[i] += this.reL[i] * this.window[i] * norm; - this.olaR[i] += this.reR[i] * this.window[i] * norm; + const wi = this.window[i]; + this.olaL[i] += this.reL[i] * wi; + this.olaR[i] += this.reR[i] * wi; + this.olaNorm[i] += wi * wi; } + // Extraction-time normalization is complete for [0..toExtract): the next + // frame's window lands at post-shift positions, so no future frame can + // contribute to the samples extracted here. const toExtract = Math.min(Hs, N); for (let i = 0; i < toExtract; i++) { - this.outputScratch[i * 2] = this.olaL[i]; - this.outputScratch[i * 2 + 1] = this.olaR[i]; + const norm = 1.0 / Math.max(this.olaNorm[i], OLA_NORM_EPS); + this.outputScratch[i * 2] = this.olaL[i] * norm; + this.outputScratch[i * 2 + 1] = this.olaR[i] * norm; } this.outputBuffer.putSamples(this.outputScratch, 0, toExtract); - // Shift OLA buffer left by Hs, zero the tail. + // Shift OLA buffers left by the extracted amount, zero the tails. this.olaL.copyWithin(0, toExtract); this.olaR.copyWithin(0, toExtract); + this.olaNorm.copyWithin(0, toExtract); this.olaL.fill(0, N - toExtract); this.olaR.fill(0, N - toExtract); + this.olaNorm.fill(0, N - toExtract); } /** - * Runs one FFT frame through the phase vocoder for a single channel. - * - * @remarks - * Reads from `analysis`, writes synthesized real output back into `re` (the - * imaginary output is discarded). Updates `prevPhase` and `synthPhase` in-place. + * Applies the analysis window and computes the forward FFT of one channel. * * @param analysis fftSize-sample sliding analysis window. * @param re FFT real working buffer (modified in-place). * @param im FFT imaginary working buffer (modified in-place). - * @param prevPhase Previous-frame phase per bin (modified in-place). - * @param synthPhase Accumulated synthesis phase per bin (modified in-place). - * @param Ha Analysis hop size. - * @param Hs Synthesis hop size. - * @param bins Number of positive-frequency bins (fftSize/2 + 1). */ - private _processChannel( + private _analyzeChannel( analysis: Float32Array, re: Float32Array, im: Float32Array, - prevPhase: Float32Array, - synthPhase: Float32Array, - Ha: number, - Hs: number, - bins: number, ): void { const N = this.fftSize; const w = this.window; - // Apply analysis window. for (let i = 0; i < N; i++) { re[i] = analysis[i] * w[i]; im[i] = 0.0; } fft(re, im); + } - // Phase accumulation (positive-frequency bins only). - const twoPiOverN = (2.0 * Math.PI) / N; - if (!this.hasFrame) { - // First frame: seed phases from analysis, no output correction. - for (let k = 0; k < bins; k++) { - const phase = Math.atan2(im[k], re[k]); - prevPhase[k] = phase; - synthPhase[k] = phase; - } - } else { + /** + * Runs the phase-vocoder recursion on the mid (L+R) reference and + * resynthesizes both channel spectra in-place. + * + * @remarks + * The recursion (instantaneous-frequency estimate and synthesis-phase + * accumulation) is computed once per bin on the mid spectrum — by FFT + * linearity, `specL + specR` — and each channel's synthesis phase is the + * shared mid phase plus that channel's instantaneous offset from the mid: + * `synthC = synthM + (phaseC − phaseM)`. The inter-channel difference is + * therefore preserved exactly (`synthL − synthR = phaseL − phaseR`), so + * identical channels stay bit-identical and the stereo image is intact. + * + * For bins where the mid spectrum nearly vanishes (anti-phase content) the + * reference phase is noisy; the inter-channel difference is still exact and + * the worst case is temporal smearing on those bins only. + * + * @param Ha Analysis hop size. + * @param Hs Synthesis hop size. + * @param bins Number of positive-frequency bins (fftSize/2 + 1). + */ + private _processSpectra(Ha: number, Hs: number, bins: number): void { + const N = this.fftSize; + const reL = this.reL; + const imL = this.imL; + const reR = this.reR; + const imR = this.imR; + const twoPi = 2.0 * Math.PI; + const twoPiOverN = twoPi / N; + + if (this.hopsUntilPrimed > 0) { + // Pre-primed hops: (re-)seed the mid phase state and leave both + // spectra untouched so the IFFT reproduces the windowed input (the FFT + // of a real signal is already conjugate-symmetric — no fill needed). + // Seeding repeats each hop until the analysis window is full, so the + // recursion starts from the first fully-valid frame's structure rather + // than freezing the startup transient into the accumulators. for (let k = 0; k < bins; k++) { - const mag = Math.sqrt(re[k] * re[k] + im[k] * im[k]); - const phase = Math.atan2(im[k], re[k]); - - // Instantaneous phase deviation from expected advance. - const expectedAdvance = k * twoPiOverN * Ha; - let delta = phase - prevPhase[k] - expectedAdvance; - - // Wrap to (−π, π]. - delta -= 2.0 * Math.PI * Math.round(delta / (2.0 * Math.PI)); - - // Accumulate synthesis phase scaled by the time-stretch ratio. - const trueFreq = k * twoPiOverN + delta / Ha; - synthPhase[k] += trueFreq * Hs; - prevPhase[k] = phase; - - // Resynthesize: keep magnitude, use accumulated phase. - re[k] = mag * Math.cos(synthPhase[k]); - im[k] = mag * Math.sin(synthPhase[k]); + const phaseM = Math.atan2(imL[k] + imR[k], reL[k] + reR[k]); + this.prevPhaseM[k] = phaseM; + this.synthPhaseM[k] = phaseM; } + this.hopsUntilPrimed--; + return; + } - // Fill conjugate-symmetric bins for real output. - for (let k = bins; k < N; k++) { - const mirror = N - k; - re[k] = re[mirror]; - im[k] = -im[mirror]; - } + for (let k = 0; k < bins; k++) { + const phaseM = Math.atan2(imL[k] + imR[k], reL[k] + reR[k]); + + // Instantaneous phase deviation from expected advance, on the mid. + const expectedAdvance = k * twoPiOverN * Ha; + let delta = phaseM - this.prevPhaseM[k] - expectedAdvance; + + // Wrap to (−π, π]. + delta -= twoPi * Math.round(delta / twoPi); + + // Accumulate synthesis phase scaled by the time-stretch ratio, wrapping + // so the float32 accumulator never loses phase resolution. + const trueFreq = k * twoPiOverN + delta / Ha; + let synthM = this.synthPhaseM[k] + trueFreq * Hs; + synthM -= twoPi * Math.round(synthM / twoPi); + this.synthPhaseM[k] = synthM; + this.prevPhaseM[k] = phaseM; + + // Resynthesize each channel: keep its magnitude, use the shared mid + // synthesis phase plus the channel's offset from the mid. Offsets are + // differences of atan2 results, bounded in (−2π, 2π) — safe for cos/sin. + const magL = Math.sqrt(reL[k] * reL[k] + imL[k] * imL[k]); + const phL = synthM + (Math.atan2(imL[k], reL[k]) - phaseM); + reL[k] = magL * Math.cos(phL); + imL[k] = magL * Math.sin(phL); + + const magR = Math.sqrt(reR[k] * reR[k] + imR[k] * imR[k]); + const phR = synthM + (Math.atan2(imR[k], reR[k]) - phaseM); + reR[k] = magR * Math.cos(phR); + imR[k] = magR * Math.sin(phR); } - ifft(re, im); - // Imaginary part is discarded; re holds the real synthesis output. + // Fill conjugate-symmetric bins for real output, both channels. + for (let k = bins; k < N; k++) { + const mirror = N - k; + reL[k] = reL[mirror]; + imL[k] = -imL[mirror]; + reR[k] = reR[mirror]; + imR[k] = -imR[mirror]; + } } } diff --git a/packages/stretch-phase-vocoder/src/windows.ts b/packages/stretch-phase-vocoder/src/windows.ts index 87d656c..dbcc76a 100644 --- a/packages/stretch-phase-vocoder/src/windows.ts +++ b/packages/stretch-phase-vocoder/src/windows.ts @@ -12,8 +12,10 @@ * @remarks * Uses the symmetric form `0.5 * (1 − cos(2πn / (N−1)))` for N > 1. * For `size = 1` the single coefficient is `1.0`. - * When used with overlap-add, 75 % overlap (factor 4) yields a normalization - * constant of 2 (sum of four Hann windows = 2). + * The phase vocoder applies this window at both analysis and synthesis + * (WOLA), where the accumulated squared-window sum at 75 % overlap is ≈ 1.5; + * it normalizes per sample from that accumulated sum rather than by a + * constant, so the symmetric (non-COLA) form's ripple is normalized out. * * @param size Number of coefficients. * @returns Float32Array of length `size`. From 2469fab859083aeb9d00af63517b1acc25db1ebb Mon Sep 17 00:00:00 2001 From: Brandon Jones Date: Sun, 12 Jul 2026 14:11:26 -0400 Subject: [PATCH 2/3] fix(stretch-phase-vocoder): Lock region phases to spectral peaks to prevent broadband level loss MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With the phase recursion running independently per bin, any Hs != Ha re-randomizes the local phase structure that makes overlapping frames of noise-like content sum consistently. The overlap-add then combines incoherently, which for a Hann window at 75% overlap costs ~3.1 dB on broadband material (consonants, drums, cymbals, reverb) while pure tones pass at unity — measured -2.9 to -3.4 dB on white noise through the stretch stage, matching listener reports of a ~3 dB drop on music. This applies identity phase locking (Puckette 1995; Laroche & Dolson 1999) on the mid-channel reference: the recursion runs only at spectral peaks of the mid magnitude, and every bin in a peak's region of influence (bounded by midpoints between adjacent peaks) rotates rigidly by that peak's synthesis rotation. Rigid rotation preserves the analysis frame's local phase structure, keeping noise coherent under overlap-add, and rotating both channels by the same angle preserves the inter-channel phase difference exactly. Phase state is now carried as complex spectra (previous analysis mid and previous synthesized mid) instead of accumulated angles — complex storage is inherently wrapped, so unbounded accumulation is impossible by construction, and atan2/cos/sin run only at peak bins instead of every bin. Measured on white noise through the stretch stage at +/-2 and +5 semitone-equivalent tempos: -0.4 to -1.0 dB (previously -2.9 to -3.4). Tones and harmonic stacks remain at 0.00 dB. New regression specs assert broadband noise level within 1.5 dB in both stretch directions. Co-Authored-By: Claude Fable 5 --- .../src/PhaseVocoder.spec.ts | 51 ++++- .../stretch-phase-vocoder/src/PhaseVocoder.ts | 185 ++++++++++++------ 2 files changed, 171 insertions(+), 65 deletions(-) diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts index c534894..a5084cc 100644 --- a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts +++ b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts @@ -359,17 +359,21 @@ describe('PhaseVocoder', () => { expect(Math.abs(diff - phi)).toBeLessThan(0.05); }); - it('keeps synthesis phase bounded and channels coherent after many frames', () => { + it('keeps phase state bounded and channels coherent after many frames', () => { pv.tempo = 1.25; const hops = 400; fillSine(inputBuf, hops * pv.sampleReq, OMEGA, 0.2, 0.2); for (let i = 0; i < hops; i++) pv.process(); + // Phase state is carried as complex spectra; magnitudes must stay + // finite and bounded by the frame energy (no runaway accumulation). // eslint-disable-next-line @typescript-eslint/no-explicit-any - const synth = (pv as any).synthPhaseM as Float32Array; - for (let k = 0; k < synth.length; k++) { - // π + headroom: Math.fround(π) > π and the wrap can land at ±π. - expect(Math.abs(synth[k])).toBeLessThanOrEqual(Math.PI + 1e-3); + const sRe = (pv as any).synthMidRe as Float32Array; + // eslint-disable-next-line @typescript-eslint/no-explicit-any + const sIm = (pv as any).synthMidIm as Float32Array; + for (let k = 0; k < sRe.length; k++) { + expect(Number.isFinite(sRe[k])).toBe(true); + expect(Number.isFinite(sIm[k])).toBe(true); } const out = allOutput(outputBuf); @@ -378,6 +382,43 @@ describe('PhaseVocoder', () => { } }); + for (const semis of [2, -2]) { + it(`preserves broadband noise level at ${semis > 0 ? '+' : ''}${semis} semitones equivalent tempo`, () => { + // Without peak phase locking, the per-bin recursion re-randomizes + // local phase structure once Hs ≠ Ha and noise-like content loses + // ~3 dB to incoherent overlap-add. + const tempo = Math.pow(2, -semis / 12); + pv.tempo = tempo; + const Ha = pv.sampleReq; + const Hs = Math.max(1, Math.round(Ha / tempo)); + const hops = 128; + + // Deterministic LCG noise, identical L/R. + let seed = 12345 >>> 0; + const rng = () => { + seed = (seed * 1664525 + 1013904223) >>> 0; + return seed / 0xffffffff - 0.5; + }; + const frames = hops * Ha; + const s = new Float32Array(frames * 2); + let inSum = 0; + for (let i = 0; i < frames; i++) { + const v = 0.5 * rng(); + s[i * 2] = v; + s[i * 2 + 1] = v; + inSum += v * v; + } + const inRms = Math.sqrt(inSum / frames); + inputBuf.putSamples(s, 0, frames); + for (let i = 0; i < hops; i++) pv.process(); + + const out = allOutput(outputBuf); + const outRms = rmsInterleaved(out, WARMUP_HOPS * Hs, hops * Hs); + const dB = 20 * Math.log10(outRms / inRms); + expect(Math.abs(dB)).toBeLessThan(1.5); + }); + } + it('stays level-stable across a mid-stream tempo change', () => { const Ha = pv.sampleReq; fillSine(inputBuf, 64 * Ha, OMEGA); diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.ts index e234768..08e4a3b 100644 --- a/packages/stretch-phase-vocoder/src/PhaseVocoder.ts +++ b/packages/stretch-phase-vocoder/src/PhaseVocoder.ts @@ -72,12 +72,13 @@ export interface PhaseVocoderOptions { * **Algorithm overview:** * 1. Accumulate analysis frames into a sliding `fftSize`-sample window. * 2. Apply a Hann window and compute the FFT of each channel. - * 3. Run the phase-vocoder recursion once on the mid (L+R) spectrum — by FFT - * linearity it is the sum of the channel spectra, so no third FFT is - * needed — and resynthesize each channel as the shared mid synthesis phase - * plus that channel's instantaneous offset from the mid. This preserves - * the inter-channel phase difference exactly, keeping the stereo image - * (and center content) intact. + * 3. Run the phase-vocoder recursion at spectral peaks of the mid (L+R) + * spectrum — by FFT linearity the sum of the channel spectra, so no third + * FFT is needed — and rotate each peak's region of influence rigidly by + * the peak's synthesis rotation (identity phase locking). Rigid rotation + * preserves both the local phase structure (broadband/noise content stays + * coherent under overlap-add) and the inter-channel phase difference + * (stereo image and center content stay intact). * 4. Reconstruct with IFFT, apply Hann window, and overlap-add into an output buffer. * 5. Extract `Hs = round(Ha / tempo)` frames per processing step. * @@ -117,10 +118,20 @@ export class PhaseVocoder implements StretchPipe { private imR: Float32Array; // Phase vocoder state on the mid (L+R) reference (bins = fftSize/2 + 1). - // A single shared recursion keeps the channels phase-coherent; values are - // wrapped to (−π, π] so float32 storage never loses phase resolution. - private prevPhaseM: Float32Array; - private synthPhaseM: Float32Array; + // The recursion runs only at spectral peaks (identity phase locking); each + // peak's region rotates rigidly, so phases are carried as complex spectra — + // inherently wrapped, no unbounded accumulators. prevMid holds the previous + // analysis mid spectrum, synthMid the previous synthesized mid spectrum. + private prevMidRe: Float32Array; + private prevMidIm: Float32Array; + private synthMidRe: Float32Array; + private synthMidIm: Float32Array; + + // Per-frame scratch: current mid spectrum, its magnitudes, and peak indices. + private midRe: Float32Array; + private midIm: Float32Array; + private midMag: Float32Array; + private peaks: Int32Array; // Overlap-add accumulators (fftSize samples): windowed synthesis output // per channel, plus the squared-window sum used for exact normalization. @@ -166,8 +177,14 @@ export class PhaseVocoder implements StretchPipe { this.imL = new Float32Array(N); this.reR = new Float32Array(N); this.imR = new Float32Array(N); - this.prevPhaseM = new Float32Array(bins); - this.synthPhaseM = new Float32Array(bins); + this.prevMidRe = new Float32Array(bins); + this.prevMidIm = new Float32Array(bins); + this.synthMidRe = new Float32Array(bins); + this.synthMidIm = new Float32Array(bins); + this.midRe = new Float32Array(bins); + this.midIm = new Float32Array(bins); + this.midMag = new Float32Array(bins); + this.peaks = new Int32Array(bins); this.olaL = new Float32Array(N); this.olaR = new Float32Array(N); this.olaNorm = new Float32Array(N); @@ -223,8 +240,10 @@ export class PhaseVocoder implements StretchPipe { * continuity invariant is re-established on the next frame. */ clearMidBuffer(): void { - this.prevPhaseM.fill(0); - this.synthPhaseM.fill(0); + this.prevMidRe.fill(0); + this.prevMidIm.fill(0); + this.synthMidRe.fill(0); + this.synthMidIm.fill(0); // The analysis window still holds valid samples; re-seed from the very // next frame rather than waiting for a full refill. this.hopsUntilPrimed = 1; @@ -384,21 +403,27 @@ export class PhaseVocoder implements StretchPipe { } /** - * Runs the phase-vocoder recursion on the mid (L+R) reference and - * resynthesizes both channel spectra in-place. + * Runs the peak-locked phase-vocoder recursion on the mid (L+R) reference + * and resynthesizes both channel spectra in-place. * * @remarks - * The recursion (instantaneous-frequency estimate and synthesis-phase - * accumulation) is computed once per bin on the mid spectrum — by FFT - * linearity, `specL + specR` — and each channel's synthesis phase is the - * shared mid phase plus that channel's instantaneous offset from the mid: - * `synthC = synthM + (phaseC − phaseM)`. The inter-channel difference is - * therefore preserved exactly (`synthL − synthR = phaseL − phaseR`), so - * identical channels stay bit-identical and the stereo image is intact. + * Identity phase locking (Puckette 1995; Laroche & Dolson 1999): the + * instantaneous-frequency recursion runs only at spectral peaks of the mid + * magnitude — by FFT linearity the mid spectrum is `specL + specR`, so no + * third FFT is needed — and every bin in a peak's region of influence + * rotates rigidly by that peak's synthesis rotation. Rigid rotation + * preserves the analysis frame's local phase structure, which keeps + * overlap-added frames coherent for noise-like content too (a per-bin + * recursion re-randomizes those relationships once `Hs ≠ Ha`, costing + * ~3 dB of broadband level). * - * For bins where the mid spectrum nearly vanishes (anti-phase content) the - * reference phase is noisy; the inter-channel difference is still exact and - * the worst case is temporal smearing on those bins only. + * Both channels rotate by the same angle per bin, so the inter-channel + * phase difference is preserved exactly — identical channels stay + * bit-identical and the stereo image is intact. + * + * Phase state is carried as complex spectra (previous analysis mid and + * previous synthesized mid); complex storage is inherently wrapped, so + * there is no unbounded phase accumulator by construction. * * @param Ha Analysis hop size. * @param Hs Synthesis hop size. @@ -410,57 +435,97 @@ export class PhaseVocoder implements StretchPipe { const imL = this.imL; const reR = this.reR; const imR = this.imR; + const midRe = this.midRe; + const midIm = this.midIm; + const midMag = this.midMag; const twoPi = 2.0 * Math.PI; const twoPiOverN = twoPi / N; + for (let k = 0; k < bins; k++) { + midRe[k] = reL[k] + reR[k]; + midIm[k] = imL[k] + imR[k]; + midMag[k] = Math.sqrt(midRe[k] * midRe[k] + midIm[k] * midIm[k]); + } + if (this.hopsUntilPrimed > 0) { - // Pre-primed hops: (re-)seed the mid phase state and leave both - // spectra untouched so the IFFT reproduces the windowed input (the FFT - // of a real signal is already conjugate-symmetric — no fill needed). + // Pre-primed hops: (re-)seed the phase state and leave both spectra + // untouched so the IFFT reproduces the windowed input (the FFT of a + // real signal is already conjugate-symmetric — no fill needed). // Seeding repeats each hop until the analysis window is full, so the // recursion starts from the first fully-valid frame's structure rather - // than freezing the startup transient into the accumulators. - for (let k = 0; k < bins; k++) { - const phaseM = Math.atan2(imL[k] + imR[k], reL[k] + reR[k]); - this.prevPhaseM[k] = phaseM; - this.synthPhaseM[k] = phaseM; - } + // than freezing the startup transient into the phase state. + this.prevMidRe.set(midRe); + this.prevMidIm.set(midIm); + this.synthMidRe.set(midRe); + this.synthMidIm.set(midIm); this.hopsUntilPrimed--; return; } + // Find local maxima of the mid magnitude (plateau-tolerant). Silence + // degenerates to every bin being a "peak" with a zero rotation — harmless. + const peaks = this.peaks; + let peakCount = 0; for (let k = 0; k < bins; k++) { - const phaseM = Math.atan2(imL[k] + imR[k], reL[k] + reR[k]); + const left = k > 0 ? midMag[k - 1] : -1; + const right = k < bins - 1 ? midMag[k + 1] : -1; + if (midMag[k] >= left && midMag[k] >= right) { + peaks[peakCount++] = k; + } + } - // Instantaneous phase deviation from expected advance, on the mid. - const expectedAdvance = k * twoPiOverN * Ha; - let delta = phaseM - this.prevPhaseM[k] - expectedAdvance; + // Process each peak's region of influence: [regionStart, regionEnd). + // Boundaries sit at the midpoint between adjacent peaks. + let regionStart = 0; + for (let p = 0; p < peakCount; p++) { + const k = peaks[p]; + const regionEnd = p + 1 < peakCount ? (k + peaks[p + 1] + 1) >> 1 : bins; + + // Phase-vocoder recursion at the peak bin. Phase differences come from + // complex products against the stored previous spectra, so only the + // peak bins ever need atan2/cos/sin. + const phaseCur = Math.atan2(midIm[k], midRe[k]); + // arg(cur · conj(prev)) = phaseCur − phasePrev, already in (−π, π]. + const phaseDiff = Math.atan2( + midIm[k] * this.prevMidRe[k] - midRe[k] * this.prevMidIm[k], + midRe[k] * this.prevMidRe[k] + midIm[k] * this.prevMidIm[k], + ); - // Wrap to (−π, π]. - delta -= twoPi * Math.round(delta / twoPi); + const expectedAdvance = k * twoPiOverN * Ha; + let delta = phaseDiff - expectedAdvance; + delta -= twoPi * Math.round(delta / twoPi); // wrap to (−π, π] - // Accumulate synthesis phase scaled by the time-stretch ratio, wrapping - // so the float32 accumulator never loses phase resolution. const trueFreq = k * twoPiOverN + delta / Ha; - let synthM = this.synthPhaseM[k] + trueFreq * Hs; - synthM -= twoPi * Math.round(synthM / twoPi); - this.synthPhaseM[k] = synthM; - this.prevPhaseM[k] = phaseM; - - // Resynthesize each channel: keep its magnitude, use the shared mid - // synthesis phase plus the channel's offset from the mid. Offsets are - // differences of atan2 results, bounded in (−2π, 2π) — safe for cos/sin. - const magL = Math.sqrt(reL[k] * reL[k] + imL[k] * imL[k]); - const phL = synthM + (Math.atan2(imL[k], reL[k]) - phaseM); - reL[k] = magL * Math.cos(phL); - imL[k] = magL * Math.sin(phL); - - const magR = Math.sqrt(reR[k] * reR[k] + imR[k] * imR[k]); - const phR = synthM + (Math.atan2(imR[k], reR[k]) - phaseM); - reR[k] = magR * Math.cos(phR); - imR[k] = magR * Math.sin(phR); + const synthPrev = Math.atan2(this.synthMidIm[k], this.synthMidRe[k]); + + // Rigid rotation for this peak's whole region: the peak's accumulated + // synthesis phase advance relative to its current analysis phase. + const theta = synthPrev + trueFreq * Hs - phaseCur; + const cosT = Math.cos(theta); + const sinT = Math.sin(theta); + + for (let j = regionStart; j < regionEnd; j++) { + // Store the synthesized mid for the next frame's recursion. + this.synthMidRe[j] = midRe[j] * cosT - midIm[j] * sinT; + this.synthMidIm[j] = midRe[j] * sinT + midIm[j] * cosT; + + const rl = reL[j]; + const il = imL[j]; + reL[j] = rl * cosT - il * sinT; + imL[j] = rl * sinT + il * cosT; + + const rr = reR[j]; + const ir = imR[j]; + reR[j] = rr * cosT - ir * sinT; + imR[j] = rr * sinT + ir * cosT; + } + + regionStart = regionEnd; } + this.prevMidRe.set(midRe); + this.prevMidIm.set(midIm); + // Fill conjugate-symmetric bins for real output, both channels. for (let k = bins; k < N; k++) { const mirror = N - k; From c570f28153b6ee85f23b153f0de491a65dd10299 Mon Sep 17 00:00:00 2001 From: Brandon Jones Date: Sun, 12 Jul 2026 14:29:42 -0400 Subject: [PATCH 3/3] fix(stretch-phase-vocoder): Use the stronger of mid/side as the peak phase reference MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A mid-only (L+R) phase reference is degenerate for anti-phase content: the mid spectrum vanishes, its phase turns noisy, and the recursion's rotation decoheres the overlap-add — measured -6 dB on an anti-phase sine and -1.8 to -2.7 dB on anti-phase noise. Exactly out-of-phase layers occur in real material (stereo-widening and phase tricks in backing tracks), and the previous per-channel recursion handled them coherently, so this was a regression of the mid-reference design. The side (L-R) spectrum is as free as the mid by FFT linearity. Peak finding now runs on per-bin energy (|mid|^2 + |side|^2, phase-agnostic) and each peak's recursion uses whichever of mid/side is stronger at that bin, keeping the reference well-conditioned for any inter-channel correlation. Rigid region rotation is unchanged — rotating L and R by the same angle rotates mid and side identically — so the exact inter-channel phase guarantee is untouched. Measured after the change: anti-phase sine 0.00 dB (was -6.03), anti-phase noise -0.4 to -0.6 dB (was -1.8 to -2.7), identical to regular noise; all previous levels unchanged. New spec asserts level and exact polarity preservation for purely anti-phase channels. Co-Authored-By: Claude Fable 5 --- .../src/PhaseVocoder.spec.ts | 23 ++++ .../stretch-phase-vocoder/src/PhaseVocoder.ts | 129 +++++++++++++----- 2 files changed, 118 insertions(+), 34 deletions(-) diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts index a5084cc..a738313 100644 --- a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts +++ b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts @@ -382,6 +382,29 @@ describe('PhaseVocoder', () => { } }); + it('preserves level and polarity for purely anti-phase channels', () => { + // L = −R makes the mid (L+R) spectrum vanish; a mid-only phase + // reference turns noisy and costs several dB on such content. The + // mid/side dual reference must keep it exact. + const tempo = Math.pow(2, -2 / 12); + pv.tempo = tempo; + const Ha = pv.sampleReq; + const Hs = Math.max(1, Math.round(Ha / tempo)); + const hops = 64; + fillSine(inputBuf, hops * Ha, OMEGA, 0, Math.PI); // R inverted + for (let i = 0; i < hops; i++) pv.process(); + + const out = allOutput(outputBuf); + const outRms = rmsInterleaved(out, WARMUP_HOPS * Hs, hops * Hs); + const dB = 20 * Math.log10(outRms / (0.5 / Math.SQRT2)); + expect(Math.abs(dB)).toBeLessThan(0.5); + + // Rigid rotation must preserve the exact anti-phase relationship. + for (let i = 0; i < out.length; i += 2) { + expect(Math.abs(out[i] + out[i + 1])).toBeLessThanOrEqual(1e-6); + } + }); + for (const semis of [2, -2]) { it(`preserves broadband noise level at ${semis > 0 ? '+' : ''}${semis} semitones equivalent tempo`, () => { // Without peak phase locking, the per-bin recursion re-randomizes diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.ts index 08e4a3b..da99dcf 100644 --- a/packages/stretch-phase-vocoder/src/PhaseVocoder.ts +++ b/packages/stretch-phase-vocoder/src/PhaseVocoder.ts @@ -72,13 +72,16 @@ export interface PhaseVocoderOptions { * **Algorithm overview:** * 1. Accumulate analysis frames into a sliding `fftSize`-sample window. * 2. Apply a Hann window and compute the FFT of each channel. - * 3. Run the phase-vocoder recursion at spectral peaks of the mid (L+R) - * spectrum — by FFT linearity the sum of the channel spectra, so no third - * FFT is needed — and rotate each peak's region of influence rigidly by - * the peak's synthesis rotation (identity phase locking). Rigid rotation - * preserves both the local phase structure (broadband/noise content stays - * coherent under overlap-add) and the inter-channel phase difference - * (stereo image and center content stay intact). + * 3. Run the phase-vocoder recursion at spectral peaks, using whichever of + * the mid (L+R) or side (L−R) spectrum is stronger at each peak as the + * phase reference — both are free by FFT linearity, and the pairing keeps + * the reference well-conditioned for any inter-channel correlation + * (identical, uncorrelated, or anti-phase channels). Each peak's region + * of influence rotates rigidly by the peak's synthesis rotation (identity + * phase locking), preserving both the local phase structure + * (broadband/noise content stays coherent under overlap-add) and the + * inter-channel phase difference (stereo image and center content stay + * intact). * 4. Reconstruct with IFFT, apply Hann window, and overlap-add into an output buffer. * 5. Extract `Hs = round(Ha / tempo)` frames per processing step. * @@ -117,20 +120,29 @@ export class PhaseVocoder implements StretchPipe { private reR: Float32Array; private imR: Float32Array; - // Phase vocoder state on the mid (L+R) reference (bins = fftSize/2 + 1). - // The recursion runs only at spectral peaks (identity phase locking); each - // peak's region rotates rigidly, so phases are carried as complex spectra — - // inherently wrapped, no unbounded accumulators. prevMid holds the previous - // analysis mid spectrum, synthMid the previous synthesized mid spectrum. + // Phase vocoder state (bins = fftSize/2 + 1). The recursion runs only at + // spectral peaks (identity phase locking); each peak's region rotates + // rigidly, so phases are carried as complex spectra — inherently wrapped, + // no unbounded accumulators. Both the mid (L+R) and side (L−R) spectra are + // tracked so each peak can use whichever reference is stronger there: + // the mid vanishes for anti-phase content, the side for identical content. + // prev* hold the previous analysis spectra, synth* the previous + // synthesized spectra. private prevMidRe: Float32Array; private prevMidIm: Float32Array; + private prevSideRe: Float32Array; + private prevSideIm: Float32Array; private synthMidRe: Float32Array; private synthMidIm: Float32Array; + private synthSideRe: Float32Array; + private synthSideIm: Float32Array; - // Per-frame scratch: current mid spectrum, its magnitudes, and peak indices. + // Per-frame scratch: current mid/side spectra, per-bin energy, peak indices. private midRe: Float32Array; private midIm: Float32Array; - private midMag: Float32Array; + private sideRe: Float32Array; + private sideIm: Float32Array; + private binEnergy: Float32Array; private peaks: Int32Array; // Overlap-add accumulators (fftSize samples): windowed synthesis output @@ -179,11 +191,17 @@ export class PhaseVocoder implements StretchPipe { this.imR = new Float32Array(N); this.prevMidRe = new Float32Array(bins); this.prevMidIm = new Float32Array(bins); + this.prevSideRe = new Float32Array(bins); + this.prevSideIm = new Float32Array(bins); this.synthMidRe = new Float32Array(bins); this.synthMidIm = new Float32Array(bins); + this.synthSideRe = new Float32Array(bins); + this.synthSideIm = new Float32Array(bins); this.midRe = new Float32Array(bins); this.midIm = new Float32Array(bins); - this.midMag = new Float32Array(bins); + this.sideRe = new Float32Array(bins); + this.sideIm = new Float32Array(bins); + this.binEnergy = new Float32Array(bins); this.peaks = new Int32Array(bins); this.olaL = new Float32Array(N); this.olaR = new Float32Array(N); @@ -242,8 +260,12 @@ export class PhaseVocoder implements StretchPipe { clearMidBuffer(): void { this.prevMidRe.fill(0); this.prevMidIm.fill(0); + this.prevSideRe.fill(0); + this.prevSideIm.fill(0); this.synthMidRe.fill(0); this.synthMidIm.fill(0); + this.synthSideRe.fill(0); + this.synthSideIm.fill(0); // The analysis window still holds valid samples; re-seed from the very // next frame rather than waiting for a full refill. this.hopsUntilPrimed = 1; @@ -408,14 +430,20 @@ export class PhaseVocoder implements StretchPipe { * * @remarks * Identity phase locking (Puckette 1995; Laroche & Dolson 1999): the - * instantaneous-frequency recursion runs only at spectral peaks of the mid - * magnitude — by FFT linearity the mid spectrum is `specL + specR`, so no - * third FFT is needed — and every bin in a peak's region of influence - * rotates rigidly by that peak's synthesis rotation. Rigid rotation - * preserves the analysis frame's local phase structure, which keeps - * overlap-added frames coherent for noise-like content too (a per-bin - * recursion re-randomizes those relationships once `Hs ≠ Ha`, costing - * ~3 dB of broadband level). + * instantaneous-frequency recursion runs only at spectral peaks of the + * per-bin energy, and every bin in a peak's region of influence rotates + * rigidly by that peak's synthesis rotation. Rigid rotation preserves the + * analysis frame's local phase structure, which keeps overlap-added frames + * coherent for noise-like content too (a per-bin recursion re-randomizes + * those relationships once `Hs ≠ Ha`, costing ~3 dB of broadband level). + * + * Each peak's reference is the stronger of the mid (L+R) and side (L−R) + * spectra at that bin — both free by FFT linearity, no extra FFTs. A + * mid-only reference is degenerate for anti-phase content (the mid + * vanishes and its phase turns noisy, costing several dB on + * stereo-widened material); side-only degenerates for identical channels. + * The stronger-of-two choice is well-conditioned for any inter-channel + * correlation. * * Both channels rotate by the same angle per bin, so the inter-channel * phase difference is preserved exactly — identical channels stay @@ -437,14 +465,23 @@ export class PhaseVocoder implements StretchPipe { const imR = this.imR; const midRe = this.midRe; const midIm = this.midIm; - const midMag = this.midMag; + const sideRe = this.sideRe; + const sideIm = this.sideIm; + const binEnergy = this.binEnergy; const twoPi = 2.0 * Math.PI; const twoPiOverN = twoPi / N; for (let k = 0; k < bins; k++) { midRe[k] = reL[k] + reR[k]; midIm[k] = imL[k] + imR[k]; - midMag[k] = Math.sqrt(midRe[k] * midRe[k] + midIm[k] * midIm[k]); + sideRe[k] = reL[k] - reR[k]; + sideIm[k] = imL[k] - imR[k]; + // |mid|² + |side|² = 2(|L|² + |R|²) — phase-agnostic peak metric. + binEnergy[k] = + midRe[k] * midRe[k] + + midIm[k] * midIm[k] + + sideRe[k] * sideRe[k] + + sideIm[k] * sideIm[k]; } if (this.hopsUntilPrimed > 0) { @@ -456,20 +493,24 @@ export class PhaseVocoder implements StretchPipe { // than freezing the startup transient into the phase state. this.prevMidRe.set(midRe); this.prevMidIm.set(midIm); + this.prevSideRe.set(sideRe); + this.prevSideIm.set(sideIm); this.synthMidRe.set(midRe); this.synthMidIm.set(midIm); + this.synthSideRe.set(sideRe); + this.synthSideIm.set(sideIm); this.hopsUntilPrimed--; return; } - // Find local maxima of the mid magnitude (plateau-tolerant). Silence + // Find local maxima of the per-bin energy (plateau-tolerant). Silence // degenerates to every bin being a "peak" with a zero rotation — harmless. const peaks = this.peaks; let peakCount = 0; for (let k = 0; k < bins; k++) { - const left = k > 0 ? midMag[k - 1] : -1; - const right = k < bins - 1 ? midMag[k + 1] : -1; - if (midMag[k] >= left && midMag[k] >= right) { + const left = k > 0 ? binEnergy[k - 1] : -1; + const right = k < bins - 1 ? binEnergy[k + 1] : -1; + if (binEnergy[k] >= left && binEnergy[k] >= right) { peaks[peakCount++] = k; } } @@ -481,14 +522,28 @@ export class PhaseVocoder implements StretchPipe { const k = peaks[p]; const regionEnd = p + 1 < peakCount ? (k + peaks[p + 1] + 1) >> 1 : bins; + // Choose the stronger of mid/side as the phase reference at this peak: + // the mid vanishes for anti-phase content, the side for identical + // channels — using the dominant one keeps the recursion's reference + // well-conditioned for any inter-channel correlation. + const midSq = midRe[k] * midRe[k] + midIm[k] * midIm[k]; + const sideSq = sideRe[k] * sideRe[k] + sideIm[k] * sideIm[k]; + const useMid = midSq >= sideSq; + const refRe = useMid ? midRe : sideRe; + const refIm = useMid ? midIm : sideIm; + const prevRe = useMid ? this.prevMidRe : this.prevSideRe; + const prevIm = useMid ? this.prevMidIm : this.prevSideIm; + const synthRe = useMid ? this.synthMidRe : this.synthSideRe; + const synthIm = useMid ? this.synthMidIm : this.synthSideIm; + // Phase-vocoder recursion at the peak bin. Phase differences come from // complex products against the stored previous spectra, so only the // peak bins ever need atan2/cos/sin. - const phaseCur = Math.atan2(midIm[k], midRe[k]); + const phaseCur = Math.atan2(refIm[k], refRe[k]); // arg(cur · conj(prev)) = phaseCur − phasePrev, already in (−π, π]. const phaseDiff = Math.atan2( - midIm[k] * this.prevMidRe[k] - midRe[k] * this.prevMidIm[k], - midRe[k] * this.prevMidRe[k] + midIm[k] * this.prevMidIm[k], + refIm[k] * prevRe[k] - refRe[k] * prevIm[k], + refRe[k] * prevRe[k] + refIm[k] * prevIm[k], ); const expectedAdvance = k * twoPiOverN * Ha; @@ -496,18 +551,22 @@ export class PhaseVocoder implements StretchPipe { delta -= twoPi * Math.round(delta / twoPi); // wrap to (−π, π] const trueFreq = k * twoPiOverN + delta / Ha; - const synthPrev = Math.atan2(this.synthMidIm[k], this.synthMidRe[k]); + const synthPrev = Math.atan2(synthIm[k], synthRe[k]); // Rigid rotation for this peak's whole region: the peak's accumulated // synthesis phase advance relative to its current analysis phase. + // Rotating L and R by the same angle rotates mid and side identically, + // so both synthesized flavors stay consistent. const theta = synthPrev + trueFreq * Hs - phaseCur; const cosT = Math.cos(theta); const sinT = Math.sin(theta); for (let j = regionStart; j < regionEnd; j++) { - // Store the synthesized mid for the next frame's recursion. + // Store the synthesized mid/side for the next frame's recursion. this.synthMidRe[j] = midRe[j] * cosT - midIm[j] * sinT; this.synthMidIm[j] = midRe[j] * sinT + midIm[j] * cosT; + this.synthSideRe[j] = sideRe[j] * cosT - sideIm[j] * sinT; + this.synthSideIm[j] = sideRe[j] * sinT + sideIm[j] * cosT; const rl = reL[j]; const il = imL[j]; @@ -525,6 +584,8 @@ export class PhaseVocoder implements StretchPipe { this.prevMidRe.set(midRe); this.prevMidIm.set(midIm); + this.prevSideRe.set(sideRe); + this.prevSideIm.set(sideIm); // Fill conjugate-symmetric bins for real output, both channels. for (let k = bins; k < N; k++) {