A second order energy stable numerical scheme is presented for the two and three dimensional Cahn-Hilliard equation, with Fourier pseudo-spectral approximation in space. The convex splitting nature assures its unique solvability and unconditional energy stability. Meanwhile, the implicit treatment of the nonlinear term makes a direct nonlinear solver not available, due to the global nature of the pseudo-spectral spatial discretization. In turn, a linear iteration algorithm is proposed 
to overcome this difficulty, in which a Douglas-Dupont-type regularization term is introduced. As a consequence, the numerical efficiency has been greatly improved, since the highly nonlinear system can be decomposed as an iteration of purely linear solvers. Moreover, a careful nonlinear analysis shows a contraction mapping property of this linear iteration, In addition, a maximum norm bound of numerical solution is also derived at a theoretical level. A few numerical examples 
are also presented in this talk.