I am aware that this is an old thread, but since I had some trouble with the exact same problem and found a different solution. looking at the output of the fft I noticed that a frequency shift from one bin to the next a phase shift occurs in the bins, correcting for this phaseshift gives quite a good result.
I would like to share it for future visitors.
Example code in python with numpy as np.
# I use an approximate confined gaussian window:
fs = 1000.0 # samplerate
N = 1024 # fft size
L = N + 1
sigma = 0.1
G = lambda x : np.exp(-np.power((x - N/2)/(2*L*sigma), 2))
w = lambda n : G(n) - (G(-0.5)*(G(n+L) + G(n-L)))/(G(-0.5+L)+G(-0.5-L))
x = np.arange(0,N)
window = w(x)
# some preparations
x_range = np.arange(0, N) / 1000.0 * 2 * math.pi
f = 10
inner_pad = np.zeros(N)
signal = np.cos(x_range * f + 2.4)
# fft calculation
mean = np.mean(signal)
windowed = (signal - mean) * window # multiply by the window function
padded = np.append(windowed, inner_pad) # add 0s to double the length of the data
spectrum = np.fft.fft(windowed) / N # take the Fourier Transform and scale by the number of samples
spectrum = spectrum[0:N//2] * 2 # remove mirrored half and scale
spectrum_abs = np.abs(spectrum)
# fft interpolation
# based on http://www.add.ece.ufl.edu/4511/references/ImprovingFFTResoltuion.pdf
max_index = np.argmax(spectrum_abs)
angle = np.angle(spectrum[max_index])
prev = spectrum_abs[max_index - 1]
maxx = spectrum_abs[max_index]
next = spectrum_abs[max_index + 1]
numerator = np.log(next / prev)
denominator = 2 * np.log(np.power(maxx, 2) / (prev * next))
# calculate bin offset from max_index
delta = numerator / denominator
# calculate exact frequency
f_exact = (max_index + delta) / (N / fs)
# calculate phase shift as result of frequency delta
delta_phase_shift = delta * math.pi
# shifting the result to be between 0 and 2pi
phase_shift = np.fmod(angle - delta_phase_shift + 2 * math.pi, 2 * math.pi)
# output result (first 100 samples)
fig, ax = plt.subplots()
ax.plot(signal[:100], "b-")
signal_recovered = np.cos(x_range * f_exact + phase_shift)
ax.plot(signal_recovered[:100], "r.")
plt.show()