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

这是 离散卷积 的公式:

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

在这种情况下,我们可以使用:conv(x,y)
现在我想实现这些操作

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

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

示例:

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

最佳答案

另一个解决方案部分受到 Dev-iL answer 的启发,以解决相同的问题

exp(sum(log(g))+conv(log(f),[1 1 1],'valid'))

或者
exp(sum(log(g))+movsum(log(f),numel(g),'Endpoints', 'discard'))

由于 exp(sum(log(x))) = prod(x)
但是这里不是一个向量,我们有两个向量 fg

所需的公式可以重新表述为:

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

Octave 计时:
f= rand(1,1000000);
g= rand(1,100);

disp('----------EXP(LOG)------')
tic
    exp(sum(log(g))+conv(log(f),ones(1,numel(g))));
toc

disp('----------BSXFUN------')
tic
    ind = bsxfun(@plus, 0:numel(f)-numel(g), (1:numel(g)).');
    x = prod(bsxfun(@times, f(ind), g(:)),1);
toc
disp('----------HANKEL------')
tic
    prod(g)*prod(hankel(f(1:numel(g)), f(numel(g):end)));
toc

disp('----------CUMPROD-----')
tic
    pf = cumprod(f);
    x = prod(g)*pf(numel(g):end)./[1 pf(1:(end-numel(g)))];
toc

结果:
----------EXP(LOG)------%rahnema1
Elapsed time is 0.211445 seconds.
----------BSXFUN--------%Luis Mendo
Elapsed time is 1.94182 seconds.
----------HANKEL--------%gnovice
Elapsed time is 1.46593 seconds.
----------CUMPROD-------%gnovice
Elapsed time is 0.00748992 seconds.

10-08 01:53