The asker has a 1D sampled signal containing periodic peaks that are narrower than the sampling interval. They want to locate each peak's true position with sub-sample accuracy better than 0.1 samples. They've already got a stable frequency estimate from autocorrelation, and they're scanning phase on a grid to refine. The trouble: with only a few cycles available and peaks narrower than Δt, naive interpolation around the discrete maximum is unreliable.

Why this is genuinely hard. The classical sub-sample peak tricks — parabolic (quadratic) interpolation across the three samples nearest the maximum, or Gaussian-log interpolation — assume the underlying peak is smooth and well-sampled. When the peak's full-width-at-half-maximum is below the sampling period, those three samples don't bracket the actual peak; they're just three points on the slope of an aliased reconstruction. The result is heavy bias that depends on where the true peak sits relative to the sample grid (the classic "picket fence" effect). Compounding this: only a few cycles means you can't lean on long-baseline phase estimation from a clean FFT either, since the frequency bin width is large and leakage dominates.

A direction that actually works. Stop treating peak localization as a peak-finding problem and treat it as a model-fitting problem on the whole signal at once:

• Assume a parametric model: s(t) = Σ A·g(t − t₀ − k/f) + noise, where g is the (known or assumed) peak shape and t₀, f, A are the unknowns.
• Seed f from autocorrelation and t₀ from the coarse phase scan, then non-linearly optimize (Levenberg–Marquardt via scipy.optimize.least_squares) over all samples — including the noisy baseline. Every cycle contributes information, so few cycles still works.
• If the peak shape is unknown, parameterize it (Gaussian width σ, or a sinc if the peak is band-limited by an anti-alias filter) and fit that too.

Gotchas. If there's no anti-alias filter before sampling, peaks narrower than Δt are aliased — the recorded waveform isn't a faithful low-passed version of the truth, and no amount of post-hoc fitting recovers the lost information. The accuracy floor is then set by the noise on the baseline, not the peak shape. Also: if peaks are truly impulse-like, the sampled amplitude depends strongly on where the peak fell within the sample interval, which biases amplitude-weighted estimators. Use phase-of-FFT-at-fundamental as a cross-check — it has the lowest variance estimator for periodic signals in white noise (CRLB) and sidesteps the peak-shape question entirely.