Issue 1:
signal.resample is broken for even-length inputs.

Reproducing code example:

>>> from scipy.signal import resample
>>> resample(resample([1, 2, 3, 4], 6), 4)
array([ 1.25, 1.75, 3.25, 3.75])

I would expect resampling to a higher number of samples and back to be the identity operation as a consequence of Fourier resampling being "perfect" for band-limited signals. This appears to be due to incorrect handling of the Nyquist frequency component.

Scipy/Numpy/Python version information:

>>> import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
('0.14.1', '1.8.2', sys.version_info(major=2, minor=7, micro=10, releaselevel='final', serial=0))

Issue 2:
Originally reported on SO here.

It seems that scipy.signal.resample() makes errors when downsampling to an even number of points. For example, if we upsample a function to a multiple of the original points and then downsample again, we should get the original function back.

from scipy import signal
import numpy as np

def test_resample(n1,n2): # upsample from n1 to n2 points and back
x1=np.arange(n1)
y1=np.sin(x1)
y2,x2=signal.resample(y1,n2,x1)
y3,x3=signal.resample(y2,n1,x2)
print np.allclose(y1,y3)

But this fails when the lower number of points is even:

test_resample(10,20)
False

test_resample(11,22)
True

test_resample(11,33)
True

The problem occurs at the downsampling step. The errors are large, at least several percent for functions I tested.

I found the bug in the code for resample(), see these lines:

sl[axis] = slice(0, (N + 1) // 2)
Y[sl] = X[sl]
sl[axis] = slice(-(N - 1) // 2, None)
Y[sl] = X[sl]
y = fftpack.ifft(Y, axis=axis) * (float(num) / float(Nx))

These lines copy the FFT of the input function X to the FFT of the output function Y. Suppose we are downsampling and the number of points in the output spectrum, N, is even. Then Y[N/2] ends up containing the component of X at frequency -N/2. However, the subsequent IFFT expects that this should hold the sum of the positive and negative frequency components, namely X[N/2]+X[-N/2]. This can be solved by adding a few lines for the special condition like this:

sl[axis] = slice(0, (N + 1) // 2)
Y[sl] = X[sl]
sl[axis] = slice(-(N - 1) // 2, None)
Y[sl] = X[sl]
if N%2==0: # low number of points is even. So far we have set Y[-N/2]=X[-N/2]
if N<Nx: # if downsampling
sl[axis] = slice(N//2,N//2+1,None) # select the component at frequency N/2
Y[sl]+=X[sl] # add the component of X at N/2
elif N<num: # if upsampling
sl[axis] = slice(num-N//2,num-N//2+1,None) # select the component at frequency -N/2
Y[sl]/=2 # halve the component at -N/2
temp=Y[sl]
sl[axis] = slice(N//2,N//2+1,None) # select the component at +N/2
Y[sl]=temp # set that equal to the component at -N/2
y = fftpack.ifft(Y, axis=axis) * (float(num) / float(Nx))

So far this seems to work fine for my private version of resample(). Some professional should fix the scipy code.

Update 4/9/17:

signal.resample() fails also on upsampling, for the same reason. However, the bug does not affect the real part of the output for real-valued input (which is why I missed it at first). By inspecting the imaginary parts the bug becomes apparent. For example:

a=np.array([1.+0.j,0.+0.j])
print signal.resample(a,4)

leads to

[ 1.0+0.j 0.5-0.5j 0.0+0.j 0.5+0.5j]

This is wrong, because a has no imaginary part, and therefore its upsampled version should not either.

Again, this behavior is fixed by the code correction suggested above.