Resampling a signal with scipy.signal.resample

Question

I was trying to resample a generated signal from 256 samples to 20 samples using this code:

import scipy.signal 
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 10, 256, endpoint=False)
y = np.cos(-x**2/6.0)
yre = signal.resample(y,20)
xre = np.linspace(0, 10, len(yre), endpoint=False)
plt.plot(x,y,'b', xre,yre,'or-')
plt.show()

Which returns this plot (apparently correct):

However, as can be noticed, the first sample is badly approximated. I believe that resample computes the average of the samples that belongs to equidistant groups of samples and, in this case, it seems that the first subgroup of samples is padded with zeros in the beggining in order to estimate the first output sample.

Thus, I consider that the first sample can be successfully estimated by telling resample function that I do not want to pad with zeroes the first subgroup.

Can somebody help me in order to achieve a correct resampling of this signal?

Thanks in advance.

Answer 1

As the reference page for scipy.signal.resample states, it uses FFT methods to perform the resampling.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html

One of the side effects is the implicit assumption (because of the underlying FFT) that the signal is periodic; hence if there is a large step from x[0] to x[-1], the resample will struggle to make them meet: the FFT thinks that the time-like axis is not a line, but a circle.

The FFT is a powerful tool, but it's a powerful tool with sharp edges that can cut you.

Answer 2

I had similar problem. Found solution on the net that seems to be also faster than scipy.signal.resample (https://github.com/nwhitehead/swmixer/blob/master/swmixer.py). It is based on np.interp function. Added also scipy.signal.resample_poly for comparison (which is not very best in this case).

import scipy.signal 
import matplotlib.pyplot as plt
import numpy as np

# DISCLAIMER: This function is copied from https://github.com/nwhitehead/swmixer/blob/master/swmixer.py, 
#             which was released under LGPL. 
def resample_by_interpolation(signal, input_fs, output_fs):

    scale = output_fs / input_fs
    # calculate new length of sample
    n = round(len(signal) * scale)

    # use linear interpolation
    # endpoint keyword means than linspace doesn't go all the way to 1.0
    # If it did, there are some off-by-one errors
    # e.g. scale=2.0, [1,2,3] should go to [1,1.5,2,2.5,3,3]
    # but with endpoint=True, we get [1,1.4,1.8,2.2,2.6,3]
    # Both are OK, but since resampling will often involve
    # exact ratios (i.e. for 44100 to 22050 or vice versa)
    # using endpoint=False gets less noise in the resampled sound
    resampled_signal = np.interp(
        np.linspace(0.0, 1.0, n, endpoint=False),  # where to interpret
        np.linspace(0.0, 1.0, len(signal), endpoint=False),  # known positions
        signal,  # known data points
    )
    return resampled_signal

x = np.linspace(0, 10, 256, endpoint=False)
y = np.cos(-x**2/6.0)
yre = scipy.signal.resample(y,20)
xre = np.linspace(0, 10, len(yre), endpoint=False)

yre_polyphase = scipy.signal.resample_poly(y, 20, 256)
yre_interpolation = resample_by_interpolation(y, 256, 20)

plt.figure(figsize=(10, 6))
plt.plot(x,y,'b', xre,yre,'or-')
plt.plot(xre, yre_polyphase, 'og-')
plt.plot(xre, yre_interpolation, 'ok-')
plt.legend(['original signal', 'scipy.signal.resample', 'scipy.signal.resample_poly', 'interpolation method'], loc='lower left')
plt.show()

CARE! This method, however, seems to perform some unwanted low-pass filtering.

x = np.linspace(0, 10, 16, endpoint=False)
y = np.random.RandomState(seed=1).rand(len(x))
yre = scipy.signal.resample(y, 18)
xre = np.linspace(0, 10, len(yre), endpoint=False)

yre_polyphase = scipy.signal.resample_poly(y, 18, 16)
yre_interpolation = resample_by_interpolation(y, 16, 18)

plt.figure(figsize=(10, 6))
plt.plot(x,y,'b', xre,yre,'or-')
plt.plot(xre, yre_polyphase, 'og-')
plt.plot(xre, yre_interpolation, 'ok-')
plt.legend(['original signal', 'scipy.signal.resample', 'scipy.signal.resample_poly', 'interpolation method'], loc='lower left')
plt.show()

Still, this is the best result I got, but I hope someone will provide something better.