У меня есть большие наборы трехмерных данных, состоящие из одномерных сигналов, полученных в двумерном пространстве. Первым этапом обработки этих данных является установление порога всех сигналов для определения прихода высокоамплитудного импульса. Этот импульс присутствует во всех сигналах и поступает в разное время. После установки порогового значения трехмерный набор данных должен быть переупорядочен таким образом, чтобы каждый сигнал начинался с момента поступления импульса, а то, что было до этого, выбрасывалось (конец сигналов не имеет значения, так как сейчас я объединяю нули до конца всех сигналы, поэтому данные остаются одинакового размера).
Теперь я реализовал это следующим образом:
Сначала я начну с вычисления номера выборки первой выборки, превышающей порог в все сигналы
M = randn(1000,500,500); % example matrix of realistic size
threshold = 0.25*max(M(:,1,1)); % 25% of the maximum in the first signal as threshold
[~,index] = max(M>threshold); % indices of first sample exceeding threshold in all signals
Далее я хочу, чтобы все сигналы были сдвинуты так, чтобы все они начинались с импульса. На данный момент я реализовал это следующим образом:
outM = zeros(size(M)); % preallocation for speed
for i = 1:size(M,2)
for j = 1:size(M,3)
outM(1:size(M,1)+1-index(1,i,j),i,j) = M(index(1,i,j):end,i,j);
end
end
Это работает нормально, и я знаю, что циклы for больше не такие медленные, но это легко занимает несколько секунд для наборов данных на моей машине. Одна итерация for-l oop занимает около 0,05-0,1 se c, что мне кажется медленным для простого копирования вектора, содержащего 500-2000 двойных значений.
Поэтому я рассмотрел лучший способ справиться с этим, но пока я не нашел ничего лучшего. Я пробовал несколько вещей: 3D-маски, линейное индексирование и параллельные циклы (parfor).
для 3D-масок, я проверил, возможны ли какие-либо улучшения. Поэтому сначала я создаю логическую маску, а затем сравниваю скорость индексации / копирования логической маски с двойной вложенностью для l oop.
%% set up for logical mask copying
AA = logical(ones(500,1)); % only copy the first 500 values after the threshold value
Mask = logical(zeros(size(M)));
Jepla = zeros(500,size(M,2),size(M,3));
for i = 1:size(M,2)
for j = 1:size(M,3)
Mask(index(1,i,j):index(1,i,j)+499,i,j) = AA;
end
end
%% speed comparison
tic
Jepla = M(Mask);
toc
tic
for i = 1:size(M,2)
for j = 1:size(M,3)
outM(1:size(M,1)+1-index(1,i,j),i,j) = M(index(1,i,j):end,i,j);
end
end
toc
Формула-l oop каждый раз быстрее даже если скопировано больше.
Далее, линейное индексирование.
%% setup for linear index copying
%put all indices in 1 long column
LongIndex = reshape(index,numel(index),1);
% convert to linear indices and store in new variable
linearIndices = sub2ind(size(M),LongIndex,repmat(1:size(M,2),1,size(M,3))',repelem(1:size(M,3),size(M,2))');
% extend linear indices with those of all values to copy
k = zeros(numel(M),1);
count = 1;
for i = 1:numel(LongIndex)
values = linearIndices(i):size(M,1)*i;
k(count:count+length(values)-1) = values;
count = count + length(values);
end
k = k(1:count-1);
% get linear indices of locations in new matrix
l = zeros(length(k),1);
count = 1;
for i = 1:numel(LongIndex)
values = repelem(LongIndex(i)-1,size(M,1)-LongIndex(i)+1);
l(count:count+length(values)-1) = values;
count = count + length(values);
end
l = k-l;
% create new matrix
outM = zeros(size(M));
%% speed comparison
tic
outM(l) = M(k);
toc
tic
for i = 1:size(M,2)
for j = 1:size(M,3)
outM(1:size(M,1)+1-index(1,i,j),i,j) = M(index(1,i,j):end,i,j);
end
end
toc
Опять же, альтернативный подход, линейное индексирование, (намного) медленнее. После того, как это не удалось, я узнал о распараллеливании, и хотя это наверняка ускорило бы мой код. Прочитав некоторую документацию по parfor и немного опробовав ее, я изменил свой код следующим образом:
gcp;
outM = zeros(size(M));
inM = mat2cell(M,size(M,1),ones(size(M,2),1),size(M,3));
tic
parfor i = 1:500
for j = 1:500
outM(:,i,j) = [inM{i}(index(1,i,j):end,1,j);zeros(index(1,i,j)-1,1)];
end
end
end
toc
Я изменил его так, чтобы «outM» и «inM» оба были нарезанными переменными, как я читаю это лучше. Тем не менее, это очень медленно, намного медленнее, чем оригинал для l oop.
Так что теперь вопрос, должен ли я отказаться от попыток улучшить скорость этой операции? Или есть другой способ сделать это? Я много искал и пока не вижу, как это ускорить. Извините за длинный вопрос, но я хотел показать, что я пытался. Заранее спасибо!