We present fast fourth-order finite difference scheme for 3D Helmholtz equation with Neumann boundary condition. We employ the discrete Fourier transform operator and divide the problem into some independent subproblems. By means of the Gaussian elimination in the vertical direction, the problem is reduced into a small system on the top layer of the domain. The procedure for solving the numerical solutions is accelerated by the sparsity of Fourier operator under the space complexity of O(M3). Furthermore, the method makes it possible to solve the 3D Helmholtz equation with large grid number. The accuracy and efficiency of the method are validated by two test examples which have exact solutions.