We present a generalized pulse-echo speed-of-sound (SoS) tomography method based on the matrix imaging formalism. Our approach begins by applying the Fourier split-step technique to back-propagate the recorded wavefield in depth. The resulting wavefields are then used to synthesize tissue's response to various combinations of emission and reception angles and to construct a distortion matrix that captures wavefront distortions. Matrix elements sharing a common mid-angle exhibit strong correlations, which we exploit via an iterative phase reversal process to estimate phase aberrations locally. Finally, we linearly relate the estimated phases to tissue SoS using the straight-ray approximation and solve a Tikhonov-regularized inversion to infer the SoS distribution. Because the Fourier split-step method supports heterogeneous media, the process is repeated iteratively to refine the SoS estimates until convergence is reached. Numerical and in vivo examples demonstrate that our method achieves higher accuracy, improved resolution, and increased robustness against the initially assumed SoS.