diff --git a/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts b/packages/stretch-phase-vocoder/src/PhaseVocoder.spec.ts index a1dd9d6..a738313 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,198 @@ 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 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 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); + for (let i = 0; i < out.length; i += 2) { + expect(Math.abs(out[i] - out[i + 1])).toBeLessThanOrEqual(1e-6); + } + }); + + 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 + // 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); + 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..da99dcf 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,24 @@ 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 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. * - * 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 +109,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 +120,53 @@ 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; - - // Overlap-add accumulator (fftSize samples). + // 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/side spectra, per-bin energy, peak indices. + private midRe: Float32Array; + private midIm: Float32Array; + private sideRe: Float32Array; + private sideIm: Float32Array; + private binEnergy: Float32Array; + private peaks: Int32Array; + + // 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 +177,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 +189,23 @@ 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.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.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); + this.olaNorm = new Float32Array(N); this.inputScratch = new Float32Array(Ha * 2); this.outputScratch = new Float32Array(N * 2); } @@ -183,7 +244,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 +258,17 @@ 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.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; } /** @@ -293,129 +363,238 @@ 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 { - 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]); + /** + * Runs the peak-locked phase-vocoder recursion on the mid (L+R) reference + * and resynthesizes both channel spectra in-place. + * + * @remarks + * Identity phase locking (Puckette 1995; Laroche & Dolson 1999): the + * 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 + * 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. + * @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 midRe = this.midRe; + const midIm = this.midIm; + 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]; + 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) { + // 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 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 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 ? binEnergy[k - 1] : -1; + const right = k < bins - 1 ? binEnergy[k + 1] : -1; + if (binEnergy[k] >= left && binEnergy[k] >= right) { + peaks[peakCount++] = k; } + } - // 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]; + // 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; + + // 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(refIm[k], refRe[k]); + // arg(cur · conj(prev)) = phaseCur − phasePrev, already in (−π, π]. + const phaseDiff = Math.atan2( + refIm[k] * prevRe[k] - refRe[k] * prevIm[k], + refRe[k] * prevRe[k] + refIm[k] * prevIm[k], + ); + + const expectedAdvance = k * twoPiOverN * Ha; + let delta = phaseDiff - expectedAdvance; + delta -= twoPi * Math.round(delta / twoPi); // wrap to (−π, π] + + const trueFreq = k * twoPiOverN + delta / Ha; + 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/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]; + 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; } - ifft(re, im); - // Imaginary part is discarded; re holds the real synthesis output. + 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++) { + 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`.