  High frequency quasi-periodic oscillations (QPOs) from weakly magnetized
neutron stars display rapid frequency variability and high coherence with
quality factors up to at least 200 at frequencies around 850 Hz. Their
parameters have been estimated so far from standard min(chi2) fitting
techniques, after combining a large number of Power Density Spectra (PDS), as
to have the powers normally distributed. Accounting for the statistical
properties of PDS, we apply a maximum likelihood method to derive the QPO
parameters in the non Gaussian regime. The method presented is general, easy to
implement and can be applied to fitting individual PDS, several PDS
simultaneously or their average, and is obviously not specific to the analysis
of kHz QPO data. It applies to the analysis of any PDS optimized in frequency
resolution and for low frequency variability or PDS containing features whose
parameters vary on short timescales, as is the case for kHz QPOs. It is
equivalent to the standard chi^2 minimization fitting when the number of PDS
fitted is large. The accuracy, reliability and superiority of the method is
demonstrated with simulations of synthetic PDS. We show that the maximum
likelihood estimates of the QPO parameters are asymptotically unbiased, and
have negligible bias when the QPO is reasonably well detected. By contrast, we
show that the standard min(chi2) fitting method gives biased parameters with
larger uncertainties. The maximum likelihood fitting method is applied to a
subset of archival Rossi X-ray Timing Explorer (RXTE) data of the neutron star
X-ray binary 4U1608-522. We show that the kHz QPO parameters can be measured on
8 second timescales and that the time evolution of the frequency is consistent
with a random walk. This enables us to estimate the intrinsic quality factor of
the QPO to be around 260, whereas previous analysis indicated a maximum value
around 200 (abridged).
