如何（高效地）计算向量的移动平均？

Question

我有一个向量，想要计算它的移动平均值（使用窗口宽度为5）。

例如，如果所涉及的向量是[1,2,3,4,5,6,7,8]，那么

- 结果向量的第一项应该是[1,2,3,4,5]中所有项的总和（即15）; - 结果向量的第二项应该是[2,3,4,5,6]中所有项的总和（即20）; - 等等。

最终，结果向量应该是[15,20,25,30]。我该怎么做？

- user3034281

---

请查看conv函数。 - jub0bs

Answer 1

conv函数非常适合你：

>> x = 1:8;
>> y = conv(x, ones(1,5), 'valid')

y =
    15    20    25    30

基准测试

三个答案，三种不同的方法... 这里是一个快速的基准测试(不同的输入大小，固定窗口宽度为5)，使用timeit; 如果您认为需要改进，请随意在评论中提出。

conv成为最快的方法； 它比coin的方法（使用filter） 快两倍左右，比Luis Mendo的方法（使用cumsum） 快四倍左右。

这里是另一个基准测试（固定输入大小为1e4，不同的窗口宽度）。 在这里，Luis Mendo的cumsum方法成为明显的赢家，因为它的复杂性主要取决于输入长度，并且对窗口的宽度不敏感。

结论

总之，您应该

• 如果您的窗口相对较小，请使用conv方法，
• 如果您的窗口相对较大，请使用cumsum方法。

代码(用于基准测试)

function benchmark

    clear all
    w = 5;                 % moving average window width
    u = ones(1, w); 
    n = logspace(2,6,60);  % vector of input sizes for benchmark
    t1 = zeros(size(n));   % preallocation of time vectors before the loop
    t2 = t1;
    th = t1;

    for k = 1 : numel(n)

        x = rand(1, round(n(k)));  % generate random row vector

        % Luis Mendo's approach (cumsum)
        f = @() luisMendo(w, x);
        tf(k) = timeit(f);

        % coin's approach (filter)
        g = @() coin(w, u, x);
        tg(k) = timeit(g);

        % Jubobs's approach (conv)
        h = @() jubobs(u, x);
        th(k) = timeit(h);
    end

    figure
    hold on
    plot(n, tf, 'bo')
    plot(n, tg, 'ro')
    plot(n, th, 'mo')
    hold off
    xlabel('input size')
    ylabel('time (s)')
    legend('cumsum', 'filter', 'conv')

end

function y = luisMendo(w,x)
    cs = cumsum(x);
    y(1,numel(x)-w+1) = 0; %// hackish way to preallocate result
    y(1) = cs(w);
    y(2:end) = cs(w+1:end) - cs(1:end-w);
end

function y = coin(w,u,x)
    y = filter(u, 1, x);
    y = y(w:end);
end

function jubobs(u,x)
    y = conv(x, u, 'valid');
end

---

function benchmark2

    clear all
    w = round(logspace(1,3,31));    % moving average window width 
    n = 1e4;  % vector of input sizes for benchmark
    t1 = zeros(size(n));   % preallocation of time vectors before the loop
    t2 = t1;
    th = t1;

    for k = 1 : numel(w)
        u = ones(1, w(k));
        x = rand(1, n);  % generate random row vector

        % Luis Mendo's approach (cumsum)
        f = @() luisMendo(w(k), x);
        tf(k) = timeit(f);

        % coin's approach (filter)
        g = @() coin(w(k), u, x);
        tg(k) = timeit(g);

        % Jubobs's approach (conv)
        h = @() jubobs(u, x);
        th(k) = timeit(h);
    end

    figure
    hold on
    plot(w, tf, 'bo')
    plot(w, tg, 'ro')
    plot(w, th, 'mo')
    hold off
    xlabel('window size')
    ylabel('time (s)')
    legend('cumsum', 'filter', 'conv')

end

function y = luisMendo(w,x)
    cs = cumsum(x);
    y(1,numel(x)-w+1) = 0; %// hackish way to preallocate result
    y(1) = cs(w);
    y(2:end) = cs(w+1:end) - cs(1:end-w);
end

function y = coin(w,u,x)
    y = filter(u, 1, x);
    y = y(w:end);
end

function jubobs(u,x)
    y = conv(x, u, 'valid');
end

- jub0bs

---

看看R2016b：故事基本相同。但是，R2016a引入了movmean内置函数。对于小窗口大小的情况，其性能与filter方法大致相当（虽然稍微有些嘈杂）。对于大窗口大小的情况，其性能与cumsum相当。 - sco1

Answer 2

另一个可能的方法是使用cumsum。这种方法可能需要比conv更少的操作：

x = 1:8
n = 5;
cs = cumsum(x);
result = cs(n:end) - [0 cs(1:end-n)];

为了节省一点时间，您可以将最后一行代码替换为：

%// clear result
result(1,numel(x)-n+1) = 0; %// hackish way to preallocate result
result(1) = cs(n);
result(2:end) = cs(n+1:end) - cs(1:end-n);

- Luis Mendo

---

@Jubobs 为什么要使用 u = ones(1, 6)? 不应该是 u = ones(1, w) 吗？你计算 y 的三个过程应该给出相同的大小。另外，为了可靠的时间测量，请使用 timeit。 - Luis Mendo

@Jubobs 如果你更新了你的基准测试（顺便说一句，已经为你的努力点赞+1），能否使用我的第二个版本？ - Luis Mendo

是的，那个 6 是打错了；我不确定它是怎么出现的。稍后我会重新运行基准测试。我现在没有访问 MATLAB 的权限，但我会在有机会时使用 timit 进行测试。 - jub0bs

@Jubobs 我明白了。 cumsum 的优势只在较大的 w 值时才明显。例如，对于 w=50，它确实是最快的方法（在我的机器上）。做得好的基准测试！ - Luis Mendo

是的，现在你已经向我解释清楚了，一切都变得很清晰。太好了！总的结论是，如果你的窗口较大，最好使用cumsum，但如果你的窗口较窄，应该使用conv。 - jub0bs

Answer 3

如果您想保留输入向量的大小，我建议使用filter

>> x = 1:8;
>> y = filter(ones(1,5), 1, x)

y =
     1     3     6    10    15    20    25     30

>> y = (5:end)

y =
     15    20    25    30

- nowox

---

请注意，您正在错误地使用 filter。语法是 filter(b,a,x)，因此您应该使用 filter(ones(1,5), 1, x)。之后，您还应该丢弃结果的前4个元素。 - jub0bs