我正在使用移动平均类型的平滑操作卷积信号块,但有问题的填充错误,影响下游的计算.

如果绘制,这是块边界上引起的问题.

enter image description here

为了解决这个问题,我将前一个块中的N个样本添加到数组的beginning个.

问题是下面的卷积操作是从原始数组的开头还是结尾删除样本?

arr = np.array([0.9, 0.9, 0.9, 0.1, 0.1, 0.1, 1.0, 0.1, 0.1, 0.1, 0.9, 0.9, 0.9, ])
print(f"length input array {len(arr)}")
result = np.convolve(arr, [1]*3, 'valid')/3
print(f"length result array {len(result)}")
print(result)
plt.plot(result)

我设置了上面的例子,试图确认将样本添加到数组的开头位置是正确的,但它没有帮助.输入长度为13,输出长度为11.

如果有重叠,卷积看起来像:

result = np.convolve(np.concatenate(previous_samples[-2:], arr), [1]*3, 'valid')/3

推荐答案

I assume (1) the reader has a basic understanding of what a convolution is and (2) a kernel with an odd number of elements is used.

Short answer:

  • 使用np.convolve(…, mode="valid").

  • 在卷积之前,

    • .在它们的开始与前一个块的最后(len(kernel) - 1) / 2个元素,
    • .在它们的结尾与下一个块的前(len(kernel) - 1) / 2个元素,

    其中kernel指的是你的滑动内核(你的问题中的[1]*3).

Long answer:

由于卷积可以被认为是一个具有滑动内核的操作,所以问题是,正如你所意识到的,边界处发生了什么.

numpy.convolve()函数有3种模式来处理边界,如图所示:

  • mode="valid"确保滑动核从不离开信号.然而,这意味着结果将比信号具有更少的元素;即少len(kernel) - 1个元素.因此,如果你有一个mode="valid"个元素的信号和一个5个元素的内核,你的结果将保持mode="valid"—(5—1)= 96个元素.操作在何处"移除"缺失元素的问题可能有点误导;然而,考虑结果的一种方式是将其值写入滑动核的中心,(这对于奇数个元素的对称核是有意义的,就像你的问题中的情况一样),所以通过这种解释,我们可以得到元素的对称移除:(len(kernel) - 1) / 2个元素在信号的开始,(len(kernel) - 1) / 2个元素在信号的结束.我们将在下面的代码示例中看到这是一个有意义的解释.
  • mode="same"确保所得信号具有与输入信号相同数量的值.然而,这意味着我们必须"弥补"值,以抵消我们在mode="valid"中看到的元素较少的影响.因此,在这种模式下,在卷积操作之前,信号在开始和结束时都用零填充(len(kernel) - 1) / 2.
  • 最后,mode="full"确保产生信号和滑动核之间的所有潜在重叠组合.这意味着(a)我们必须"弥补"更多的值,(b)结果将大于输入信号.这里的"合成"值再次为零,结果将包含len(signal) + len(kernel) - 1个元素,因此对于前面的示例来说,mode="full" + 5—1 = 104个元素.

如果你以块处理你的信号,mode="valid"是一种方法:你可以确保你得到的卷积块具有与输入块相同的长度,而不需要"弥补"值(即,不需要零填充),通过填充前一个块的尾部值和下一个块的前导值.除了第一个和最后一个块之外,这对所有块都是如此,在这种情况下,你必须自己决定如何处理边界情况.在所有其他情况下,填充前一个和后一个块的值小于内核长度的一半.

下面的代码演示了所描述的块填充产生的结果与根本不对信号进行块处理完全相同.请注意,我将图中的卷积信号移位了(len(kernel) - 1) / 2个值,以解释所描述的缺失值,信号的最开始为mode="valid".

import matplotlib.pyplot as plt
import numpy as np

len_signal, len_kernel = 100, 5
num_chunks = 10

rand = np.random.default_rng(seed=42)
signal = np.cumsum(rand.normal(size=len_signal))
kernel = np.ones(len_kernel) / len_kernel

# Produce convolution result without chunking for reference
signal_convolved = np.convolve(signal, kernel, mode="valid")

# Split into chunks, then pad them: (1) start with `(len_kernel - 1) / 2`
# last values from preceding chunk (except first chunk), (2) end with
# `(len_kernel - 1) / 2` first values from following chunk (except last chunk)
len_arm = (len_kernel - 1) // 2  # Length of a "kernel arm"
chunks_convolved = []
for i, chunk in enumerate(chunks := np.split(signal, num_chunks)):
    p1 = [] if i == 0 else chunks[i - 1][-len_arm:]  # Start padding
    p2 = [] if i == len(chunks) - 1 else chunks[i + 1][:len_arm]  # End padding
    chunk_conv = np.convolve(np.r_[p1, chunk, p2], kernel, mode="valid")
    chunks_convolved = np.r_[chunks_convolved, chunk_conv]  # Append

assert np.allclose(signal_convolved, chunks_convolved)  # Check: same result?
# Show: convolved w/o chunking (left) vs chunked, padded, and convolved (right)
x_c = np.arange(len(signal_convolved)) + len_arm   # Right-shift conv. result
plt.subplot(121)
plt.plot(signal, "--"); plt.plot(x_c, signal_convolved); plt.title("no chunks")
plt.subplot(122)
plt.plot(signal, "--"); plt.plot(x_c, chunks_convolved); plt.title("chunks")
plt.show()

As you can see (and as we checked with assert), both results are the same: convolution result without (left) and with (right) chunking and padding

作为最后一个侧记:scipy's convolve1d()提供了除零填充之外的更多选项,用于处理总信号开始和结束时的剩余边界情况,因此可能是一个替代方案,具体取决于用例.

Python相关问答推荐

添加包含中具有任何值的其他列的计数的列

Python 3.12中的通用[T]类方法隐式类型检索

在内部列表上滚动窗口

为什么我的Python代码在if-else声明中的行之前执行if-else声明中的行?

try 在树叶 map 上应用覆盖磁贴

为什么符号没有按顺序添加?

PyQt5,如何使每个对象的 colored颜色 不同?'

Stacked bar chart from billrame

创建可序列化数据模型的最佳方法

多处理队列在与Forking http.server一起使用时随机跳过项目

如何合并两个列表,并获得每个索引值最高的列表名称?

如何从列表框中 Select 而不出错?

try 检索blob名称列表时出现错误填充错误""

如何创建引用列表并分配值的Systemrame列

在二维NumPy数组中,如何 Select 内部数组的第一个和第二个元素?这可以通过索引来实现吗?

用两个字符串构建回文

判断Python操作:如何从字面上得到所有decorator ?

Python Mercury离线安装

根据过滤后的牛郎星图表中的数据计算新系列

极点:在固定点扩展窗口