This research develops efficient solution methods for linear systems with scalar and multi-level Toeplitz structure. Toeplitz systems are common in one-dimensional signal-processing applications, and typically correspond to temporal- or spatial-invariance in the underlying physical phenomenon. Over time, a number of algorithms have been developed to solve these systems economically by exploiting their structure. These developments began with the Levinson-Durbin recursion, a classical fast method for solving Toeplitz systems that has become a standard algorithm in signal processing. Over time, more advanced routines known as superfast algorithms were introduced that are capable of solving Toeplitz systems with even lower asymptotic complexity. For multi-dimensional signals, temporally- and spatially-invariant systems have linear-algebraic descriptions characterized by multi-level Toeplitz matrices, which exhibit Toeplitz structure on multiple levels. These matrices lack the same algebraic properties and structural simplicity of their scalar analogs. As a result, it has proven exceedingly difficult to extend the existing scalar Toeplitz algorithms for their treatment. This research presents algorithms to solve scalar and two-level Toeplitz systems through a constructive approach, using methods devised for specialized cases to build more general solution methods. These methods extend known scalar Toeplitz inversion results to more general scalar least-squares problems and to multi-level Toeplitz problems. The resulting algorithms have the potential to provide substantial computational gains for a large class of problems in signal processing, such as image deconvolution, non-uniform resampling, and the reconstruction of spatial volumes from non-uniform Fourier samples.