Numerical solutions of the multi-dimensional spatial fractional Allen-Cahn equations are studied. After the semi-discretization of the spatial fractional Riesz derivative, a system of nonlinear ordinary differential equations with the Toeplitz structure is obtained. In order to reduce the computational complexity, a two-level Strang splitting method is proposed, where the Toeplitz matrix in the system is represented as the sum of circulant and skew-circulant matrices. Therefore, the method can be quickly implemented by the fast Fourier transform, avoiding expensive Toeplitz matrix exponential calculations. It is shown that the discrete maximum principle of the method is unconditionally preserved. Moreover, the analysis of errors in the infinite norm with second-order accuracy is carried out in both time and space. Numerical tests support the theoretical findings and show the efficiency of the method.