1
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – hotpaw2
    Jun 17, 2019 at 9:01
  • $\begingroup$ 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) $\endgroup$
    – audioboy
    Jun 17, 2019 at 9:29
  • $\begingroup$ Update: dsprelated.com/showarticle/1284.php $\endgroup$
    – Cedron Dawg
    Jul 17, 2019 at 15:32

3 Answers 3

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ More details in the comments to answers to this question: dsp.stackexchange.com/questions/3301/… $\endgroup$
    – hotpaw2
    Jun 17, 2019 at 9:19
  • $\begingroup$ 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? $\endgroup$
    – audioboy
    Jun 17, 2019 at 9:37
0
$\begingroup$

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

Improve this answer
$\endgroup$
0
$\begingroup$

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()
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.