To investigate the origins of S0 galaxies, we present a new method of analysing their stellar kinematics from discrete tracers such as planetary nebulae. This method involves binning the data in the radial direction so as to extract the most general possible non-parametric kinematic profiles, and using a maximum-likelihood fit within each bin in order to make full use of the information in the discrete kinematic tracers. Both disc and spheroid kinematic components are fitted, with a two-dimensional decomposition of imaging data used to attribute to each tracer a probability of membership in the separate components. Likelihood clipping also allows us to identify objects whose properties are not consistent with the adopted model, rendering the technique robust against contaminants and able to identify additional kinematic features. The method is first tested on an N-body simulated galaxy to assess possible sources of systematic error associated with the structural and kinematic decomposition, which are found to be small. It is then applied to the S0 system NGC 1023, for which a planetary nebula catalogue has already been released and analysed by Noordermer et al. The correct inclusion of the spheroidal component allows us to show that, contrary to previous claims, the stellar kinematics of this galaxy are indistinguishable from those of a normal spiral galaxy, indicating that it may have evolved directly from such a system via gas stripping or secular evolution. The method also successfully identifies a population of outliers whose kinematics are different from those of the main galaxy; these objects can be identified with a stellar stream associated with the companion galaxy NGC 1023A.