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

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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

• Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

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