Numpy/Scipy irfft strange behaviour

Question

I am resampling a real signal, and since I have at my disposal its fft from rfft, I want to use irfft(signal, new_length). But I can't seem to get it working.

This is a working code snippet that resamples a signal of length 4 using complex fft:

from numpy.fft import fft,ifft
p=array([1.,2.2,4.,1.])
pk=fft(p)
pnew=ifft(pk,8)*(8./4.)

where the factor (8./4.) rescales from the original to the new length. You can check that pnew[::2]==p.

Now, when I try to apply the same strategy with real Fourier transform, I get the wrong result at the original points:

from numpy.fft import rfft,irfft
p=array([1.,2.2,4.,1.])
pk=rfft(p)
pnew=irfft(pk,8)*(8./4.)

and I have pnew[::2]=[ 1.45, 1.75, 4.45, 0.55]!=p.

Does anybody have a clue of what is going on? I have tried using the routines from scipy, with the same result. The documentation itself discusses briefly how to do this, see here, bottom of the page

Answer 1

The documentation you like to says:

In other words, irfft(rfft(a), len(a)) == a to within numerical accuracy.

This is not the case if you do irfft(pk, 8)! The problem is due the odd samples and the symmetry of the Fourier transform in addition to your padding. Note that there are no problems at all all if len(p) is odd.

For a better understanding consider this:

>>> p = np.array([1.,2.2,4.,1.])
>>> np.fft.fft(p)
array([ 8.2+0.j , -3.0-1.2j,  1.8+0.j , -3.0+1.2j])
>>> np.fft.fftfreq(len(p))
array([ 0.  ,  0.25, -0.5 , -0.25]) # 0.5 only occurs once negative
>>> np.fft.rfft(p)
array([ 8.2+0.j , -3.0-1.2j,  1.8+0.j ])
>>> np.fft.rfftfreq(len(p)) # (not available in numpy 1.6.)
array([ 0.  ,  0.25,  0.5 ]) # 0.5 occurs, here positive, it does not matter

# also consider the odd length FFT
>>> np.fft.fftfreq(len(p)+1)
array([ 0. ,  0.2,  0.4, -0.4, -0.2]) # 0.4 is in there twice.

# And consider that this gives the result you expect:
>>> symmetric_p = np.fft.rfft(p)
>>> symmetric_p[-1] /= 2
>>> np.fft.irfft(symmetric_p, 8)[::2]*(8./4.)
array([ 1. ,  2.2,  4. ,  1. ])

Which means if you look closely. The FFT frequencies calculated are not symmetric if the input samples are even, instead there is an extra negative frequency (which actually could just as well be a positive frequency, since it always has no phase shift).

Because you are padding (for no real reason?) to a different frequency, the RFFT suddenly has extra "room" for this frequency. So if you look at it from the FFT point of view, you add this normally only once occuring negative frequency also as a positive frequency (which basically means it goes in double). If you look above the symmetric_p halving this frequency gives the expected result with padding (it will not give the expected result without padding).