Я хочу совет о том, как оптимизировать мой код. Это занимает слишком много времени для исполнения

Question

Я написал код MATLAB для поиска сейсмического сигнала (например, P-волны) из файла SAC(сейсмического) (который читается через другой код). Этот алгоритм называется алгоритмом запуска STA/LTA (на самом деле это не так важно для моего вопроса)

Важно то, что на самом деле этот код работает хорошо, но поскольку мой сейсмический файл слишком большой (1 ГБ, то есть за два месяца), для просмотра результата требуется почти 40 минут. Таким образом, я чувствую необходимость оптимизировать код.

Я слышал, что замена циклов расширенными функциями поможет, но, поскольку я новичок в MATLAB, я не могу понять, как это сделать, поскольку целью кода является сканирование каждого временного ряда. Кроме того, я слышал, что предварительное распределение может помочь, но у меня есть простая идея о том, как на самом деле это сделать.

Поскольку этот код посвящен сейсмологии, его может быть трудно понять, но мои заметки вверху могут помочь. Я надеюсь, что смогу получить полезный совет здесь. Ниже мой код.

function[pstime]=classic_LR(tseries,ltw,stw,thresh,dt)

% This is the code for "Classic LR" algorithm
% 'ns' is the number of measurement in STW-used to calculate STA 
% 'nl' is the number of measurement in LTW-used to calculate LTA
% 'dt' is the time gap between measurements i.e. 0.008s for HHZ and 0.02s for BHZ
% 'ltw' and 'stw' are long and short time windows respectively
% 'lta' and 'sta' are long and short time windows average respectively
% 'sra' is the ratio between 'sta' and 'lta' which will be each component
% for a vector containing the ratio for each measurement point 'i'
% Index 'i' denotes each measurement point and this will be converted to actual time

nl=fix(ltw/dt);
ns=fix(stw/dt);
nt=length(tseries);
aseries=abs(detrend(tseries));
sra=zeros(1,nt);

for i=1:nt-ns
    if  i>nl 
        lta=mean(aseries(i-nl:i));
        sta=mean(aseries(i:i+ns));
        sra(i)=sta/lta;
    else
        sra(i)=0;
    end
end

[k]=find(sra>thresh);
if ~isempty(k)
    pstime=k*dt;
else 
    pstime=0;
end

return;

Answer 1

Если у вас есть MATLAB 2016a или более поздняя версия, вы можете использовать movmean вместо вашего цикла (это означает, что вам также не нужно ничего предварительно выделять):

lta = movmean(aseries(1:nt-ns),nl+1,'Endpoints','discard');
sta = movmean(aseries(nl+1:end),ns+1,'Endpoints','discard');
sra = sta./lta;

Единственная разница здесь в том, что вы получите sra без начальных и конечных нулей. Скорее всего, это будет самый быстрый способ. Если, например, aseries Это "всего" 8 МБ, чем этот метод занимает менее 0,02 секунды, в то время как оригинальный метод занимает почти 6 секунд!

---

Однако, даже если у вас нет Matlab 2016a, учитывая ваш цикл, вы все равно можете сделать следующее:

1. Удалить else заявление - sta(i) уже ноль от предварительного выделения.
2. Начать цикл с nl+1 вместо проверки когда i больше, чем nl,

Итак, ваш новый цикл будет:

for i=nl+1:nt-ns
    lta = mean(aseries(i-nl:i));
    sta = mean(aseries(i:i+ns));
    sra(i)=sta/lta;
end

Но это не будет так быстро.