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).