计算离散卷积"乘积"的有效方法

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

这首先是累计f使用的产品cumprod.通过将每个元素除以之前的累积乘积3个元素,我们得到每个numel(g)滑动窗口的乘积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))以便缩放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个解决方案(在你的问题中循环,从2个解决方案rahnema(conv(log)和movsum(log))时,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