This paper presents a new consistent splitting scheme for the numerical solution of incompressible Navier-Stokes flows; allowing to consistently decouple the computation of velocity and pressure. The scheme is not a pressure-correction or velocity-correction scheme, and does not display the splitting error in pressure associated with these fractional step methods. The (linearized) momentum equations are first solved based on an explicit treatment of the pressure, resulting in an advection-diffusion problem for each velocity component. A least-squares projection is used to numerically solve the advection-diffusion problem. Next, a div-curl problem is solved to make the velocity field solenoidal. This step is also handled by a least-squares projection. Finally, a pressure Poisson problem is solved to obtain the pressure field induced by the solenoidal velocity field. This is done by solving the weak Poisson problem by a Galerkin projection or alternatively by solving the strong Poisson problem by a least-squares projection. At each stage we only see coefficient matrices with a symmetric positive definite structure, and use matrix-free (preconditioned) conjugate gradient methods to numerically solve for the velocity and pressure fields. High-order C 0 spectral basis are used to span the finite element spaces. A verification benchmark shows optimal algebraic convergence rates in time for the velocity, pressure, and vorticity. The scheme is further verified by simulating the two-dimensional unsteady flow past a circular cylinder up to moderately high Reynolds numbers.