MATLAB：寻找更快/更优雅的方法来计算非常长的时间序列的移动中值绝对偏差

Posted 2023-03-11

【中文标题】MATLAB：寻找更快/更优雅的方法来计算非常长的时间序列的移动中值绝对偏差

【英文标题】：MATLAB: Looking for faster/more elegant way to compute a moving median absolute deviation for very long time series

【发布时间】：2017-12-26 18:06:37

【问题描述】：

我有一个非常大的数组，有五个通道和大约 600 万个条目 (5 x 6000000)。我的目标是用 7 点窗口扫描阵列并移除“尖峰”，这些“尖峰”被定义为大于中值绝对偏差 (MAD) 的某个缩放量。

我只运行时间序列的 10000 个初始点来测试代码。目前，我运行前 10,000 个点大约需要 3 秒。我在一台相对较旧的 32 位戴尔笔记本电脑上运行，该笔记本电脑配备 2.30 GHz 处理器和 4 GB RAM。显然，如果我使用的是较新的计算机，我可以很快完成任务。例如，我功能更强大的桌面在 0.7 秒内完成相同的任务。但是，我需要在笔记本电脑上运行代码，并且每次需要运行代码时都不能等待 35 - 40 分钟。我正在寻求帮助，以寻找效率低下的地方以及可以使代码更快的地方。

下面是代码。任何有关如何提高速度的建议表示赞赏。我注意到“MAD”的计算是最耗时的（大约需要 1.9 秒，或超过总时间的一半）。

load('data.mat') % data (approx 5 channels x 6000000 data points (int32))
nscans = length(data); %number of data points in each channel

nwide = 7; %number of data points in the window

% Rejection parameters (not so important for the question)
iscale = 50; %scale factor for MAD
minmad = 2;
mincrit = [100 100 100 500 500];

nfixed=zeros(1,5);
L = floor(nwide/2); %half of window (odd window length only)

%Padding for beginning and end of data
data = [repmat(data(:,1),[1 L]) data repmat(data(:,end),[1 L])];

nfixed = zeros(1,5); %initialize counter
tic 
for n = L+1:10000
    idata = data(:,n-L:n+L)'; % temporary window

    % compute median of window
    med=median(idata);

    %compute median absolute deviation (MAD)
    % Note: mad = median(abs(X - median(X)))
    mad = median(abs(idata-repmat(median(idata),[nwide 1])));
    mad = max([mad;minmad*ones(1,5)]); %minmad threshold added

    %compute rejection threshold
    icrit=max([iscale*mad;mincrit]);

    for i = 1:5 %loop over channels
        if abs(data(i,n)-med(i)) > icrit(i) %if threshold is exceeded
            data(i,n)=med(i); %then replace with median value
            nfixed(i)=nfixed(i)+1; %count number of replacements
        end
    end

end
toc

data = data(:,L+1:end-L)'; %remove padding

我觉得可能有一种更优雅的方式来执行“repmat”命令。

感谢任何想法。

干杯

【问题讨论】：

当你移除一个尖峰时，这会影响未来窗口中的计算吗？它看起来如此形成您的代码。既然如此，就很难避免循环了

我必须考虑更多，但第一件事是你可以做array-scalar 而不必重新设置它。这可能会为您节省一些时间。

【参考方案1】：

感谢任何有关如何提高速度的建议。

您可以通过不再重复您的median(idata) 调用来稍微收紧您的代码。

改变这个：

mad = median(abs(idata-repmat(median(idata),[nwide 1])));

到这里：

mad = median(abs(idata-repmat(med,[nwide 1])));

或者，您可能会从 MATLAB's mad 函数中获得更多里程，它出现在 2006 年之前。不过，您需要更改变量名称。

例如，您可以从此更改您的代码：

mad = median(abs(idata-repmat(median(idata),[nwide 1])));
mad = max([mad;minmad*ones(1,5)]); %minmad threshold added

到

madV = max(mad(idata);[2 2 2 2 2]);

我只是将 2 的向量放在那里，因为代码中没有任何内容显示 minmad 正在更新。

【讨论】：

matlab 的mad 命令不太可能有帮助。我查看了源代码，这正是他的实现，开销更大。可能需要更长的时间，除非 2017b 之前的版本将其混合，并且他们出于某种原因将其改回。

@Durkee 我明白你对mad() 的意思。希望调用median(idata) 一次而不是两次仍然会有所帮助。

【参考方案2】：

你能用movmedian得到移动中位数吗？

med = movmedian(data,7,2,'Endpoints','discard');

尽可能避免使用 for 循环。

例如：

data = randn([5,1E6]);tic;
med = movmedian(data,7,2); %moving median
dev = abs(data-med); %deviation from median
thres = median(dev,2); %threshold
rep = dev>thres; %points to replace
data(rep) = med(rep); %replace data with median
toc
>>Elapsed time is 0.285828 seconds.
memory
>>Memory used by MATLAB:      1629 MB (1.708e+09 bytes)