我需要有关如何优化我的代码的建议。执行时间太长
I want advice about how to optimize my code. It takes too long for execution
我编写了一个 MATLAB 代码,用于从 SAC(地震)文件(通过另一个代码读取)中查找地震信号(例如 P 波)。这个算法叫做STA/LTA触发算法(其实对我的问题没那么重要)
重要的是其实这段代码运行良好,但是由于我的地震文件太大(1GB,两个月),执行了将近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;
如果您有 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
是 'only' 8 MB,则此方法需要不到 0.02 秒,而原始方法需要将近 6 秒!
但是,即使您没有 Matlab 2016a,考虑到您的循环,您仍然可以执行以下操作:
- 删除
else
语句 - sta(i)
在预分配中已经为零。
- 从
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
但不会那么快。
我编写了一个 MATLAB 代码,用于从 SAC(地震)文件(通过另一个代码读取)中查找地震信号(例如 P 波)。这个算法叫做STA/LTA触发算法(其实对我的问题没那么重要)
重要的是其实这段代码运行良好,但是由于我的地震文件太大(1GB,两个月),执行了将近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;
如果您有 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
是 'only' 8 MB,则此方法需要不到 0.02 秒,而原始方法需要将近 6 秒!
但是,即使您没有 Matlab 2016a,考虑到您的循环,您仍然可以执行以下操作:
- 删除
else
语句 -sta(i)
在预分配中已经为零。 - 从
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
但不会那么快。