我正在寻找一种优雅的方法来计算离散卷积的“乘积”，而不是求和。
这是一个离散卷积的公式：
在这种情况下我们可以使用：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
【问题讨论】：
f,g) 始终包含整数，则无论算法如何，您都可以通过将它们转换为 uint16 等来获得一些速度（假设您可以排除溢出）。
标签： arrays matlab multiplication convolution sliding-window
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))，以便缩放后的 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
【讨论】：
f 包含诸如4 5 10 而不是rand 之类的值，cumprod(f) 的结果将不会那么可靠
另一个解决方案部分受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)
但这里不是一个向量，而是两个向量 f 和 g。
所需的公式可以重新表述为：
八度音阶：
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.
【讨论】：
..., 'valid') 添加到conv 调用，并将..., 'Endpoints', 'discard'); 添加到movsum 调用。
这是一种避免循环的方法：
ind = bsxfun(@plus, 0:numel(f)-numel(g), (1:numel(g)).');
x = prod(bsxfun(@times, f(ind), g(:)), 1);
它的工作原理如下。 ind 表示索引的滑动窗口，每一列对应一个位移。例如，如果g 的大小为3，则矩阵ind 将是
1 2 3 4 ...
2 3 4 5 ...
3 4 5 6 ...
这用于索引f。结果乘以g（通过广播）作为一列。最后计算每列元素的乘积。
【讨论】：