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

Question

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

我通过仅运行时间序列的 10000 个初始点来测试代码。目前，我运行前 10,000 个点大约需要 3 秒。我在一台相对较旧的 32 位戴尔笔记本电脑上运行，配备 2.30 GHz 处理器和 4 GB 内存。显然，如果我使用较新的计算机，我可以很快完成任务。例如，我的功能更强大的台式机在 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”命令可能有更优雅的方式。

欢迎任何想法。

干杯

Answer 1

Any suggestions for how to improve the speed is appreciated.

您可以通过不重复您的 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。