Helmholtz equation is widely applied in the scientific and engineering problem. For the solution of the threedimensional Helmholtz equation, the computational efficiency of the algorithm is especially important. In this paper, in order to solve the contradiction between accuracy and efficiency, a fast high order finite difference method is proposed for solving the three-dimensional Helmholtz equation. First, a traditional fourth order method is constructed. Then, fast Fourier transformation are introduced to generate a block-tridiagonal structure which can easily divide the original problem into small and independent subsystems. For large 3D problems, the computation of traditional discrete Fourier transformation is less efficient, and the memory requirements increase rapidly. To fix this problem, the transformed coefficient matrix is constructed as a sparse structure. In light of the sparsity, the algorithm presented in this paper requires less memory space and computational cost. This sparse structure also leads to independent solving procedure of any plane in the domain. Therefore, parallel method can be used to solve the Helmholtz equation with large grid number. In the end, three numerical experiments are presented to verify the effectiveness of the fast fourth-order algorithm, and the acceleration effect to use the parallel method has been demonstrated.