We consider in this paper numerical approximations of a phase field model for twophase ferrofluids, which consists of the Navier-Stokes equations, the Cahn-Hilliard equation, the magnetostatic equations, and the magnetic field equation. By combining the projection method for the Navier-Stokes equations and some subtle implicit-explicit treatments for coupled nonlinear terms, we construct a linear, decoupled, fully discrete finite element scheme to solve the highly nonlinear and coupled multiphysics system efficiently. The scheme is provably unconditionally energy stable and leads to a series of decoupled linear equations to solve at each time step. Through numerous numerical examples in simulating benchmark problems such as the Rosensweig instability and droplet deformation, we demonstrate the stability and accuracy of the numerical scheme.