具有两个条件的移动窗口-Matlab

Question

我有一个时间序列数据，其值的范围在1到100之间。我需要使用移动窗口（持续时间为10秒），并标记/突出显示阈值范围为91到200的数据段。 97至少10秒钟或更长时间。该数据可在以下链接中找到

https://www.dropbox.com/s/fx7z9qzg8gxb4x3/data.mat?dl=0

在这方面的任何帮助都将受到高度赞赏。

我使用了以下代码：

`x= load('data.mat');
time = 1:length(x);
window = 1000; % 10 second window - data was sampled at 100Hz.
[TF,L,U,C] = isoutlier(x); % I tried to find outlier
figure;
plot(time,x,time(TF),x(TF),'*')
legend('Data','Outlier')
figure;
plot(time,x,time(TF),x(TF),'x',time,93*ones(1,length(x)),time,
97*ones(1,length(x)))`

我得到下图[https://i.stack.imgur.com/fZa9m.jpg]但不确定。如何将其用作阈值窗口。

提前感谢

Answer 1

可能有一种方法可以使用实际的滑动窗口来执行此操作，但是由于我不得不以不同的方式对您执行类似的操作，因此我将重用我几年前完成的功能，该功能非常适合解决此问题另一种方式。

基本过程是：

• 查找所有验证条件的索引（x> 91＆x <97）。
• 查找在上面的​​步骤中找到的每个连续组/间隔的开始和停止索引。
• 每个间隔的计算长度（上面给出了开始和停止的难度）
• 丢弃间隔太短

到现在，您已经有了一个间隔列表（开始和停止索引），可以验证所有条件（信号电平和持续时间）。这样，您就可以：

• 重建有效索引的单个向量
• 创建仅包含标记数据（验证所有条件的x数据）的辅助向量。
• 显示：-）

在代码中看起来像这样：

%% Your inputs
load('data.mat');
time = 1:length(x);

%% Setup
winSize = 1000 ;    % Fs * Duration = 100 * 10 = 1000 points
lvlmin = 91 ;       % Minimum level to flag
lvlmax = 97 ;       % Maximum level to flag

%% Find which interval to flag
% find all the indices where the condition is true
idx = ( x>lvlmin ) & ( x<lvlmax ) ;
% get the start and stop index of each group of consecutive indices
itBounds = get_interval_boundaries( idx ) ;
% get the length of each interval/group
itLenght = diff(itBounds,1,2)+1 ; 
% only consider intervals >= winSize
it2flag = itLenght >= winSize ; 
nint = sum(it2flag) ;       % found 241 valid intervals out of 596.

%% Clear [idx] of the short intervals
itbad = itBounds( ~it2flag , : ) ; % keep only the intervals to discard
for k=1:size(itbad,1)
    idx(itbad(k,1):itbad(k,2)) = false ;
end

%% Display
flaggedTime = time(idx) ;
flaggedData = x(idx) ;

figure
plot(time,x)
hold on
plot(flaggedTime,flaggedData,'.')

lx = [time(1) time(end)] ;
plot( lx , [lvlmin lvlmin], '-.k')
plot( lx , [lvlmax lvlmax], '-.k')

%% OR, alternatively, keep vectors the same lenght by adding NaNs
flaggedData = x ;
flaggedData(~idx) = NaN ;

figure
plot(time,x)
hold on
plot(time,flaggedData)

以及有关如何标记数据的预览：

---

您需要get_interval_boundaries.m的代码。我可以只用较少的代码编写所需的功能，但是既然可以使用并且可以正常工作，那么就无需重新发明轮子了：

function itbound = get_interval_boundaries(vec)
% function itbound = get_interval_boundaries(vec)
%
% This function takes a vector of index as input (or a logical index array)
% It returns a nx2 table containing on each line the first and last index
% of consecutive intervals of indexes.
% ex:
% A = [1 2 3 5 7 8 9] ;
% [itbound] = get_interval_boundaries(A)
% itbound =
%      1     3
%      5     5
%      7     9
% 
% 09-Oct-2011 - Hoki:  creation
% 15-Sep-2012 - Hoki: Corrected last index special case
%                       (return last idx instead of 0)
% 01-Sep-2014 - Hoki: Corrected first index special case
%                       (return [1 1] in case vec is a scalar)

%% Check vec type (logical array or direct indexing)

% Return empty vector if input is empty
itbound = [] ;
if isempty(vec)
    return
end

% Check the type of input vector
if islogical(vec)
    idxDirect = find(vec) ;

elseif isnumeric(vec)
    idxDirect = vec ;
else
    errordlg('bad type for ''vec''. Variable should be numeric or logical',mfilename,'modal')
    return
end

%% Detect intervals
Npts = length(idxDirect) ;

% return [] in case vec is all [0]
if Npts == 0 ; return ; end

itbound(1,1) = idxDirect(1) ;
% return [x x] in case vec is a scalar value [x]
if Npts == 1
    itbound(1,2) = idxDirect(1) ;
    return
end

j=1 ;
for k = 2:Npts
    if idxDirect(k)==idxDirect(k-1)+1
        if k~=Npts
            % Cycle the loop
            continue
        else
            % Last point: Assign closing boundary of last interval
            itbound(j,2) = idxDirect(k) ;
        end
    else
        % Assign closing boundary of current interval
        itbound(j,2) = idxDirect(k-1) ;
        % Assign opening boundary of next interval
        j = j + 1 ;
        itbound(j,1) = idxDirect(k) ;
        % If we're on the very last index, close the interval.
        if k==Npts
            itbound(j,2) = idxDirect(k) ;
        end
    end
end