信号间隙宽度的直方图（Matlab）答案

【问题标题】：histogram of signals gaps width (Matlab)信号间隙宽度的直方图（Matlab）
【发布时间】：2021-05-25 09:44:20
【问题描述】：

我正在寻找算法（有效 + 向量化）如何通过以下方式找到间隙 (NaN) 宽度的直方图：

信号由 (Nsamples x Nsig) 数组表示
信号中的间隙由 NaN 编码
间隙宽度：是信号中连续 NaN 的数量
间隙宽度直方图：是信号中具有特定宽度的间隙的频率

并且满足以下条件：

[Nsamples,Nsig ]= size(signals)
isequal(size(signals),size(gapwidthhist)) % true
isequal(sum(gapwidthhist.*(1:Nsamples)',1),sum(isnan(signals),1)) % true

当然，gapwidthhist 的压缩形式（由两个单元格表示：“gapwidthhist_compressed_widths”和“gapwidthhist_compressed_freqs”）也是必需的。

例子：

signals = [1.1 NaN NaN NaN  -1.4 NaN 8.3 NaN NaN NaN  NaN 1.5 NaN NaN; % signal No. 1
           NaN 2.2 NaN 4.9   NaN 8.2 NaN NaN NaN NaN  NaN 2.4 NaN NaN]' % signal No. 2
gapwidthhist = [1 1 1 1 0 0 0 0 0 0 0 0 0 0; % gap histogram for signal No. 1
                3 1 0 0 1 0 0 0 0 0 0 0 0 0]' % gap histogram for signal No. 2

其中整数直方图 bin（间隙宽度）为 1:Nsamples (Nsamples=14)。

对应的压缩间隙直方图如下所示：

gapwidthhist_compressed_widths =  cell(1,Nsig)
gapwidthhist_compressed_widths =
  1×2 cell array
    {[1 2 3 4]}    {[1 2 5]}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gapwidthhist_compressed_freqs = cell(1, Nsig)
gapwidthhist_compressed_freqs =
  1×2 cell array
    {[1 1 1 1]}    {[3 1 1]}

典型问题维度：

Nsamples = 1e5 - 1e6
Nsig = 1e2 - 1e3.

提前感谢您的帮助。

添加备注：我目前最好的解决方案是以下代码：

signals = [1.1 NaN NaN NaN  -1.4 NaN 8.3 NaN NaN NaN  NaN 1.5 NaN NaN; % signal No. 1
           NaN 2.2 NaN 4.9   NaN 8.2 NaN NaN NaN NaN  NaN 2.4 NaN NaN; % signal No. 2
           1 NaN NaN NaN NaN  NaN NaN NaN NaN  NaN NaN NaN NaN NaN]' % signal No. 3
[numData, numSignals] = size(signals)
gapwidthhist = zeros(numData, numSignals);
for column = 1 : numSignals
    thisSignal = signals(:, column); % Extract this column.
    % Find lengths of all NAN runs
    props = regionprops(isnan(thisSignal), 'Area');
    allLengths = [props.Area]
    edges = [1:max(allLengths), inf]
    hc = histcounts(allLengths, edges)
    % Load up gapwidthhist
    for k2 = 1 : length(hc)
        gapwidthhist(k2, column) = hc(k2);
    end
end
% What it is:
gapwidthhist'

但我主要在寻找没有任何内置 matlab 函数的纯 Matlab 代码（如 Image Processing Toolbox 中的“regionprops”）！！！

【问题讨论】：

为什么不循环遍历列的元素；如果元素是 NaN，则增加一个计数器；如果元素不是 NaN，则保存计数并重置计数器。
@beaker 你能详细说明你的方法吗？

标签： matlab signals histogram

【解决方案1】：

这是更简单的 Matlab 实现，但仍然不是最优的（+ 未矢量化）：

signals = [1.1 NaN NaN NaN  -1.4 NaN 8.3 NaN NaN NaN  NaN 1.5 NaN NaN; % signal No. 1
    NaN 2.2 NaN 4.9   NaN 8.2 NaN NaN NaN NaN  NaN 2.4 NaN NaN; % signal No. 2
        1 NaN NaN NaN NaN  NaN NaN NaN NaN  NaN NaN NaN NaN NaN]';  % signal No. 3
signals
[numData, numSignals] = size(signals);
gapwidthhist = zeros(numData, numSignals);
gaps = zeros(numData+1,numSignals);
auxnan = isnan(signals);
for i = 1:numSignals
    c = 0;
    for j = 1:numData
        if auxnan(j,i)
            c = c + 1;
        else
            gaps(j,i) = c;
            c = 0;
        end
    end
    gaps(numData+1,i) = c;
    gapwidthhist(:,i) = histcounts(gaps(:,i),1:numData+1);
end
gapwidthhist

感谢@breaker 的帮助。

知道如何优化（矢量化）此代码以更有效吗？

【讨论】：

这在什么方面不是最佳的？时间复杂度为 O(n)，这是可能的最小值，因为您必须检查数组的所有元素。您的实际数据的大小是多少？您的脚本在该数据上的执行时间是多少？
从良好的 Matlab 编程的角度来看，这段代码仍然不是最优的：1）gaps 数组使用的内存非常无效，2）每个信号列处理是完全独立的，因此外部循环与 parfor 的潜在并行化-loop 应该是可能的，但现在不行。 3) 应消除在外循环调用 histcounts 函数，并应直接在内循环计算 gapwidthhist(:,i)。 4) numData = 1e6 numSignals = 1e3.
5) 代码是串行的（仅运行一个核心）。执行时间：numData = 1e5 numSignals = 1e3：3.7sec，执行时间：numData = 1e6 numSignals = 1e2：6.7sec
执行时间很大程度上取决于 NaN 信号的数量！
啊，我错过了你还在使用histcounts。是的，不要那样做。将内部else 子句更改为gapwidthhist(c,i) = gapwidthhist(c,i) + 1; c = 0; 现在不需要gaps 数组。

【解决方案2】：

这是一个稍微矢量化的版本，可能会更快一些。我用的是 Octave，所以不知道 MATLAB 的 JIT 编译器在另一种方法中会优化多少内循环。


% Set up the data
signals = [1.1 NaN NaN NaN  -1.4 NaN 8.3 NaN NaN NaN  NaN 1.5 NaN NaN; % signal No. 1
    NaN 2.2 NaN 4.9   NaN 8.2 NaN NaN NaN NaN  NaN 2.4 NaN NaN; % signal No. 2
        1 NaN NaN NaN NaN  NaN NaN NaN NaN  NaN NaN NaN NaN NaN]';  % signal No. 3
signals
[numData, numSignals] = size(signals);
gapwidthhist = zeros(numData, numSignals);
gaps = zeros(numData+1,numSignals);
auxnan = ~isnan(signals);   % We want non-NaN values to be 1

for i = 1:numSignals
   difflist = diff(find([1; auxnan(:,i); 1])) - 1;   % get the gap lengths
   gapList = difflist(find(difflist));   % keep only the non-zero gaps
   for c = gapList.'                     % need row vector to loop over elements
      gapwidthhist(c,i) = gapwidthhist(c,i) + 1;   % each gap length increments the histogram
   end
end

gapwidthhist

这是程序流程：

首先，取反 auxnan 数组，使 NaN 为 0，非 NaN 为 1。
在外部循环中，用 1 在顶部和底部填充每列，以捕获信号开头和结尾处的 NaN 字符串。
使用 find 获取 1（非 NaN）元素的索引。
获取索引的diff。
diff 为 1 表示没有间隙，大于 1 的 diff 表示间隙长度加 1，因此从 diff 结果中减去 1。
使用find 的结果（索引）来获取非零元素的值。这些是间隙宽度。
现在循环遍历这些值并将结果累积到直方图中。您可以尝试用 accumarray 替换这个内部循环，看看是否可以加快速度。

【讨论】：

【解决方案3】：

可能是最终解决方案：

signals = [1.1 NaN NaN NaN  -1.4 NaN 8.3 NaN NaN NaN  NaN 1.5 NaN NaN; % signal No. 1
    NaN 2.2 NaN 4.9   NaN 8.2 NaN NaN NaN NaN  NaN 2.4 NaN NaN; % signal No. 2
        1 NaN NaN NaN NaN  NaN NaN NaN NaN  NaN NaN NaN NaN NaN]';  % signal No. 3

signals

[numData, numSignals] = size(signals);
gapwidthhist = zeros(numData, numSignals);
auxnan = isnan(signals);
for i = 1:numSignals
    c = 0;
    for j = 1:numData
        if auxnan(j,i)
            c = c + 1;
        else
            if c > 0
                gapwidthhist(c,i) = gapwidthhist(c,i) + 1;
                c = 0;
            end
        end
    end
    if c > 0
        gapwidthhist(c,i) = gapwidthhist(c,i) + 1;
    end
end

gapwidthhist

悬而未决的问题：如何修改外部for-loop应该可以使用parfor-loop的代码？

【讨论】：