# Interpolation of estimated frequency between fft bins: what about phase?

I am developing an additive re-synthesier, which involves identifying partials from FFT frames.

Partial identification requires amplitude peaks. The simplest way to do this is to just look at the amplitude of each bin, but this of course can miss partials that happen to lay between bins, and there energy is spread over the two adjacent bins.

I understand there are various interpolation methods (of varying complexity) to estimate the actual amplitude peak and frequency... however, if taking such an approach, how to estimate the phase offset? Is it even possible?

• Interpolation for analysis-resynthesis is often done phase first, rather than magnitude first. The interpolated spectrum peaks are often estimated by using phase vocoder methods, looking at the phase change between adjacent overlapped windows. Jun 17, 2019 at 9:01
• Ok thanks. I already know of the method to estimate the bin frequency by comparing adjacent frame phase offsets (and comparing it to the expected change)... however, with that approach, I cant then estimate the corresponding amplitude... or can I? (That is why I was looking into the interpolation approach above) Jun 17, 2019 at 9:29
• Jul 17, 2019 at 15:32

If you do an fftshift before the FFT (to center the phase reference in the data frame), you can use the same interpolation methods on the real (even) and imaginary (positive odd) components in the complex FFT result, and thus get interpolated phase estimates (referenced to the original data frame’s center) between FFT result bins.

• More details in the comments to answers to this question: dsp.stackexchange.com/questions/3301/… Jun 17, 2019 at 9:19
• Ah ok, so do you mean then to do first do regular short term FFTs, and estimate my exact frequencies by comparing phase offsets of adjacent phases, and then secondly, do another round of DFTs at those exact frequencies, to get the actual phase and amplitude? Alternatively, a cheaper idea would be to do my stFFTs, estimate exact frequencies, then use those estimated frequencies to interpolate between my bins and get the corresponding amplitudes.. does that make sense? But In this case I would just use the phase from the FFT, and not adjust it? Jun 17, 2019 at 9:37

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
signal = np.cos(x_range * f + 2.4)

# fft calculation
mean = np.mean(signal)
windowed = (signal - mean) * window                       # multiply by the window function
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
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()


For a single pure real noiseless tone, the technique given in my blog article

will give you an exact answer.

For tone peaks that aren't too close to other tones, (e.g. a bin or two), and not too much noise, it will give you a very good estimate.

Mathematically, it is equivalent to constructing two basis vectors at your estimated frequency and solving for the phase like a DFT bin does, but it does it from the neighboring bin values rather than the source values so it is extremely efficient in comparison. (And frame sample count independent)

Ced