The Wigner functions are quantum analogues of the classical phase space density. They are functions of both momenta and coordinates and are typically obtained from wavefunctions that are usually given as functions of only either coordinates or momenta. The motivation for this work was the question whether using the Wigner functions one can perform the formalism of quantum mechanics in the full phase space without referring to the actual wavefunction. While the formal equations for such a procedure are well known in the form of Moyal brackets, these can hardly be used directly as, in the general case, they represent a differential equation of an infinite order. While a naive semiclassical expansion of the Wigner function does make the equations in the Moyal approach to become of finite order for any order of the Planck's constant, such an approach fails immediately after the lowest order solution, which itself is the delta distribution on the invariant component of the classical phase space of the system. The Moyal bracket can, however, be rewritten in the form of an integral equation. In analogy to the usual WKB semiclassical approach for the wavefunction we expand the phase of the Wigner function in terms of the powers of the Planck's constant. If we now perform the analysis of the integral equation around the stationary point as given by the lowest order contribution to the phase of the Wigner function, we obtain a set of recursive equations that, at each order of the Planck's constant, contain only a finite number of terms. Unlike in the previously mentioned naive Moyal approach this set of equations allows us to find a solution to all orders of the Planck's constant. One should also note that such a separation does not happen if one inserts the phase expansion directly into the Moyal bracket. We apply this approach to the harmonic oscillator. We find the energy eigenvalues to all orders of the Planck's constant and hence obtain them exactly, while for the Wigner functions themselves we show that even the low order solutions match the exact ones very nicely and that they also agree with certain analytical approximations for Legendre polynomials.