测试一个滑动区间内的中点是否为最大值

Question

lengthData = 1500;
data = rand(lengthData,1);
result = zeros(floor(lengthData/2),lengthData);
for i = 1:floor(lengthData/2)-1
    for j = 1+i:length(data)-i
        if data(j) == max(data(j-i:j+i))
            result(i,j)=1;
        end
    end
end

我有一个指定lengthData 的数据点存储在变量data 中。现在我试图在指定的间隔中找到最大值，其中间隔长度正在增加。

一个快速的观察是，随着间隔大小的增加，result 中特定行中的条目数正在减少。以下是sum(result,2) 的图表，用于验证代码是否按预期工作。

但是，此代码需要很长时间才能执行。对于较大的 lengthData 值（约 6000），所用时间几乎为 22 秒（高于 0.4 秒，lengthData 为 1500）。

是否有另一种方法可以实现我的逻辑或对其进行矢量化，以某种方式加快它的速度？

Answer 1

基于this comment：

在我的代码实现中，Maxima 出现在窗口的中间。可能是我在问题描述中不清楚。我正在寻找各种最大值。只有一个全局最大值，但有许多局部最大值。现在，只有当它们大于其附近的数字时，才会存在局部最大值。我认为这个附近是移动窗口的长度。当滑动窗口为 5 时，一个特定的值可能是局部最大值，但当滑动窗口增加到 10 时，它可能不是局部最大值。(...)

我建议看一下findpeaks()，它正是这样做的：找到局部最大值。然后，您的“滑动窗口长度”将被合并到 'MinPeakDistance' 名称-值对中，您在问题中提到的最小高度将由 'MinPeakHeight' 给出。

Answer 2

我已经使用vectorization solution 和findpeaks() 运行了代码。此功能为我提供了具有更好时间性能的预期结果，请参阅以下基准：

tic;
lengthData=6000;
data=rand(lengthData,1);
result=zeros(floor(lengthData/2),lengthData);
for i=1:floor(lengthData/2)-1
    for j=1+i:length(data)-i
        if(data(j)==max(data(j-i:j+i)))
            result(i,j)=1;
        end
    end
end
toc;

tic;
result1=zeros(floor(lengthData/2),lengthData);
for i=1:floor(lengthData/2)-1
    del = 2*i+1;
    ind_arr = [];
    ind_arr(:,1) = 1:lengthData+1-del;
    ind_arr(:,2:del) = 1;
    ind_arr = cumsum(ind_arr,2);
    data_arr = data(ind_arr);
    [~,max_ind] = max(data_arr,[],2);
    result1(i,1+i:length(data)-i) = max_ind==(i+1);
end
toc;


tic;
result2=zeros(floor(lengthData/2),lengthData);

for i=1:floor(lengthData/2)-1

    [~,locs]=findpeaks(data,'MinPeakDistance',i);
    result2(i,1:length(locs))=locs;
end
toc;

结果：

Elapsed time is 23.317170 seconds.   % Original
Elapsed time is 312.340804 seconds.  % Vectorized
Elapsed time is 3.053548 seconds.    % findpeaks()

我使用 MATLAB R2018a，在 Xeon 3.4 处理器上使用 32GB RAM。

Answer 3

我能够使用以下方法对内部循环进行矢量化...

lengthData=1500;
data=rand(lengthData,1);
result=zeros(floor(lengthData/2),lengthData);

for i=1:floor(lengthData/2)-1

    del = 2*i+1;
    ind_arr = [];
    ind_arr(:,1) = 1:lengthData+1-del;
    ind_arr(:,2:del) = 1;
    ind_arr = cumsum(ind_arr,2);

    data_arr = data(ind_arr);

    [~,max_ind] = max(data_arr,[],2);

    result(i,1+i:length(data)-i) = max_ind==(i+1);

end

可能有一种更简洁的方法来构建ind_arr，但这似乎会为您的原始代码提供匹配的结果。

通过 R2014b 中的分析器运行它会显示

  time   calls  line
             1    1 lengthData=1500; 
             1    2 data=rand(lengthData,1); 
             1    3 result=zeros(floor(lengthData/2),lengthData); 
             1    4 result2 = result; 
                  5 
             1    6 t1ID = tic; 
                  7 
             1    8 for i=1:floor(lengthData/2)-1 
           749    9     for j=1+i:length(data)-i 
  3.66  561750   10         if(data(j)==max(data(j-i:j+i))) 
< 0.01    4276   11             result(i,j)=1; 
          4276   12         end 
  0.70  561750   13     end 
           749   14 end 
                 15 
             1   16 fprintf('Original Time: %g\n',toc(t1ID)); 
                 17 
             1   18 t2ID = tic; 
                 19 
             1   20 for i=1:floor(lengthData/2)-1 
                 21     
           749   22     del = 2*i+1; 
  0.02     749   23     ind_arr = []; 
< 0.01     749   24     ind_arr(:,1) = 1:lengthData+1-del; 
  0.45     749   25     ind_arr(:,2:del) = 1; 
  0.51     749   26     ind_arr = cumsum(ind_arr,2); 
                 27 
  1.42     749   28     data_arr = data(ind_arr); 
                 29 
  0.52     749   30     [~,max_ind] = max(data_arr,[],2); 
                 31 
  0.05     749   32     result2(i,1+i:length(data)-i) = max_ind==(i+1); 
                 33     
           749   34 end 
                 35 
             1   36 fprintf('Single Loop Time: %g\n',toc(t2ID));

Original Time: 4.53376
Single Loop Time: 2.81267