scipy.信号：筛选可变时间数据

Question

我有一个[t, y(t)]形式的大型数据集，我想使用scipy.signal（特别是scipy.filter.butter和scipy.filter.filtfilt）对其应用IIR低通滤波器（一阶或二阶巴特沃斯就足够了）。问题是t不是规则间隔的，这似乎是scipy.signal中函数的一个要求。你知道吗

对于任何“缺失”的点，我知道我的信号与其之前的值保持不变（因此给定我的t-数据中的两个连续点t1和t2，以及数据中的一个点T，这样t1<T<t2，我正在采样的“实”函数Y(t)将取值Y(T)=Y(t1)）。t是整数值，因此我可以简单地添加缺失的点，但这会导致数据集的大小增加约10倍，这是有问题的，因为它已经非常大了。你知道吗

所以问题是，有没有一种（足够简单和低开销）的方法来过滤我的数据集而不插入所有的“缺失”点？你知道吗

Answer 1

您可以有效地将数据“包装”到函数中。你知道吗

如果您的数据是列表形式的，那么您需要将其转换为dict，并创建t值的排序列表。然后，您可以使用^{}模块中的列表对分算法对缺少的值进行插值。你知道吗

这里有一些用python2编写的演示代码，但是如果需要的话，可以直接将其转换为python3。你知道吗

from random import seed, sample
from bisect import bisect

#Create some fake data
seed(37)
data = dict((u, u/10.) for u in sample(xrange(50), 25))

keys = data.keys()
keys.sort()
print keys

def interp(t):
    i = bisect(keys, t) 
    k = keys[max(0, i-1)]
    return data[k]

for i in xrange(50):
    print i, interp(i)    

输出

[2, 4, 8, 10, 14, 15, 19, 21, 22, 23, 26, 27, 29, 30,
 32, 33, 34, 35, 37, 38, 39, 42, 43, 44, 48]
0 0.2
1 0.2
2 0.2
3 0.2
4 0.4
5 0.4
6 0.4
7 0.4
8 0.8
9 0.8
10 1.0
11 1.0
12 1.0
13 1.0
14 1.4
15 1.5
16 1.5
17 1.5
18 1.5
19 1.9
20 1.9
21 2.1
22 2.2
23 2.3
24 2.3
25 2.3
26 2.6
27 2.7
28 2.7
29 2.9
30 3.0
31 3.0
32 3.2
33 3.3
34 3.4
35 3.5
36 3.5
37 3.7
38 3.8
39 3.9
40 3.9
41 3.9
42 4.2
43 4.3
44 4.4
45 4.4
46 4.4
47 4.4
48 4.8
49 4.8

（我手动包装了keys的输出，以便在没有水平滚动的情况下更易于阅读）。你知道吗

通过在一行中重新写入插值函数体，您将获得一个微小的加速：

def interp(t):
    return data[keys[max(0, bisect(keys, t)-1)]]

它的可读性要差得多，IMHO，但是如果函数被频繁调用，那么速度差可能是值得的。你知道吗

Answer 2

pm2ring给出的答案是可行的，但是假设您的数据已经由t排序，那么它的效率就比可能的低。它需要对数线性时间和线性附加空间。您可以编写一个生成器，以在线性时间和恒定的附加空间中生成具有规则采样间隔的转换数据集：

# Assumes that dataset rows are lists as described in the question:
# [[t1, Y(t1)], [t2, Y(t2)], [t3, Y(t3)], ..., [tz, Y(tz)]]
# If this assumption is wrong, just extract t and Y(t) in another way.

# The generated range starts at t1 and ends directly after tz.
# Warning: will overgenerate points if the data are more densely sampled
# than the requested sampling interval.

def step_interpolate(dataset, interval):
    left = next(dataset)   # [t1, Y(t1)]
    right = next(dataset)  # [t2, Y(t2)]
    t_regular = left[0]
    while True:
        if left is right:          # same list object
            right = next(dataset)  # iteration stops when dataset stops
        if right[0] <= t_regular:
            left = right
        yield [t_regular, left[1]]
        t_regular += interval

测试：

data = [[1, 10], [15, 2], [50, 100], [55, 17]]
for item in step_interpolate(iter(data), 10):
    print item[0], item[1]

输出：

1 10
11 10
21 2
31 2
41 2
51 100
61 17