An O(n log n) direct solution method, based on circulant and skew-circulant representation (CSR) for Toeplitz matrix inversion, where n is the total number of spatial grid points, is proposed for solving one-dimensional (1D) and two-dimensional (2D) tempered fractional diffusion equations with constant coefficients. It is proved that the CSR for the inverse of the coefficient matrix for 1D problems can be obtained in O(n log n) operations with circulant preconditioner. For 2D problems, where the total number of spatial grid points n = nxny, the operation cost to obtain the CSR is O(p log p), where p = max{nx, ny}, which is negligible in the whole process of solving the tempered fractional diffusion equation. Numerical experiments reveal that the proposed method is second order accurate and the computations are very efficient.