将函数应用于滚动 window
Apply function to rolling window
假设我有一长串 A 个值（比如长度 1000），我想计算 100 对中的 std，即我想计算 std(A(1:100)), std(A(2:101)), std(A(3:102)), ..., std(A(901:1000))。
在 Excel/VBA 中，可以通过编写例如=STDEV(A1:A100) 一个单元格，然后一口气往下填。现在我的问题是，如何在 Matlab 中有效地完成这一任务而不必使用任何昂贵的 for 循环。
编辑：是否也可以对时间序列列表执行此操作，例如当 A 的维度为 1000 x 4（即 4 个长度为 1000 的时间序列）？输出矩阵的尺寸应为 901 x 4。
你所做的基本上是过滤操作。
如果您有权访问图像处理工具箱，
stdfilt(A,ones(101,1)) %# assumes that data series are in columns
会成功（无论 A 的维数如何）。请注意，如果您还可以访问并行计算工具箱，则可以在 GPU 上进行类似 运行 的过滤操作，尽管您的问题可能太小而无法产生明显的加速。
注意：有关最快的解决方案，请参阅
所以首先为此使用 for 循环（特别是如果这些是您的实际尺寸）确实不会很昂贵。除非您使用的是非常旧版本的 Matlab，否则 JIT 编译器（当然还有预分配）使 for 循环成本低廉。
其次 - 你试过 for 循环了吗？因为在过早地开始优化之前，您真的应该先尝试简单的实现。
第三 - arrayfun 可以使它成为一个单行，但它基本上只是一个带有额外开销的 for 循环，如果你真的关心速度，它很可能比 for 循环慢。
最后是一些代码：
n = 1000;
A = rand(n,1);
l = 100;
for 循环（体积不大，可能很高效）：
S = zeros(n-l+1,1); %//Pre-allocation of memory like this is essential for efficiency!
for t = 1:(n-l+1)
S(t) = std(A(t:(t+l-1)));
end
矢量化（内存效率低下！）解决方案：
[X,Y] = meshgrid(1:l)
S = std(A(X+Y-1))
一个可能更好的矢量化解决方案（和一个单行）但仍然内存效率低下：
S = std(A(bsxfun(@plus, 0:l-1, (1:l)')))
请注意，使用所有这些方法，您可以将 std 替换为任何函数，只要它适用于矩阵的列（这是 Matlab 中的标准）
进入二维：
要进入 2D，我们需要进入 3D
n = 1000;
k = 4;
A = rand(n,k);
l = 100;
ind = bsxfun(@plus, permute(o:n:(k-1)*n, [3,1,2]), bsxfun(@plus, 0:l-1, (1:l)')); %'
S = squeeze(std(A(ind)));
M = squeeze(mean(A(ind)));
%// etc...
或
[X,Y,Z] = meshgrid(1:l, 1:l, o:n:(k-1)*n);
ind = X+Y+Z-1;
S = squeeze(std(A(ind)))
M = squeeze(mean(A(ind)))
%// etc...
或
ind = bsxfun(@plus, 0:l-1, (1:l)'); %'
for t = 1:k
S = std(A(ind));
M = mean(A(ind));
%// etc...
end
OR（取自 - 在他的回答中请注意，他展示了一个比这个简单循环更快的替代方法）
S = zeros(n-l+1,k);
M = zeros(n-l+1,k);
for t = 1:(n-l+1)
S(t,:) = std(A(k:(k+l-1),:));
M(t,:) = mean(A(k:(k+l-1),:));
%// etc...
end
几个函数可以在 Matlab 中高效地完成这项工作。
一方面，您可以使用 colfilt 或 nlfilter 等函数，对滑动块执行计算。 colfilt 比 nlfilter 更有效，但只能在块内元素的顺序无关紧要时使用。以下是如何在您的数据上使用它：
S = colfilt(A, [100,1], 'sliding', @std);
或
S = nlfilter(A, [100,1], @std);
在你的例子中，你可以清楚地看到性能的差异。但是有一个技巧：两个函数都对输入数组进行填充，以便输出向量与输入数组具有相同的大小。要仅获取输出向量的相关部分，您需要跳过前 floor((100-1)/2) = 49 个元素，并取 1000-100+1 个值。
S(50:end-50)
但还有另一种解决方案，接近colfilt，效率更高。 colfilt 调用 col2im 将输入向量重塑为一个矩阵，在该矩阵上它在每个不同的列上应用给定的函数。这会将大小为 [1000,1] 的输入向量转换为大小为 [100,901] 的矩阵。但是 colfilt 用 0 或 1 填充输入数组，您不需要它。所以你可以 运行 colfilt 没有填充步骤，然后在每一列上应用 std 这很容易，因为 std 应用于矩阵 returns 行向量列的标准。最后，如果需要，将其转置以获得列向量。简而言之，一行：
S = std(im2col(X,[100 1],'sliding')).';
备注：如果要应用更复杂的功能，请参见colfilt的代码，第144和147行（适用于v2013b）。
如果您关心 for 循环的速度，您可以通过将向量折叠成一个数组（使用 reshape）来大大减少循环迭代的次数，其中的列数为您要应用函数的元素。
这将使 Matlab 和 JIT 通过在数组的每一列上计算您的函数来执行优化（在大多数情况下他们比我们做得更好）。
然后您重塑数组的偏移版本并执行相同的操作。您仍然需要一个循环，但迭代次数仅为 l（因此在您的示例中为 100），而不是经典 for 循环中的 n-l+1=901（一个window一次）。
完成后，将结果数组重塑为向量，然后您仍然需要手动计算最后一个window，但总体而言仍然快得多。
采用与 Dan 相同的输入符号：
n = 1000;
A = rand(n,1);
l = 100;
它会变成这样：
width = (n/l)-1 ; %// width of each line in the temporary result array
tmp = zeros( l , width ) ; %// preallocation never hurts
for k = 1:l
tmp(k,:) = std( reshape( A(k:end-l+k-1) , l , [] ) ) ; %// calculate your stat on the array (reshaped vector)
end
S2 = [tmp(:) ; std( A(end-l+1:end) ) ] ; %// "unfold" your results then add the last window calculation
如果我 tic ... toc 完整的循环版本和折叠的版本，我得到这个平均结果：
Elapsed time is 0.057190 seconds. %// windows by window FOR loop
Elapsed time is 0.016345 seconds. %// "Folded" FOR loop
我知道 tic/toc 不是完美计时的方法，但我的 matlab 版本没有 timeit 函数。此外，差异显着足以表明存在改进（尽管无法通过此方法精确量化）。我当然删除了第一个 运行 并检查了结果是否与不同的矩阵大小一致。
关于您的 "one liner" 请求，我建议您将这段代码包装成一个函数，如下所示：
function out = foldfunction( func , vec , nPts )
n = length( vec ) ;
width = (n/nPts)-1 ;
tmp = zeros( nPts , width ) ;
for k = 1:nPts
tmp(k,:) = func( reshape( vec(k:end-nPts+k-1) , nPts , [] ) ) ;
end
out = [tmp(:) ; func( vec(end-nPts+1:end) ) ] ;
在你的主代码中哪一个允许你在一行中调用它：
S = foldfunction( @std , A , l ) ;
这种格式的另一个好处是，您可以将完全相同的子函数用于其他统计函数。例如，如果您想要 windows 的 "mean"，只需更改 func 参数即可调用相同的方法：
S = foldfunction( @mean , A , l ) ;
唯一的限制，因为它只适用于向量作为输入，但通过一些返工，它也可以将数组作为输入。
为了最大程度地减少操作次数，您可以利用以下事实，即标准偏差可以计算为涉及第二个和第一个时刻的差异 、
并且通过累积和（使用cumsum）有效地获得滚动window上的力矩：
A = randn(1000,4); %// random data
N = 100; %// window size
c = size(A,2);
A1 = [zeros(1,c); cumsum(A)];
A2 = [zeros(1,c); cumsum(A.^2)];
S = sqrt( (A2(1+N:end,:)-A2(1:end-N,:) ...
- (A1(1+N:end,:)-A1(1:end-N,:)).^2/N) / (N-1) ); %// result
基准测试
这是与基于循环的解决方案的比较，使用 timeit. The loop approach is as in but adapted to the 2D case, exploting the fact that std 以矢量化方式沿每一列工作。
%// File loop_approach.m
function S = loop_approach(A,N);
[n, p] = size(A);
S = zeros(n-N+1,p);
for k = 1:(n-N+1)
S(k,:) = std(A(k:(k+N-1),:));
end
%// File bsxfun_approach.m
function S = bsxfun_approach(A,N);
[n, p] = size(A);
ind = bsxfun(@plus, permute(0:n:(p-1)*n, [3,1,2]), bsxfun(@plus, 0:n-N, (1:N).')); %'
S = squeeze(std(A(ind)));
%// File cumsum_approach.m
function S = cumsum_approach(A,N);
c = size(A,2);
A1 = [zeros(1,c); cumsum(A)];
A2 = [zeros(1,c); cumsum(A.^2)];
S = sqrt( (A2(1+N:end,:)-A2(1:end-N,:) ...
- (A1(1+N:end,:)-A1(1:end-N,:)).^2/N) / (N-1) );
%// Benchmarking code
clear all
A = randn(1000,4); %// Or A = randn(1000,1);
N = 100;
t_loop = timeit(@() loop_approach(A,N));
t_bsxfun = timeit(@() bsxfun_approach(A,N));
t_cumsum = timeit(@() cumsum_approach(A,N));
disp(' ')
disp(['loop approach: ' num2str(t_loop)])
disp(['bsxfun approach: ' num2str(t_bsxfun)])
disp(['cumsum approach: ' num2str(t_cumsum)])
disp(' ')
disp(['bsxfun/loop gain factor: ' num2str(t_loop/t_bsxfun)])
disp(['cumsum/loop gain factor: ' num2str(t_loop/t_cumsum)])
结果
我正在使用 Matlab R2014b，Windows 7 64 位，双核处理器，4 GB 内存：
4 列大小写：
loop approach: 0.092035
bsxfun approach: 0.023535
cumsum approach: 0.0002338
bsxfun/loop gain factor: 3.9106
cumsum/loop gain factor: 393.6526
Single-column 案例：
loop approach: 0.085618
bsxfun approach: 0.0040495
cumsum approach: 8.3642e-05
bsxfun/loop gain factor: 21.1431
cumsum/loop gain factor: 1023.6236
所以基于 cumsum 的方法似乎是最快的：比 4 列情况下的循环快 400 倍，在 single-column 情况下快 1000 倍。
假设我有一长串 A 个值（比如长度 1000），我想计算 100 对中的 std，即我想计算 std(A(1:100)), std(A(2:101)), std(A(3:102)), ..., std(A(901:1000))。
在 Excel/VBA 中，可以通过编写例如=STDEV(A1:A100) 一个单元格，然后一口气往下填。现在我的问题是，如何在 Matlab 中有效地完成这一任务而不必使用任何昂贵的 for 循环。
编辑：是否也可以对时间序列列表执行此操作，例如当 A 的维度为 1000 x 4（即 4 个长度为 1000 的时间序列）？输出矩阵的尺寸应为 901 x 4。
你所做的基本上是过滤操作。
如果您有权访问图像处理工具箱，
stdfilt(A,ones(101,1)) %# assumes that data series are in columns
会成功（无论 A 的维数如何）。请注意，如果您还可以访问并行计算工具箱，则可以在 GPU 上进行类似 运行 的过滤操作，尽管您的问题可能太小而无法产生明显的加速。
注意：有关最快的解决方案，请参阅
所以首先为此使用 for 循环（特别是如果这些是您的实际尺寸）确实不会很昂贵。除非您使用的是非常旧版本的 Matlab，否则 JIT 编译器（当然还有预分配）使 for 循环成本低廉。
其次 - 你试过 for 循环了吗？因为在过早地开始优化之前，您真的应该先尝试简单的实现。
第三 - arrayfun 可以使它成为一个单行，但它基本上只是一个带有额外开销的 for 循环，如果你真的关心速度，它很可能比 for 循环慢。
最后是一些代码：
n = 1000;
A = rand(n,1);
l = 100;
for 循环（体积不大，可能很高效）：
S = zeros(n-l+1,1); %//Pre-allocation of memory like this is essential for efficiency!
for t = 1:(n-l+1)
S(t) = std(A(t:(t+l-1)));
end
矢量化（内存效率低下！）解决方案：
[X,Y] = meshgrid(1:l)
S = std(A(X+Y-1))
一个可能更好的矢量化解决方案（和一个单行）但仍然内存效率低下：
S = std(A(bsxfun(@plus, 0:l-1, (1:l)')))
请注意，使用所有这些方法，您可以将 std 替换为任何函数，只要它适用于矩阵的列（这是 Matlab 中的标准）
进入二维：
要进入 2D，我们需要进入 3D
n = 1000;
k = 4;
A = rand(n,k);
l = 100;
ind = bsxfun(@plus, permute(o:n:(k-1)*n, [3,1,2]), bsxfun(@plus, 0:l-1, (1:l)')); %'
S = squeeze(std(A(ind)));
M = squeeze(mean(A(ind)));
%// etc...
或
[X,Y,Z] = meshgrid(1:l, 1:l, o:n:(k-1)*n);
ind = X+Y+Z-1;
S = squeeze(std(A(ind)))
M = squeeze(mean(A(ind)))
%// etc...
或
ind = bsxfun(@plus, 0:l-1, (1:l)'); %'
for t = 1:k
S = std(A(ind));
M = mean(A(ind));
%// etc...
end
OR（取自
S = zeros(n-l+1,k);
M = zeros(n-l+1,k);
for t = 1:(n-l+1)
S(t,:) = std(A(k:(k+l-1),:));
M(t,:) = mean(A(k:(k+l-1),:));
%// etc...
end
几个函数可以在 Matlab 中高效地完成这项工作。
一方面，您可以使用 colfilt 或 nlfilter 等函数，对滑动块执行计算。 colfilt 比 nlfilter 更有效，但只能在块内元素的顺序无关紧要时使用。以下是如何在您的数据上使用它：
S = colfilt(A, [100,1], 'sliding', @std);
或
S = nlfilter(A, [100,1], @std);
在你的例子中，你可以清楚地看到性能的差异。但是有一个技巧：两个函数都对输入数组进行填充，以便输出向量与输入数组具有相同的大小。要仅获取输出向量的相关部分，您需要跳过前 floor((100-1)/2) = 49 个元素，并取 1000-100+1 个值。
S(50:end-50)
但还有另一种解决方案，接近colfilt，效率更高。 colfilt 调用 col2im 将输入向量重塑为一个矩阵，在该矩阵上它在每个不同的列上应用给定的函数。这会将大小为 [1000,1] 的输入向量转换为大小为 [100,901] 的矩阵。但是 colfilt 用 0 或 1 填充输入数组，您不需要它。所以你可以 运行 colfilt 没有填充步骤，然后在每一列上应用 std 这很容易，因为 std 应用于矩阵 returns 行向量列的标准。最后，如果需要，将其转置以获得列向量。简而言之，一行：
S = std(im2col(X,[100 1],'sliding')).';
备注：如果要应用更复杂的功能，请参见colfilt的代码，第144和147行（适用于v2013b）。
如果您关心 for 循环的速度，您可以通过将向量折叠成一个数组（使用 reshape）来大大减少循环迭代的次数，其中的列数为您要应用函数的元素。
这将使 Matlab 和 JIT 通过在数组的每一列上计算您的函数来执行优化（在大多数情况下他们比我们做得更好）。
然后您重塑数组的偏移版本并执行相同的操作。您仍然需要一个循环，但迭代次数仅为 l（因此在您的示例中为 100），而不是经典 for 循环中的 n-l+1=901（一个window一次）。
完成后，将结果数组重塑为向量，然后您仍然需要手动计算最后一个window，但总体而言仍然快得多。
采用与 Dan 相同的输入符号：
n = 1000;
A = rand(n,1);
l = 100;
它会变成这样：
width = (n/l)-1 ; %// width of each line in the temporary result array
tmp = zeros( l , width ) ; %// preallocation never hurts
for k = 1:l
tmp(k,:) = std( reshape( A(k:end-l+k-1) , l , [] ) ) ; %// calculate your stat on the array (reshaped vector)
end
S2 = [tmp(:) ; std( A(end-l+1:end) ) ] ; %// "unfold" your results then add the last window calculation
如果我 tic ... toc 完整的循环版本和折叠的版本，我得到这个平均结果：
Elapsed time is 0.057190 seconds. %// windows by window FOR loop
Elapsed time is 0.016345 seconds. %// "Folded" FOR loop
我知道 tic/toc 不是完美计时的方法，但我的 matlab 版本没有 timeit 函数。此外，差异显着足以表明存在改进（尽管无法通过此方法精确量化）。我当然删除了第一个 运行 并检查了结果是否与不同的矩阵大小一致。
关于您的 "one liner" 请求，我建议您将这段代码包装成一个函数，如下所示：
function out = foldfunction( func , vec , nPts )
n = length( vec ) ;
width = (n/nPts)-1 ;
tmp = zeros( nPts , width ) ;
for k = 1:nPts
tmp(k,:) = func( reshape( vec(k:end-nPts+k-1) , nPts , [] ) ) ;
end
out = [tmp(:) ; func( vec(end-nPts+1:end) ) ] ;
在你的主代码中哪一个允许你在一行中调用它：
S = foldfunction( @std , A , l ) ;
这种格式的另一个好处是，您可以将完全相同的子函数用于其他统计函数。例如，如果您想要 windows 的 "mean"，只需更改 func 参数即可调用相同的方法：
S = foldfunction( @mean , A , l ) ;
唯一的限制，因为它只适用于向量作为输入，但通过一些返工，它也可以将数组作为输入。
为了最大程度地减少操作次数，您可以利用以下事实，即标准偏差可以计算为涉及第二个和第一个时刻的差异 、
并且通过累积和（使用cumsum）有效地获得滚动window上的力矩：
A = randn(1000,4); %// random data
N = 100; %// window size
c = size(A,2);
A1 = [zeros(1,c); cumsum(A)];
A2 = [zeros(1,c); cumsum(A.^2)];
S = sqrt( (A2(1+N:end,:)-A2(1:end-N,:) ...
- (A1(1+N:end,:)-A1(1:end-N,:)).^2/N) / (N-1) ); %// result
基准测试
这是与基于循环的解决方案的比较，使用 timeit. The loop approach is as in std 以矢量化方式沿每一列工作。
%// File loop_approach.m
function S = loop_approach(A,N);
[n, p] = size(A);
S = zeros(n-N+1,p);
for k = 1:(n-N+1)
S(k,:) = std(A(k:(k+N-1),:));
end
%// File bsxfun_approach.m
function S = bsxfun_approach(A,N);
[n, p] = size(A);
ind = bsxfun(@plus, permute(0:n:(p-1)*n, [3,1,2]), bsxfun(@plus, 0:n-N, (1:N).')); %'
S = squeeze(std(A(ind)));
%// File cumsum_approach.m
function S = cumsum_approach(A,N);
c = size(A,2);
A1 = [zeros(1,c); cumsum(A)];
A2 = [zeros(1,c); cumsum(A.^2)];
S = sqrt( (A2(1+N:end,:)-A2(1:end-N,:) ...
- (A1(1+N:end,:)-A1(1:end-N,:)).^2/N) / (N-1) );
%// Benchmarking code
clear all
A = randn(1000,4); %// Or A = randn(1000,1);
N = 100;
t_loop = timeit(@() loop_approach(A,N));
t_bsxfun = timeit(@() bsxfun_approach(A,N));
t_cumsum = timeit(@() cumsum_approach(A,N));
disp(' ')
disp(['loop approach: ' num2str(t_loop)])
disp(['bsxfun approach: ' num2str(t_bsxfun)])
disp(['cumsum approach: ' num2str(t_cumsum)])
disp(' ')
disp(['bsxfun/loop gain factor: ' num2str(t_loop/t_bsxfun)])
disp(['cumsum/loop gain factor: ' num2str(t_loop/t_cumsum)])
结果
我正在使用 Matlab R2014b，Windows 7 64 位，双核处理器，4 GB 内存：
4 列大小写：
loop approach: 0.092035 bsxfun approach: 0.023535 cumsum approach: 0.0002338 bsxfun/loop gain factor: 3.9106 cumsum/loop gain factor: 393.6526Single-column 案例：
loop approach: 0.085618 bsxfun approach: 0.0040495 cumsum approach: 8.3642e-05 bsxfun/loop gain factor: 21.1431 cumsum/loop gain factor: 1023.6236
所以基于 cumsum 的方法似乎是最快的：比 4 列情况下的循环快 400 倍，在 single-column 情况下快 1000 倍。