<p dir="auto">signal.resample is broken for even-length inputs.</p>
<h3 dir="auto">Reproducing code example:</h3>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; from scipy.signal import resample
&gt;&gt;&gt; resample(resample([1, 2, 3, 4], 6), 4)
array([ 1.25,  1.75,  3.25,  3.75])"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; from scipy.signal import resample
&gt;&gt;&gt; resample(resample([1, 2, 3, 4], 6), 4)
array([ 1.25,  1.75,  3.25,  3.75])
</code></pre></div>
<p dir="auto">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.</p>
<h3 dir="auto">Scipy/Numpy/Python version information:</h3>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; 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))"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; 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))
</code></pre></div>
<p dir="auto">Originally reported on SO <a href="http://stackoverflow.com/questions/43283845/scipy-signal-resample-malfunction-with-even-number-of-points" rel="nofollow">here</a>.</p>
<p dir="auto">It seems that <strong>scipy.signal.resample()</strong> 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.</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="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)"><pre class="notranslate"><code class="notranslate">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)
</code></pre></div>
<p dir="auto">But this fails when the lower number of points is even:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="test_resample(10,20)
    False

test_resample(11,22)
    True

test_resample(11,33)
    True"><pre class="notranslate"><code class="notranslate">test_resample(10,20)
    False

test_resample(11,22)
    True

test_resample(11,33)
    True
</code></pre></div>
<p dir="auto">The problem occurs at the downsampling step. The errors are large, at least several percent for functions I tested.</p>
<p dir="auto">I found the bug in the code for <strong>resample()</strong>, see <a href="https://github.com/scipy/scipy/blob/v0.19.0/scipy/signal/signaltools.py#L2231-L2235">these lines</a>:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="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))"><pre class="notranslate"><code class="notranslate">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))
</code></pre></div>
<p dir="auto">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:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="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&lt;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&lt;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))"><pre class="notranslate"><code class="notranslate">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&lt;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&lt;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))
</code></pre></div>
<p dir="auto">So far this seems to work fine for my private version of resample(). Some professional should fix the scipy code.</p>
<hr>
<p dir="auto">Update 4/9/17:<br>
<strong>signal.resample()</strong> 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:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="a=np.array([1.+0.j,0.+0.j])
print signal.resample(a,4)"><pre class="notranslate"><code class="notranslate">a=np.array([1.+0.j,0.+0.j])
print signal.resample(a,4)
</code></pre></div>
<p dir="auto">leads to</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="[ 1.0+0.j   0.5-0.5j  0.0+0.j   0.5+0.5j]"><pre class="notranslate"><code class="notranslate">[ 1.0+0.j   0.5-0.5j  0.0+0.j   0.5+0.5j]
</code></pre></div>
<p dir="auto">This is wrong, because <strong>a</strong> has no imaginary part, and therefore its upsampled version should not either.<br>
Again, this behavior is fixed by the code correction suggested above.</p>