First a demonstration that the squares of both
$$\begin{align}&[\dots, 0, 0, 1,\hphantom{-}1, 0, 0, \dots] \text{ and}\\
&[\dots, 0, 0, 1, -1, 0, 0, \dots]\end{align}$$
equal
$$[\dots, 0, 0, 1, \hphantom{-}1, 0, 0, \dots]\hphantom{\text{ and}}$$
but the squares of their sinc interpolations differ (Fig. 1):
Figure 1. Squares of sinc interpolations of $[1, 1]$ (blue) and $[1, -1]$ (red).
This demonstrates that it is not in general possible to recover the square of a band-limited signal from its uniform samples taken at the critical sampling frequency of the band-limited signal.
Let's test different interpolation approaches in Octave. The gold standard of band-limited interpolation is DFT–zero-pad–DFT:
>> format free
>> x = [1 2 3];
>> y = [1 2 3 0 0 0];
>> abs(fft(x))
ans =
6 1.73205 1.73205
>> abs(fft(y))
ans =
6 4.3589 1.73205 2 1.73205 4.3589
The last set of numbers are the magnitude values calculated from perfectly interpolated frequency domain bins. Let's try interpolating magnitude instead:
>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)))), [0])))
ans =
6 4.57735 1.73205 0.309401 1.73205 4.57735
Seems quite a bit off. Now let's try interpolating squared magnitude:
>> fft(fftshift(horzcat([0 0], fftshift(ifft(abs(fft(x)).^2)), [0])))
ans =
36 25 3 -8 3 25
>> sqrt(ans)
ans =
(6,0) (5,0) (1.73205,0) (0,2.82843) (1.73205,0) (5,0)
It's not only off, but also sqrt() returned a complex number for a negative interpolated value. So is the only valid way to interpolate the bin values?
Let's give this one more chance, by trying to interpolate frequency domain data upsampled by a factor of 2. Because of the even length transform, this requires duplicating the "Nyquist sample", so apologies if the code is getting hard to read.
>> z = [1 2 3 0 0 0 0 0 0 0 0 0];
>> abs(fft(z))
ans =
6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485
The above is what we want. Let's try interpolating 2x upsampled magnitude:
>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)))), [0 0 0]))
ans =
3.36365 1.10447 0.318175 0 0 0 0 0 0 -0.208949 0.318175 1.10447
>> ans(4) = ans(length(ans)-2)
ans =
3.36365 1.10447 0.318175 -0.208949 0 0 0 0 0 -0.208949 0.318175 1.10447
>> fft(ans)
ans =
5.79105 5.59483 4.56785 2.7273 1.5231 1.76882
2.20895 1.76882 1.5231 2.7273 4.56785 5.59483
That is still off. Let's try interpolating 2x upsampled squared magnitude:
>> fftshift(horzcat([0 0 0], fftshift(ifft(abs(fft(y)).^2)), [0 0 0]))
ans =
14 8 3 0 0 0 0 0 0 1.18424e-15 3 8
>> ans(4) = ans(length(ans)-2)
ans =
14 8 3 1.18424e-15 0 0 0 0 0 1.18424e-15 3 8
>> sqrt(fft(ans))
ans =
6 5.55485 4.3589 2.82843 1.73205 1.77302 2 1.77302 1.73205 2.82843 4.3589 5.55485
Now it works perfect! The take home message is to upsample (time domain zero pad) at least by a factor of two before trying to interpolate in the frequency domain, and to interpolate squared magnitude rather than magnitude. It works because taking the squared magnitude is the same as multiplying each bin value by its complex conjugate. Complex conjugation preserves the bandwidth of the bandlimited function represented by the data, so the multiplication doubles the "time domain bandwidth" because it is equivalent to time domain convolution. Note that when choosing the interpolation method the 2x upsampled squared magnitude is still critically sampled, so further oversampling should make accurate interpolation much easier.
I about answered my own question, but will accept further insight as an answer!
P.S. Just found out there is also interpft, which does the interpolation with less syntax.
Utilizing additional information
Interpolation becomes easier or even exact with additional information about the data, for example that it is a critically sampled squared time-shifted sinc. In that case, given the two samples $\alpha$ just before and $\gamma$ just after the largest sample, the time $-1<d<1$ of the peak can be calculated by:
$$d = \begin{cases}0&\text{if }\alpha = \gamma,\\
\frac{1 - \sqrt{1 - p^2}}{p},\, p = \frac{γ - α}{α + γ}&\text{otherwise,}\end{cases}$$
with $d = 0$ meaning that the sinc is shifted to exactly the time of the largest sample. The same interpolation can be done at half sample rate, which is the critical sampling frequency of the underlying function the magnitude of which equals that of a time-shifted sinc, with the two successive largest-valued samples $\alpha$ and $\beta$ of the square of the absolute value of the underlying function:
$$0 \le d \le 1\\
d = \begin{cases}0&\text{if }\beta = 0,\\
1&\text{if }\alpha = 0,\\
\frac{1 + p - \sqrt{1 - p^2}}{2p},\, p = \frac{\beta - α}{α + \beta}&\text{otherwise,}\end{cases}$$
The formulas are not affected by amplitude scaling of the data. For frequency estimation you'd rarely have a purely real time-shifted sinc, but if you do, there are exact interpolation formulas for it.
