Upsample data using FFTs. How is this exactly done?

Question

I am basing my question on this post here, because I would like additional details on it, as I have not had any success in re-creating it.

I would like to simply learn, in no uncertain terms, how to do upsampling via the FFT. There seem to be a lot of details regarding if the signal is even or odd, where exactly to zero-pad, how much exactly, etc, and I am not getting those correct numbers for some reason.

Question: So let us say I have an even or odd length signal, of length $N$. I would like to upsample it to, say, $2N$, using FFTs.

How is this done exactly?

What I tried:

I tried following the post and others on the web but there seems to be very little consistency. It basically amounts to "zero-padding the frequency domain and then IFFTing" but I do not get proper results. In my case, I have a signal of length $100$, and I would like to upsample it to length $200$. So I do the following command:

ifft(fftshift([zeros(1,50) fftshift(fft(signal, 100)) zeros(1,50)]));

This however gives me complex data back, and not a signal that is simply upsampled by a factor of $2$. I am not sure where the problem is.

Thank you.

In short, you have been misinterpreting how the spectrum is stored. Thanks for the answer, Peter K. — user7358, Mar 10 '14 at 19:37

Peter K. · Accepted Answer · 2014-03-10T18:48:34.827

One trick, for even-length signals, is what to do with the "middle" sample. Here, I've split it half and half between each side of the FFT.

The other trick is to ensure that you have the right amplitudes in the resampled signal. Here's it's a factor of 2.

Try this in scilab:

x = rand(1,100,'normal');
X = fft(x);
XX = 2*[X(1:50) X(51)/2 zeros(1,99) X(51)/2 X(52:100)];
xx = ifft(XX);

clf;
plot([0:199]/200,xx)
plot([0:199]/200,xx,'go')
plot([0:99]/100,x,'r.')

which appears to generate the right thing. The red dots are the original samples, and the green circles (and blue connecting lines) are the resampled data.

enter image description here

UPDATE

For the case of an odd number of samples, there are fewer fiddly bits (no splitting the last coefficient):

x2 = rand(1,101,'normal');
X2 = fft(x2);
XX2 = 2*[X2(1:51) zeros(1,101)  X2(52:101)];
xx2 = ifft(XX2);

clf;
plot([0:201]/202,xx2)
plot([0:201]/202,xx2,'gx')
plot([0:100]/101,x2,'r.')

Thanks Peter. May you also please place an addition to your answer for the odd-cases for completion, based on the same example? Much appreciated! — TheGrapeBeyond, Mar 10 '14 at 14:02

score -1 · Answer 2 · answered Jul 11 '22 at 04:23

To extend PeterK's answer and refer a derivation, having a Nyquist bin = aliasing! This means no general-purpose recovery (or upsampling) exists, but PeterK's approach for real signals at worst loses an imaginary part of one RFFT bin, which is perhaps the best we can get.

This has an overlookable implication: in some applications, we keep all FFTs as powers of 2 for speed, and do filtering/conv/etc with pow2. The processing pipeline may have subsampling and upsampling steps, or, we may wish to test the quality of our lowpass filtering by comparing recovery against original. What we may miss is, if the downsampled signal has a non-zero Nyquist, we've already aliased. This is the only bin with such property.

Python version, with example in the linked answer:

xf, N = fft(x), len(x)
xupf = 2 * hstack([xf[:N//2], 
                   xf[N//2]/2, zeros(N - 1), xf[N//2]/2, 
                   xf[-(N//2 - 1):]])
xup = ifft(xupf)

Upsample data using FFTs. How is this exactly done?

2 Answers2

Linked