arrays - 计算离散卷积的 "product"的有效方法

Question

我正在寻找一种优雅的方法来计算离散卷积的“乘积”，而不是总和。

这是一个离散卷积的公式:

在这种情况下我们可以使用:conv(x,y)

现在我想实现这些操作

当然我可以使用循环，但我正在寻找一种技巧来线性化此操作。

示例:

f = [2 4 3 9 7 1]
g = [3 2 1]

dist = length(g)-1;

for ii = 1:length(f)-dist
    x(ii) = prod(f(ii:ii+dist).*g)
end

x =

144    648   1134    378

Answer 1

cumprod解决方案:(非常高效)

>> pf = cumprod(f);
>> x = prod(g).*pf(numel(g):end)./[1 pf(1:(end-numel(g)))]

x =

         144         648        1134         378

这首先使用 cumprod 获取 f 的累积乘积.通过将每个元素除以它之前的 3 个元素的累积乘积，我们得到每个 numel(g)-wide 滑动窗口沿 f 的乘积。然后乘以 g 的元素的乘积。

注意:当 f 有很多元素或极值(大或小)时，您可能会在执行累积乘积时遇到准确性或下溢/溢出问题.缓解这种情况的一种可能方法是在累积乘积之前对 f 应用缩放，然后在之后撤消它:

c = ...set a scaling factor...
pf = cumprod(f./c);
x = prod(c.*g).*pf(numel(g):end)./[1 pf(1:(end-numel(g)))];

c 的选择可能类似于 mean(abs(f)) 或 max(abs(f)) 这样scaled f 给出的结果有更好的界限(即值更接近 1)。这不会明显改变下面的计时结果。

---

hankel 解决方案:(效率不高但仍然有趣)

>> x = prod(g).*prod(hankel(f(1:numel(g)), f(numel(g):end)))

x =

         144         648        1134         378

调用hankel创建一个矩阵，其中每一列都包含其中一个 numel(g) 范围的滑动窗口的内容。将每列的乘积乘以 g 的元素的乘积即可得到答案。但是，对于大型向量 f 和/或 g，这可能涉及大量额外计算并使用大量内存。

---

对结果进行计时:

我测试了 6 个解决方案(你问题中的循环，来自 rahnema 的 2 个解决方案(conv(log) 和 movsum(log))，来自 Luis 的 bsxfun 解决方案，以及我的 cumprod 和 hankel 解决方案)使用 f = rand(1,1000000); 和g = rand(1,100); 并平均超过 40 次迭代。这是我得到的(运行 Windows 7 x64、16 GB RAM、MATLAB R2016b):

solution    | avg. time (s)
------------+---------------
loop        |       1.10671
conv(log)   |       0.04984
movsum(log) |       0.03736
bsxfun      |       1.20472
cumprod     |       0.01469
hankel      |       1.17704