Unmix

An MP3 never stores the waveform. The encoder cuts each channel into granules of 576 samples (two per frame at 32–48 kHz, one per frame at 16–24 kHz, as in the test songs here). It turns each granule into 576 frequency coefficients, rounds those coefficients as coarsely as the ear allows, and stores the rounded integers. The decoder only undoes the transform, so everything it plays is a linear function of those integers.

1. The transform

A bank of 32 filters first splits the signal into 32 equal-width subbands (250 Hz each at 16 kHz). Each subband signal is then split again, into 18 frequencies, by a modified discrete cosine transform (MDCT) over 36 samples that overlap the previous granule by half:

That gives 32 × 18 = 576 coefficients per granule and channel, the same number as samples: the MDCT is critically sampled. Its outputs alias, but with the right window w the aliasing of neighbouring granules cancels on overlap-add (time-domain aliasing cancellation), so without rounding the decoder would rebuild the input exactly. Around sharp attacks the encoder switches to three short MDCTs with N = 6, so that rounding noise doesn't smear audibly before the attack (pre-echo).

2. The rounding

The coefficients are grouped into about 21 scalefactor bands, narrow at low frequencies and wide at high ones, like the ear's critical bands. A psychoacoustic model estimates how much noise each band can hide under the louder sounds around it, and picks a quantiser step Δ_b for the band to stay under that level. The file stores one integer q_k per coefficient, and the decoder turns it back into

where G is the granule's global gain, s_b the band's scalefactor, π_b a fixed pre-emphasis boost for high bands (when the granule turns it on), and m ∈ {½, 1} a per-granule scale. Short blocks add a gain 2^−2g_w per window. The 4/3 power makes the step grow with the amplitude, so loud coefficients are rounded more coarsely in absolute terms, which the ear also tolerates.

3. Packaging

The integers are Huffman-coded, which is lossless. Stereo is usually coded as mid and side, M = (L + R)/√2 and S = (L − R)/√2, switched per granule, because the side channel of most music is small and cheap. Frequencies above the encoder's lowpass are simply zero. All of the loss happens in step 2.

4. What this means for separation

Decoding is linear, x̂ = S X̂, with S the synthesis (IMDCT, window, overlap-add, polyphase filterbank). So the file doesn't pin down one waveform. For every coefficient it says the original coefficient lay in the interval of values that round to q_k. LAME, the encoder used here, rounds to the nearest decoded level, so that interval runs between the midpoints to the neighbouring levels:

mirrored for q_k < 0, and |(Ax)_k| ≤ Δ_b/2 for q_k = 0. Here A is the encoder's analysis (filterbank then MDCT). Round 1 uses only the decoded waveform, which is the centre of every cell. Round 2 (section 07) uses the whole cell as the measurement. A coded zero above the lowpass, or in a band the encoder dropped entirely, carries no information, so those cells are left out.

The measurement

The recording has two channels g_L, g_R of T samples (16 kHz). It is modelled as a voice v, panned to the centre so the same signal appears in both channels, plus a harmonic accompaniment h_c (chords, bass) and percussion p_c in each channel:

That is 5T unknowns and 2T equations. Moving any signal δ from one part to another, (h_c + δ, p_c − δ), leaves g unchanged, so the null space has dimension 3T and the data alone decide nothing about the split. Everything else comes from the priors below.

The voice lives on a comb

F is a short-time Fourier transform (1024-sample window, hop R = 256) scaled to be a tight frame, F^TF = I. Given the tracked pitch f₀(n) in frame n (section 03), the comb M marks the bins near its harmonics, and the voice is built only from coefficients z on the comb (synthesis form, section 06):

counting harmonics j·f₀(n) up to 5 kHz, and M_n,· = 0 in frames whose pitch is outside 75–1100 Hz, the range of human singing. Δf = 15.6 Hz is the bin spacing.

Round 1: the convex problem

• Voice, α_V‖z‖₁: few active coefficients on the comb.
• Harmony, ‖D_tEFh_c‖₁: sustained notes barely change from frame to frame. D_t is the forward difference over frames, (D_tY)_n,k = Y_n+1,k − Y_n,k. A steady tone's coefficient still rotates in phase by its frequency times the hop, so E = diag(e^−iθ_n,k) first removes that rotation, using the mixture's instantaneous frequency (section 04):

• Percussion, Σ_n‖(Fp_c)_n,:‖₂: a group norm over each frame's whole spectrum, so only a few frames (the hits) are active, but broadband when they are.
• The constraint: the three parts must add up to the recording exactly, in both channels. Round 1 also tried a tolerance per frequency band, ‖B_bF(v + h_c + p_c − g_c)‖₂ ≤ ε_b; it made no difference (section 05).

Solving it: PDHG

Every term is convex but non-smooth. Write x = (z, h, p), and give each term i its linear operator K_i: K_Hx = D_tEFh_c, K_Px = Fp_c, K_Dx = F(F^T(M ⊙ z) + h_c + p_c), with f_D the indicator of {Fg_c}. The problem becomes a saddle point,

and the primal–dual hybrid gradient method (Chambolle–Pock) alternates a step on each dual variable, a step on the primal, and an extrapolation:

Each prox is one line. The dual of an ℓ1 term is a projection onto a box, y ↦ y / max(1, |y|/α). The group norm projects each frame onto a ball, y_n,: ↦ y_n,: / max(1, ‖y_n,:‖₂/α_P). By Moreau's identity the data constraint gives y ↦ y − σ_DFg_c. The primal step soft-thresholds the voice and keeps it on the comb,

while h and p take plain gradient steps. The step sizes are σ_i = γ_i/L_i² and τ = 0.99/Σγ_i, with L_i ≥ ‖K_i‖ (‖D_tEF‖ ≤ 2, ‖F‖ = 1). Then τ Σ σ_iL_i² = 0.99 < 1, the condition under which PDHG converges.

The refit

The ℓ1 penalty shrinks the voice, so the parts are used only as soft masks to re-split the recording. With S_j ∈ {Fv, Fh_c, Fp_c},

The masks sum to one in every bin, so the new parts still add up to the recording exactly.

Round 2: the MP3 cells as the data term

Write y_c = v + h_c + p_c for the model's mixture. The equality y_c = g_c is replaced by the cells from part A. The encoder's analysis is A = D⁻¹S^T, with D = diag(S^TS), because S is orthogonal up to a scale per coefficient. Rescaled by Ψ = D^−1/2S^T (so ΨΨ^T = I), the constraint over both coded channels is a box:

The second formula is its dual update, again by Moreau's identity. Only this one line of the solver changes.

Round 2: smoothness that follows the pitch

A constant-Q transform C (36 bins per octave) shifts every harmonic by the same number of bins when the pitch moves. The new voice penalty α_T‖D_θCv‖₁ takes its differences along that shift:

with κ an integer, 0 ≤ β < 1, and |d_t| ≤ 3 (d_t = 0 in unvoiced frames). The unit-modulus factor ρ_k,t removes the expected phase advance, so a steady sung partial has zero difference. In the outer loop it comes from the first solve's voice estimate v̂: with a_k,t = (Cv̂)_k+dt,t+1 · conj((Cv̂)_k,t), ρ_k,t = conj(a_k,t)/|a_k,t|. The full round-2 problem is

The range gate, the third idea, acts before the solve: a pitch frame is kept only if it lies within the line's own register, P₅/1.25 ≤ f₀(t) ≤ 1.25·P₉₅, where P₅ and P₉₅ are percentiles of the tracked pitch.

Why the cells didn't help

The box is only as wide as the codec error. At 128 kbps that is −26.7 dB, about 5% of the signal's amplitude, while separation errors are 4–10 dB down, 30–60% of the voice's amplitude. Relaxing the equality to the box moves the feasible set by much less than the null space the priors have to resolve, so the extra freedom mostly lets the ℓ1 terms shrink the parts a little more. See Results for the numbers.

Each part of the song is placed in the stereo field with a left gain and a right gain, then everything is added. That is a linear map from the parts f to the recording g, the same shape as a tomography scanner turning a body into projections. Here there are only two “views”, left and right, for K parts.

One difference matters. In tomography the trouble is conditioning: some singular values are tiny. Here A is perfectly conditioned on its row space, but its null space has dimension (K − 2)T. So the measurements fix almost nothing, and the priors carry nearly all the information.

All the priors live in the time-frequency plane, so the solver keeps applying a short-time Fourier transform F and its adjoint F^T. PDHG only converges if F^T really is the adjoint and the operator norms are known. Building F as a tight (Parseval) frame makes both easy: ‖Fx‖ = ‖x‖, so its norm is exactly 1.

The strongest prior is the voice's pitch: a sung note only has energy at its harmonics. But a pitch tracker run on the whole mix follows whatever is loudest, and on the evaluation song it found the singer in only 3% of frames. Three cues fix most of that. The first is that lead vocals sit dead centre, so bins where left and right agree are kept and the rest are suppressed.

The register also gives a voice-type estimate: it said “alto” and “tenor” for the two test voices, both right. (The median of the evaluation song's line sits near the alto/soprano border, so the browser port, whose track differs in under 1% of frames, says “soprano” there.) Building this also exposed a bug in Pianoify: its octave correction pulled high voices down an octave, because the even harmonics of a high note make the octave below look salient. It now moves down only when the lower pitch shows its own fundamental or third harmonic. On an isolated high voice that took pitch accuracy from 33% to 95%. Both the Python and browser versions of Pianoify are fixed.

Since the data cannot tell the parts apart, each part gets a penalty that describes what it looks like. Sustained chords barely change over time. Drums are short and broadband, so only a few frames are active. The voice only exists on a comb M at the harmonics of the tracked pitch, and is sparse there.

The harmonic penalty needs a trick. A steady tone's STFT coefficient still rotates in phase from frame to frame, so its time difference is not small. The operator E undoes that rotation first, as in Masuyama, Yatabe and Oikawa (2019).

The problem is convex but every term is non-smooth, so it is solved with the primal-dual hybrid gradient method (Chambolle–Pock). Each penalty and the data constraint gets its own dual variable, and each dual update is a cheap projection. The method converges when the step sizes satisfy τ Σ σ_i‖K_i‖² < 1, where K_i is the operator inside term i.

The first version said “the voice's STFT is zero off the comb”. That returned an almost silent voice: in a 4× redundant frame, onsets and window side-lobes always leak off the comb, so almost no real signal qualifies. The fix is to build the voice from coefficients that live on the comb, v = F^T(M ⊙ z). Any z gives a valid voice.

The ℓ1 penalty also shrinks the voice, which came out 3–4 dB too quiet. A final step uses the three recovered parts only as soft masks and re-splits the measured recording with them.

Round 1 treated the decoded audio as exact. But an MP3 doesn't store a waveform. It stores rounded MDCT coefficients, so each one really says “the original coefficient was somewhere in this interval”. That is a much more honest measurement, and closer to the tomography setting: g is no longer one point but a box. Reading it needs two things: the decoder's synthesis as a linear operator S, and the interval each stored integer stands for.

The decoder is minimp3 (public domain) with one hook added after its Huffman decoding. For every granule and channel the hook writes the integer quantised value q and the quantiser step of every coefficient, with the block type and the mid/side flag. That is everything the file says about the signal.

Layer III synthesis is linear: short-block reorder, anti-alias butterflies, a 36- or 12-point IMDCT with a window per block type and overlap-add, a sign flip on odd subbands, then the 32-band polyphase filterbank. The IMDCT/window/overlap stage and the polyphase stage are measured as impulse responses of minimp3 itself, one set per pair of consecutive block types. So S reproduces minimp3's output to 3·10⁻⁷, and its transpose is exact to 10⁻¹⁵. The encoder's analysis has to be S⁻¹. S turns out to be orthogonal up to a per-coefficient scale (long and short blocks have different IMDCT gains, and mid/side frames carry a factor of 2), so the analysis is a weighted adjoint, A = W S^T, with W read off the measured responses. It inverts S to 84 dB.

The test that matters is whether the clean mix, before encoding, analysed with A, lands in the cells the file specifies. It exposed two facts about the encoder:

• The encoder scaled its input by 0.950. Every decoded coefficient is 0.950 times the clean analysis, with correlation 1.0000 in every band. This is also why round 1 measured the codec as only 23.6 dB clean: most of that error was this gain, not quantisation.
• LAME rounds to the nearest dequantised level. So the cell for q runs between the midpoints of neighbouring levels |q|^4/3·step. At q = 0 and 1 this is exactly the ISO offset 0.4054 = 1 − 0.5^3/4. That is the code above.

Not every cell is information. A zero above the encoder's lowpass, or in a band the encoder dropped as a whole, means “not stored”, not “close to zero”. Those coefficients are left free.

Leaving the uninformative cells out, the clean mix lands inside its cell for almost every coefficient:

In the solver this only changes the data term's prox. Equality with g becomes a box, Gram^½·lo ≤ Ψ(v + h_c + p_c) ≤ Gram^½·hi, with Ψ = Gram^−½S^T (so ΨΨ^T = I and ‖Ψ‖ = 1): the dual update clips the encoder-side analysis of the current mix into its cells.

The voice's round-1 prior is only “sparse on the comb”. A sung note is also smooth over time, but not along a fixed frequency: vibrato and glides bend every harmonic. In a constant-Q frame C (a painless non-stationary Gabor frame, 36 bins per octave from 80 Hz to 4 kHz, a common 3 ms hop, Parseval on its band) a pitch change moves every harmonic by the same number of bins. So the total variation ‖D_θCv‖₁ can step along the pitch slope, dk(t) = 36·Δlog₂f₀, instead of straight across. Both operators check out numerically: the adjoints are exact to 10⁻¹⁵, and ‖D_θC‖ = 1.93.

That phase prediction is the hard part. On a synthetic vibrato tone over a steady chord, following the pitch with the phase predicted from the harmonics h·f₀ makes the voice 3.8× cheaper than the chord. On the test songs it doesn't, and the reason is precision. Over a 3 ms hop the predicted phase is only right if the instantaneous frequency is right to a few cents. A pitch tracker gives 16 ms frames on a ⅓-semitone grid and smooths vibrato away, so the prior ends up rewarding the accompaniment, whose steady partials predict their own phase well. The fix is an outer loop: solve once without the TV, then use the recovered voice's own pitch and phase to define the TV, and solve again.

The third idea, a range gate, is one function in run2.py: pitch frames outside the line's own register (5th–95th percentile, widened by a quarter) are dropped before building the comb. What all three were worth is under Results, round 2.

The tool on the project page runs the round-1 method: centre extraction, vibrato-weighted pitch tracking and register folding from voiceline.py, then separate and refit from unmix.py, with the same frame (1024-sample sqrt-Hann, hop 256), the same comb and the same tuned weights (α_H = 3s, α_P = 9s·N^¼, α_V = 0.3s, with s the median magnitude of the mixture's spectrum). It is plain JavaScript with no dependencies: web/unmix.js is the maths, web/worker.js runs it off the main thread, web/app.js is the page. Nothing is uploaded.

The browser decodes whatever file it can play, and its own resampler brings it to 16 kHz, the rate the method was tuned at.

Speed comes from two changes that leave the maths exactly as it was. Every signal is real, so every spectrum the solver touches is Hermitian: only bins 0–512 are stored, and two real frames share one complex FFT (a radix-4 FFT on Float64Arrays, results kept in Float32Arrays like Python's complex64).

The second change: with exact data equality (ε = 0, the default), the data dual only ever accumulates F(v + h_c + p_c − g_c), and F^TF = I. So it can live in the time domain, which removes three transforms per iteration. Round 2's separate2 does the same.

The pitch tracker had to be rebuilt without librosa. Its constant-Q transform is emulated directly at 16 kHz: each of the 252 bins gets a Hann filter of length Q·sr/f (Q = 51.9, librosa's), applied as a sparse kernel on one 16384-point FFT per pair of columns, on librosa's 22050/512 time grid with its √length scaling. The 4–8 Hz band-pass for the vibrato score uses scipy's Butterworth coefficients, with the same odd padding and initial conditions as sosfiltfilt; the median and box filters use scipy's “reflect” edges; and the Viterbi search is the full search librosa does, ties included. The worker runs the three steps in order and posts progress, so the page can show it and stop it.

Both versions were run on the same 16 kHz excerpts of the two synthetic songs (20 s from 2.5 s, the MP3 decoded by ffmpeg and aligned as in run.py), scored against the true parts with the metrics in run.py. Given the same pitch track, the JavaScript solver and refit reproduce Python's output to 126 dB SDR, so the remaining differences come from the pitch tracker, whose track matches Python's in 99.3% (evaluation) and 99.4% (tuning) of frames.

Tests use synthetic songs whose parts are known exactly. The voice is formant-synthesised (glottal pulses through vowel resonances, with vibrato, glides, jitter and breath noise). The chords are panned across the stereo field, the bass is centred, and the drums are spread. The mix goes through a real 128 kbps MP3 encoder. Every setting was chosen on a tuning song: a fifth lower, tenor voice, different vowels, panning and random seed. The numbers below are from a separate evaluation song (alto voice), 20 seconds, 300 iterations. The demo song on the project page is that same 20-second excerpt.

What helped, and what didn't

• The voice cues matter most (section 03): −4.8 dB voice SI-SDR with Pianoify's plain tracker, +4.4 dB with the cues.
• Stereo is worth about 1–2 dB. Treating the song as mono costs 1.3 dB on the voice and 2 dB on the accompaniment.
• The reconstruction beats masking with the same pitch: +1.2 dB voice SDR and +1.6 dB accompaniment. It also splits the accompaniment into harmonic (7.4 dB) and drums (1.5 dB), and its voice is much cleaner to pitch-track (92% vs. 67% with true pitch).
• The split data fidelity did nothing here (section 05): exact equality, bands and tolerances all landed within 0.1 dB.
• Pitch tracking is the ceiling. With the true pitch, quality roughly doubles in dB. Around 5–6 s in the spectrograms below, the tracker locks onto chord tones and they leak into the voice.

Round 2: the MP3 as the measurement, and a prior that follows the pitch

The first round left three ideas open, each from the tomography analogy: the MP3's quantiser cells as the data term (section 07), directional TV that follows the pitch (section 08), and range as a constraint, dropping pitch frames outside the singer's own register from the comb. All three were tested on the same held-out song, with settings chosen on the tuning song, 20 s, 300 iterations.

The conclusion of round 1 stands, sharper: pitch tracking is the bottleneck. Every extension that helps does so mostly when the pitch is right.

Python with numpy, scipy and librosa, plus ffmpeg for the MP3. It uses ../pianoify/pianoify.py for the pitch salience, so keep the two folders side by side.

The browser code also runs in Node (require('./web/unmix.js') exposes track, separate and refit), which is how it was checked against the Python on the same excerpts.