使用 scipy.signal.resample 重新采样信号答案

【问题标题】：Resampling a signal with scipy.signal.resample使用 scipy.signal.resample 重新采样信号
【发布时间】：2018-12-27 11:45:15
【问题描述】：

我试图使用此代码将生成的信号从 256 样本重新采样为 20 样本：

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()

返回这个情节（显然是正确的）：

但是，可以注意到，第一个样本的近似值很差。我相信resample 计算属于等距样本组的样本的平均值，在这种情况下，似乎第一个样本子组在开始时用零填充，以估计第一个输出样本。

因此，我认为可以通过告诉resample 函数我不想用零填充第一个子组来成功估计第一个样本。

有人可以帮我正确重采样这个信号吗？

提前致谢。

【问题讨论】：

标签： python python-3.x scipy

【解决方案1】：

我有类似的问题。在网上找到了似乎也比scipy.signal.resample (https://github.com/nwhitehead/swmixer/blob/master/swmixer.py) 更快的解决方案。它基于np.interp 函数。还添加了scipy.signal.resample_poly 进行比较（在这种情况下不是最好的）。

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()

小心！然而，这种方法似乎执行了一些不需要的低通滤波。

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()

尽管如此，这是我得到的最好结果，但我希望有人能提供更好的东西。

【讨论】：

np.interp 方法实际上不会对您的信号进行低通滤波，而是在数据点之间进行线性插值。如果您查看估计的数据点，它们通常正好位于原始点的线性连接处。在第二个示例中，您的数据点在高值和低值之间交替。因此，插值点自然会位于它们之间的某个位置，从而呈现出低通滤波信号的（错误）外观。
降低采样率将始终“低通滤波”您的信号。想想奈奎斯特定理：en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem
@N4ppeL 如果没有高于新奈奎斯特频率的频率，那么良好的重采样不应该“低通滤波器”任何东西

【解决方案2】：

作为 scipy.signal.resample 状态的参考页面，它使用 FFT 方法执行重采样。

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

其中一个副作用是信号是周期性的隐含假设（由于底层 FFT）；因此，如果从 x[0] 到 x[-1] 有很大的步长，重采样将难以使它们相遇：FFT 认为类时间轴不是一条线，而是一个圆。

FFT 是一种强大的工具，但它是一种具有锋利边缘的强大工具，可能会割伤你。

【讨论】：