Abstract Binary black hole (BBH) systems detected via gravitational-wave emission are a recently opened astrophysical frontier with many unknowns and uncertainties. Accurate reconstruction of the binary distribution with as few assumptions as possible is desirable for inference of formation channels and environments. Most population analyses have, though, assumed a power law in binary mass ratio q, and/or assumed a universal q distribution regardless of primary mass. Methods based on kernel density estimation allow us to dispense with such assumptions and directly estimate the joint binary mass distribution. We deploy a self-consistent iterative method to estimate this full BBH mass distribution, finding local maxima in primary mass consistent with previous investigations and a secondary mass distribution with a partly independent structure, inconsistent both with a power law and with a constant function of q. We find a weaker preference for near-equal-mass binaries than in most previous investigations; instead, the secondary mass has its own “spectral lines” at slightly lower values than the primary, and we observe an anticorrelation between primary and secondary masses around the ∼10 M ⊙ peak.