This thesis concerns preconditioning for Toeplitz-related systems. Specifically, we consider functions of Toeplitz matrices, i.e. h(T_n)x_n=b_n, where h(z) is an analytic function and T_nϵℂ^n×n is a Toeplitz matrix. We propose the absolute value circulant matrix |h(C_n)| as a preconditioner for h(T_n), where C_nϵℂ^n×n is the optimal circulant preconditioner, the superoptimal circulant preconditioner, or Strang's circulant preconditioner derived from T_n, and show that |h(C_n)|^-1h(T_n) has clustered spectra that account for the effectiveness of such preconditioner. When h(T_n) is a real matrix, we can first premultiply it by the anti-identity matrix Y_nϵℝ^n×n to obtain a (real) symmetric matrix Y_nh(T_n) without normalizing the original matrix. To ensure |h(C_n)| is an effective preconditioner for Y_nh(T_n), we show that |h(C_n)|^-1Y_nh(T_n) has clustered spectra around ±1. As Y_nhT_n) is symmetric yet possibly indefinite, we can use the minimal residual method for the corresponding linear system with guaranteed convergence that depends only on its eigenvalues. We further show that the ideas of symmetrization and absolute value preconditioning for Toeplitz systems can be extended to the block Toeplitz matrix case. An application on time-stepping methods for evolutionary ordinary/partial differential equation problems is also discussed. Numerical results are given to demonstrate the effectiveness of our proposed preconditioners.