从可变速率重新采样到固定速率的算法是什么？

#1 楼

最简单的方法是进行某种样条插值，如吉姆·克莱（Jim Clay）建议的（线性或其他方式）。但是，如果您有大量的批处理功能，尤其是如果您有一组不确定的不均匀样本，那么有一种“完美重构”算法非常优雅。由于数量上的原因，它可能并非在所有情况下都可行，但是至少值得在概念上有所了解。我首先在本文中读过它。

诀窍在于将您的一组非均匀样本视为已经通过正弦插值从均匀样本中重构出来。遵循本文中的注释：

$$
y（t）= \ sum_ {k = 1} ^ {N} {y（kT）\ frac {\ sin（\ pi （t-kT）/ T）} {\ pi（t-kT）/ T}} = \ sum_ {k = 1} ^ {N} {y（kT）\ mathrm {sinc}（\ frac {t-kT } {T}）}。
$$

请注意，这提供了一组线性方程，每个非均匀样本$ y（t）$都具有一个线性方程，其中未知数等于间隔的样本$ y（kT）$，例如：

$$
\ begin {bmatrix} y（t_0）\\ y（t_1）\\ \ cdots \\ y（t_m ）\ end {bmatrix} = \ begin {bmatrix} \ mathrm {sinc}（\ frac {t_0-T} {T}）＆\ mathrm {sinc}（\ frac {t_0-2T} {T}）＆\ cdots ＆\ mathrm {sinc}（\ frac {t_0-nT} {T}）\\ \ mathrm {sinc}（\ frac {t_1-T} {T}）＆\ mathrm {sinc}（\ frac {t_1-2T } {T}）＆\ cdots＆\ mathrm {sinc}（\ frac {t_1-nT} {T}）\\ \ cdots＆\ cdots＆\ cdots＆\ cdots \\ \ mathrm {sinc}（\ frac { t_m-T} {T}）和\ mathrm {sinc}（\ frac {t_m-2T} {T}）＆\ cdots和\ mathrm {sinc}（\ frac {t_m-nT} {T}）\ end { bmatrix} \ begin {bmatrix} y（T）\\ y（2T）\\ \ cdots \\ y（nT）\ end {bmatrix}。
$$

在上式中，$ n $是未知均匀样本的数量，$ T $是均匀样本速率的倒数，$ m $是非均匀样本的数量（可能大于$ n $）。通过计算该系统的最小二乘解，可以重建均匀样本。从技术上讲，仅需要$ n $个非均匀样本，但是根据它们在时间上的“分散”程度，插值矩阵可能会变得异常不适。在这种情况下，使用更多不均匀的样本通常会有所帮助。

作为一个玩具示例，这是上述方法与轻度抖动网格上三次样条插值的比较​​（使用numpy）：


（此答案的末尾包含用于复制上述图的代码）

对于高质量，鲁棒的方法，从以下论文之一开始，可能会更合适：


A. Aldroubi和Karlheinz Grochenig，非不变采样和不变不变空间中的重构，SIAM 。，2001，no。 4，

K. Grochenig和H. Schwab，《不变移位空间中非均匀采样的快速局部重建方法》，SIAM J.矩阵分析。
应用。 ，24（2003），899-
913.


-

import numpy as np
import pylab as py

import scipy.interpolate as spi
import numpy.random as npr
import numpy.linalg as npl

npr.seed(0)

class Signal(object):

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y ,'bo-')
        py.ylim([-1.8,1.8])
        py.plot(hires.x,hires.y, 'k-', alpha=.5)

    def _plot(self, title):
        py.grid()
        py.title(title)
        py.xlim([0.0,1.0])

    def sinc_resample(self, xnew):
        m,n = (len(self.x), len(xnew))
        T = 1./n
        A = np.zeros((m,n))

        for i in range(0,m):
            A[i,:] = np.sinc((self.x[i] - xnew)/T)

        return Signal(xnew, npl.lstsq(A,self.y)[0])

    def spline_resample(self, xnew):
        s = spi.splrep(self.x, self.y)
        return Signal(xnew, spi.splev(xnew, s))

class Error(Signal):

    def __init__(self, a, b):
        self.x = a.x
        self.y = np.abs(a.y - b.y)

    def plot(self, title):
        self._plot(title)
        py.plot(self.x, self.y, 'bo-')
        py.ylim([0.0,.5])

def grid(n): return np.linspace(0.0,1.0,n)
def sample(f, x): return Signal(x, f(x))

def random_offsets(n, amt=.5):
    return (amt/n) * (npr.random(n) - .5)

def jittered_grid(n, amt=.5):
    return np.sort(grid(n) + random_offsets(n,amt))

def f(x):
    t = np.pi * 2.0 * x
    return np.sin(t) + .5 * np.sin(14.0*t)

n = 30
m = n + 1

# Signals
even   = sample(f, np.r_[1:n+1] / float(n))
uneven = sample(f, jittered_grid(m))
hires  = sample(f, grid(10*n))

sinc   = uneven.sinc_resample(even.x)
spline = uneven.spline_resample(even.x)

sinc_err   = Error(sinc, even)
spline_err = Error(spline, even)

# Plot Labels
sn = lambda x,n: "%sly Sampled (%s points)" % (x,n)
r  = lambda x: "%s Reconstruction" % x
re = lambda x: "%s Error" % r(x)

plots = [
    [even,       sn("Even", n)],
    [uneven,     sn("Uneven", m)],
    [sinc,       r("Sinc")],
    [sinc_err,   re("Sinc")],
    [spline,     r("Cubic Spline")],
    [spline_err, re("Cubic Spline")]
]

for i in range(0,len(plots)):
    py.subplot(3, 2, i+1)
    p = plots[i]
    p[0].plot(p[1])

py.show()

#2 楼

这听起来像是异步采样率转换的问题。为了从一种采样率转换为另一种采样率，我们可以通过执行Sinc插值来计算信号的连续时间表示，然后以新的采样率进行重新采样。你在做什么并没有太大的不同。您需要重新采样信号以使采样时间固定。

连续时间信号可以通过将每个采样与Sinc函数卷积来计算。由于sinc函数可以无限期地继续，因此我们使用更实用的方法，例如具有实际有限长度的窗口sinc。棘手的部分是，因为您的样本会随时间移动，所以在重采样时可能需要对每个样本使用具有不同相位偏移的正弦信号。

来自采样信号的连续时间信号：

$ x（t）= \ sum \ limits_ {n =-\ infty} ^ \ infty x [n] sinc（{t-nT_s \ over T_s}）$

其中$ T_s $是您的采样时间。但是，根据您的情况，采样时间不是固定的。因此，我认为您需要用该示例的采样时间替换它。

$ x（t）= \ sum \ limits_ {n =-\ infty} ^ \ infty x [n] sinc（ {t-nT_s [n] \ over T_s [n]}）$

您可以从中重新采样信号：

$ y [n] = x（nT_ { ns} $）

其中$ T_ {ns} $是所需的采样时间。

将它们全部放在一起即可得到：

$ y [m] = \ sum \ limits_ {n =-\ infty} ^ \ infty x [n] sinc（{mT_ {ns} -nT_s [n] \ over T_s [n]}）$

因为这不是因果关系或难于处理的，所以可以用有限支持函数代替sinc函数，并相应地调整总和的限制。长度为2k的则为：

$ y [m] = \ sum \ limits_ {n = -k} ^ kx [n] kernel（{mT_ {ns} -nT_s [n] \ over T_s [ n]}）$

我希望这会有所帮助...，但我可能在此过程中犯了一个错误，并且可能需要大量数学运算。我建议研究采样率转换以获取更多信息。也许这里的其他人也可以给出更好的解释或解决方案。

#3 楼

我认为Jacob的答案非常可行。我将使用线性插值（简单，信号性能不佳）或三次样条曲线（仍然不太硬，信号性能更好）在任意时间采样任意时间生成采样。

#4 楼

（一个月后），任何插值方法都有两个主要选择：
1）最接近要使用的缺失点的数据点数，$ Nnear $ 2 4 6 ...
2）类使用的基本函数：线性，多项式，正弦余弦（傅立叶），分段三次（B样条或插值样条），类似Sinc的...
将直线拟合到$ Nnear $点很容易：
2点[-1，$ y_ { -1} $]，[1，$ y_1 $]：
$ \ qquad $估计$ y_0 \ sim（y _ {-1} + y_1）/ 2 $
点$ [x_i，y_i] $平均$ x_i = 0 $：
$ \ qquad y_0 \ sim $平均$ y_i $
一般$ [x_i，y_i] $：
$ \ qquad $参见例如>食谱p。 781：
拟合一条线$ a + bx $
并估计$ y_0 \ sim a $。
$ \ qquad $一个可以拟合二次方，三次方，正弦余弦...

我理解您的数据均匀分布，但缺少一些点，
对吗？
这种情况下线性插值的效果如何？
/>那么，让我们尝试使用$ f $ = 0.25的cos $ 2 \ pi ft $：1 0 -1 0 1 0 -1 0 ...
2个任意点的邻居平均为0，太糟糕了。 /> 4个邻居：[1 0（缺失-1）0 1]的平均值= 1/2，太糟了。
（尝试使用4个邻居过滤器[-1 3 3 -1] /4。）

与4或6或8个邻居进行线性插值可能对您的数据足够好
。
我建议您从一种您已经理解的方法入手
跳入样条线，类似于Sinc ...虽然它们也很有趣。



另一种完全不同的方法是反距离加权。
易于实现
（请参见idw-interpolation-with-python on SO），工作在2d和3d以上也是如此，
但是在理论上很难分析。

（显然，没有一种插值方法可能适合
现实中出现的成千上万的
[信号，噪声，误差度量，测试函数]的组合。
世界上有比测试功能更多的方法，带有更多的旋钮。
尽管如此，方法和测试功能的库还是有用的。）

#5 楼

如果您使用matlab，则可以通过时间序列来实现。

time  % is your starting vector of time

data % vector of data you want to resample 

data_TS = timeseries(data,time); % define the data as a timeseries 

new_time = time(0):dt:time(end); % new vector of time with fixed dt=1/fs

data_res = resample(data_TS,new_time); % data resampled at constant fs

#6 楼

在进行一些特殊处理之前，您可以尝试类似以下的简单操作（伪代码-无插值，但可以添加）

TimeStamp[]  //Array of Sensor TimeStamps -NULL terminated – TimeStamp[i] corresponds to Reading[i]
Reading[]      //Array of Sensor Readings       -NULL terminated

AlgorithmOut   //Delimited file of of readings in fixed sample time (5ms) 
CurrentSavedReading = Reading[0]

SampleTime=TimeStamp[0] //ms virtual sample time, 5ms fixed samples

i = 0 // loop index
While(TimeStamp[i] != NULL)
{
   FileWrite (CurrentSavedReading, AlgorithmOut)//write value to file
   SampleTime = SampleTime + 5//ms
   if(SampleTime > TimeStamp[i])
   {
      i++
      CurrentSavedReading = Reading[i]
   }
}

#7 楼

恕我直言，Datageist的答案是正确的，Jacob的答案是不正确的。一种简单的验证方法是，保证datageist的建议算法可以对原始样本进行插值（假设数值精度无限大），而Jacob的答案则不能。 Sinc函数的集合是正交的：如果每个移位的Sinc函数在输入样本上离散，则它们将形成一个无限的单位矩阵。这是因为对于除n = 0外的所有n，sin（n pi）/（n pi）均为零。在输入样本上，将产生一个非平凡的矩阵。因此，周围样本的贡献将不会为零，并且重构将不再通过输入样本进行插值。