【问题标题】:Resampling (upsampling, interpolating) a series of numbers重采样(上采样、插值)一系列数字
【发布时间】:2018-06-22 03:37:12
【问题描述】:

我有一个逗号分隔的整数值系列,我想重新采样,这样我就有两倍的整数值,其中一个新值被添加到每个现有值的中间。例如,如果这是我的来源:

1,5,11,9,13,21

结果是:

1,3,5,8,11,10,9,11,13,17,21

如果不清楚,我会尝试在源系列中的每个值之间添加一个数字,如下所示:

1   5   11    9    13    21
1 3 5 8 11 10 9 11 13 17 21

我已经搜索了很多,似乎 scipy.signal.resample 或 panda 之类的东西应该可以工作,但我对此完全陌生,我无法让它工作。例如,这是我对 scipy 的尝试之一:

import numpy as np
from scipy import signal
InputFileName = "sample.raw"
DATA250  = np.loadtxt(InputFileName, delimiter=',', dtype=int);
print(DATA250)
DATA500 = signal.resample(DATA250, 11)
print(DATA500)

哪些输出:

[ 1  5 11  9 13 21]
[ 1.         -0.28829461  6.12324489 10.43251996 10.9108191   9.84503237
  8.40293529 10.7641676  18.44182898 21.68506897 12.68267746]

显然,我错误地使用了 signal.resample。有没有办法用 signal.resample 或 panda 做到这一点?我应该使用其他方法吗?

此外,在我的示例中,所有源数字之间都有一个整数。在我的实际数据中,情况并非如此。因此,如果其中两个数字是 10,15,则新数字将是 12.5。但是我想让所有的结果数字都是整数。所以插入的新数字需要是 12 或 13(对我来说它是什么并不重要)。

请注意,一旦我完成这项工作,源文件实际上将是一个逗号分隔的 2,000 个数字列表,输出应该是 4,000 个数字(或者技术上说是 3,999,因为不会在末尾添加一个)。此外,这将用于处理类似于 ECG 记录的内容 - 目前 ECG 以 250 Hz 的频率采样 8 秒,然后将其传递给单独的进程以分析记录。然而,这个单独的过程需要以 500 Hz 的频率对记录进行采样。所以工作流程是我每 8 秒进行一次 250 Hz 的记录并将其上采样到 500 Hz,然后将结果输出传递给分析过程。

感谢您提供的任何指导。

【问题讨论】:

  • 您只是将这些内容读入然后将它们写出来,除了插值之外没有任何处理?如果是这样,首先不使用numpy 可能更简单。只需对值进行循环,记住最后一个值。然后,每次将if last is not None: 写入(last + value) // 2,然后写入value 并设置last = value。如果这对你来说更容易理解和编码,谁在乎它在理论上是否不那么优雅?
  • 谢谢,我更新了描述以表明我所做的只是读入字符串,然后获取输出并将其传递给不同的进程以进行进一步分析。

标签: python pandas numpy interpolation resampling


【解决方案1】:

由于插值简单,你可以手动完成:

import numpy as np
a = np.array([1,5,11,9,13,21])
b = np.zeros(2*len(a)-1, dtype=np.uint32)
b[0::2] = a
b[1::2] = (a[:-1] + a[1:]) // 2

你也可以这样使用scipy.signal.resample

import numpy as np
from scipy import signal
a = np.array([1,5,11,9,13,21])
b = signal.resample(a, len(a) * 2)
b_int = b.astype(int)

诀窍是元素的数量正好是两倍,这样奇数点就可以匹配你的初始点。另外我认为scipy.signal.resample 完成的傅立叶插值比您要求的线性插值更适合您的心电图信号。

【讨论】:

  • 谢谢,太好了。我想我将使用 scipy.signal.resample - 你是对的,它似乎确实产生了更平滑的曲线,尽管在每个图表的开头和结尾都有一些波动。无论如何,这会很好用。
  • scipy.signal.resample 假设信号是周期性的,因此它会尝试“匹配”数据的开头和结尾。看文档中的例子,很清楚:docs.scipy.org/doc/scipy/reference/generated/…
【解决方案2】:

虽然我可能会在这里只使用 NumPy,与 J. Martinot-Lagarde's answer 非常相似,但实际上你不必这样做。


首先,您可以仅使用 csv 模块读取一行逗号分隔的数字:

with open(path) as f:
    numbers = map(int, next(csv.reader(f))

…或者只是字符串操作:

with open(path) as f:
    numbers = map(int, next(f).split(','))

然后您可以轻松地进行插值:

def interpolate(numbers):
    last = None
    for number in numbers:
        if last is not None:
            yield (last+number)//2
        yield number
        last=number

如果您希望它完全通用且可重用,只需使用function 参数和yield function(last, number),并将None 替换为sentinel = object()


现在,您需要做的就是join 结果和write 他们:

with open(outpath, 'w') as f:
    f.write(','.join(map(str, interpolate(numbers))))

这个解决方案有什么优势吗?好吧,除了读取/拆分和加入/写入之外,它纯粹是懒惰的。我们可以很容易地编写惰性拆分和连接函数(或者只是手动编写)。因此,如果您不得不处理十亿个逗号分隔的数字而不是一千个,那么您只需要更改即可。

这是一个懒惰的split

def isplit(s, sep):
    start = 0
    while True:
        nextpos = s.find(sep, start)
        if nextpos == -1:
            yield s[start:]
            return
        yield s[start:nextpos]
        start=nextpos+1

您可以使用 mmap 作为延迟读取的字符串(嗯,bytes,但我们的数据是纯 ASCII,所以没关系):

with open(path, 'rb') as f:
    with mmap.mmap(inf.fileno(), 0, access=mmap.ACCESS_READ) as mm:
        numbers = map(int, isplit(mm, b','))

让我们为懒惰的写作使用不同的解决方案,只是为了多样化:

def icsvwrite(f, seq, sep=','):
    first = next(seq, None)
    if not first: return
    f.write(first)
    for value in seq:
        f.write(sep)
        f.write(value)

所以,把它们放在一起:

with open(inpath, 'rb') as inf, open(outpath, 'w') as outf:
    with mmap.mmap(inf.fileno(), 0, access=mmap.ACCESS_READ) as mm:
        numbers = map(int, isplit(mm, b','))
        icsvwrite(outf, map(str, interpolate(numbers)))

但是,即使我能够很快地将它们拼凑在一起,并且所有部件都可以很好地重复使用,我仍然可能会使用 NumPy 来解决您的特定问题。您不会阅读一行十亿个数字。您已经在唯一一台运行该脚本的机器上安装了 NumPy。每 8 秒导入一次的成本(您可以通过让脚本在运行之间休眠来解决)。因此,很难击败优雅的 3 行解决方案。

【讨论】:

  • 哇,感谢您提供详细信息!我最终选择了一个不同的解决方案,但这是我可以用来学习更多关于 python 的东西。
【解决方案3】:

既然您提出了 pandas 解决方案,这里有一种可能性:

import pandas as pd
import numpy as np

l = [1,4,11,9,14,21]
n = len(l)

df = pd.DataFrame(l, columns = ["l"]).reindex(np.linspace(0, n-1, 2*n-1)).interpolate().astype(int)

print(df)

不过,这感觉不必要的复杂。我在 pandas 中标记,所以更熟悉 pandas 功能的人可以看到它。

【讨论】:

  • 感谢熊猫示例。我同意你的评论——这可能比仅仅计算数学更复杂,这就是我最终要做的。
猜你喜欢
  • 2019-02-25
  • 1970-01-01
  • 1970-01-01
  • 2020-12-28
  • 2023-03-06
  • 2013-07-25
  • 2022-01-04
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多