Context. Oscillations in rotating stars with frequency (sigma) over bar of the same order or smaller than the rotation rate Omega cannot be described by a single spherical harmonic due to the effect of the Coriolis force. This is a serious complication which is usually treated by writing the eigenfunctions as a (truncated) sum of different spherical harmonic degrees for a given m-value, or by neglecting the theta part of the rotation vector (Traditional Approximation). Aims. We aim for a more adequate treatment of the coupling with higher angular degrees for low frequency oscillations of rotating stars by taking the Coriolis force fully into account, so that the coupling is included in the analysis (up to the grid resolution) and to compare the results with analyses based on the above mentioned approximations. Methods. To this end a new, more efficient version of a 2D(r, theta) implicit oscillation code was developed in which no a priori assumptions about the theta variation of the eigenfunctions is made, enabling a better treatment of the rotational truncation of g-modes near the stellar poles, and for which similar to 150 gridpoints in theta are feasible. We test the code by comparing the simulated oscillation spectrum with that obtained by the MOST photometry for the similar or equal to 6 M-circle dot Be star HD 163868. Results. We find both prograde and retrograde overstable modes (although more prograde than retrograde modes) and confirm the existence of low degree odd r-modes, destabilised by the kappa-mechanism. The ultra-low frequency modes that could not be explained in a previous analysis are interpretated as (high degree) retrograde m = 1 modes with vertical bar(sigma) over bar vertical bar similar to Omega i.e. sigma similar to 0. A reasonably good fit the observed oscillation spectrum is possible if we assume that only even modes are observed (no unstable r-modes visible). This requires a nearly equator-on view of the observed star, consistent with the measured high nu sin i value of 250 km s(-1).