# Interpolation of magnitude of discrete Fourier transform (DFT)

For example for peak frequency finding, it seems valid to use band-limited interpolation methods on the complex DFT bins, or separately on their real and imaginary parts and to calculate the magnitudes or squared magnitudes of the results. But how about band-limited interpolation of the magnitudes of the bins (I don't think that's valid), or their squared magnitudes (maybe valid)? By valid I mean that perfectly interpolated values should be equal to those found by calculating them from a larger DFT of a zero-padded version of the time-domain signal.

The first approach guarantees a non-negative result, unlike the others if the interpolation is not perfect, see this question about nonnegative or positive band-limited interpolation.

• hey Olli, is this a question? – robert bristow-johnson Oct 23 '16 at 17:15
• @robertbristow-johnson Yes it is... the "how about" part. I have my guesses there though. – Olli Niemitalo Oct 23 '16 at 17:17
• BTW, the gaussian is AN eigenfunction of the Fourier Transform. there are actually an infinite number of eigenfunctions and it's pretty easy to construct an eigenfunction. – robert bristow-johnson Oct 23 '16 at 17:55
• @robertbristow-johnson Getting somewhat tangential, but what else is? – Olli Niemitalo Oct 23 '16 at 20:11
• any even-symmetry function plus the fourier transform of itself (with $f$ replaced with $t$) is an eigenfunction of the Fourier Transform. – robert bristow-johnson Oct 24 '16 at 7:06

Interpolated points of the DFT can be computed using a dot-product of a few samples around the peak region with a pre-computed interpolation vector. The interpolation vector is determined by the location of the desired interpolated sample, taking into consideration the amount of zero-padding required, etc., etc.

This technique and the methodology to compute the interpolating vectors is covered in Appendix B of this document:

http://ericjacobsen.org/FTinterp.pdf

I hope that helps a bit.

• Thanks Eric, it is nice to see that you also have noted on magnitude interpolation that "The non-linearity in the magnitude detection results in a signal which is essentially undersampled". I show in my answer that squared magnitude is not undersampled when the time domain signal has been zero padded to 2x length. Unfortunately not all equations in the pdf display correctly in Chrome, for example Eq. 3.5. – Olli Niemitalo Oct 23 '16 at 20:50
• It's not undersampled until the magnitude is detected, so interpolating the complex coefficients to get the complex interpolated coefficient and then computing magnitude (or squared magnitude) avoids that problem. – Eric Jacobsen Oct 23 '16 at 21:10
• Agreed! The pdf works fine in Edge. – Olli Niemitalo Oct 23 '16 at 21:45

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.

P.S. Just found out there is also interpft, which does the interpolation with less syntax.

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}$$