Unexpected results for fourier resampling (Python 3)

Question

I am trying to resample a signal using the fourier method in python using np.fft.rfft() and np.fft.irfft(). (I'm not using scipy.signal.resample() because in my final application I need to do more to my signal in the frequency domain than what can be accomplished by the window argument of scipy.signal.resample().) I have reduced the problems I see to this simple test case (ipython notebook):

%pylab inline

N = 16
P = 8

S = [0.49, 1.68, 0.78, 0.05, 0.21, 0.29, 1.34, 1.17, 16.73,
     48.78, 16.90, 0.62, 0.40, 1.60, 0.86, 0.57]
assert len(S) == N

figure(figsize=(14, 5))
plot(np.arange(N), S, '.')

dft = np.fft.rfft(S)
plot(np.arange(N), np.fft.irfft(dft), '-')

dft *= P
dft.resize(P*N/2+1)
plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')

xlim(6.5, 11.5)
ylim(-5, 55);

And here is what I get back:

As you can see, the result looks almost correct, but is nevertheless a bit off. (The resampled red curve does not go through the original sample points, even though there are sample points at the same positions in time as in the original signal.)

The greeen line shows that everything is still ok before extending and scaling the dft array. So the problem seems to be in the actual resampling:

dft *= P
dft.resize(P*N/2+1)
plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')

I don't see how this result is mathematically possible because other than scaling of time and amplitude axis, adding zeroes to the high-frequency end of the dft array should only add zero terms to the fourier series expansion of the signal. So the signal should still pass through the original sample values.

But the small error I see does not seem to be related to scaling: sometimes the resampled waveform lags behind the original (t=8) and sometimes it is faster (t=10), sometimes it is larger in magnitude (t=7,9,11) and sometimes it is smaller (t=8,10).

I'm playing with this for hours now and can't find what is wrong with my code.. Maybe I just don't see it because I'm already staring at it for so long now..

(see np.fft documentation for the definitions used by the python FFT implementations)

Edit (in response to comment): When I replace S with an impulse at t=9, or a fast impulse train, I get the following:

Answer 1

If you perform the full DFT of a signal with odd number of samples, the frequency values in the spectrum are a DC coefficient and (N-1)/2 conjugate pairs of sinusoid coefficients.

If you perform the full DFT of a signal with even number of samples, the frequency values in the spectrum are a DC coefficient, (N-1)/2 conjugate pairs of sinusoid coefficients, and a single coefficient corresponding to frequency 0.5.

rfft only keeps the positive frequency part of the spectrum, and the inverse, irfft, presumably recreates the negative part of the spectrum before performing an inverse-DFT.

For example, for an odd size signal with rfft of [0 1 2 3], the full spectrum would be [3 2 1 0 1 2 3].

For an even size signal with rfft of [0 1 2 3 4] it presumably has a flag set somewhere in the code that indicates an even sized signal, and so the full spectrum would be [3 2 1 0 1 2 3 4]. Note the 4 only appears once.

When you pad the rfft result though, e.g. [0 1 2 3 4 0 0 0 0], it probably recreates the spectrum as [0 0 0 4 3 2 1 0 1 2 3 4 0 0 0 0]. That is, the 4 is now repeated twice, and therefore that frequency component is twice as strong.

Therefore, before padding you must halve the last rfft coefficient iff the signal has an even number of samples.