understanding a newly proposed STA/LTA algorithm

given below is the modified STA/LTA equation.

from my understanding i have formulated the code as follows.Kindly go through the code and correct me if my notion is wrong

clc
clear all
%% Data Inputs
Acc_EW = importdata("ADIB.HHE.dat")
Acc_NS = importdata("ADIB.HHN.dat");
Acc_ver = importdata("ADIB.HHZ.dat");
Fs = 100; %sampling frequency
%% Signal Pre-Processing
%Filter Design
digfilt = designfilt(‘lowpassiir’, ‘PassbandFrequency’, 20, ‘StopbandFrequency’, 25, ‘PassbandRipple’, 1, ‘StopbandAttenuation’, 60, ‘SampleRate’, 200);
% Filtering Data
Acc_EW_filt = filter(digfilt,Acc_EW);
Acc_NS_filt = filter(digfilt,Acc_NS);
Acc_ver_filt = filter(digfilt,Acc_ver);
Fhp = 0.8; % high pass filter cutofff frequency
[b1,a1] = butter(3,Fhp/Fs,’high’); %3rd order high pass butterworth filter
fildat = filter(b1,a1,Acc_ver); % filtered acceleration data
vel = cumtrapz(fildat)./Fs; % Integrating acceleration data for velocity
[b2,a2] = butter(3,Fhp/Fs,’high’); %3rd order high pass butterworth filter
fildat1 = filter(b2,a2,vel); % filtered velocity data
dis = cumtrapz(fildat1)./Fs; % Integrating velocity data for displacement
peakToPeakRange = max(fildat) – min(fildat);
dt = 1/Fs; %sampling time
nt = length(fildat); % length of the input signal
time = (1:nt).*dt; % time duration of the input signal

%% STA-LTA Algorithm gor P-Wave detection
stw = 0.2; %short time window length
ltw = 70; %long time window length
thresh = 3; % Threshold
thresh1 = 4;
%t = 1;
nl = fix(ltw / dt); %no. of data points in the long time window
ns = fix(stw / dt); %no. of data points in the short time window
nt = length(fildat);
sra = zeros(1, nt);

%%this where i have modified accroding to excerpt from the paper-‘Framework
%%for Automated Earthquake Event Detection Based On Denoising by Adaptive
%%Filter’
for k = nl+1:nt
staz(k,1) = (1/ns)* trapz(abs(fildat(k-ns:k)));
ltaz(k,1) = (1/nl)* trapz(abs(fildat(k-nl:k)));
sta(k,1) = (1/ns)* [(staz(k-1)*ns)-fildat(k-ns)-fildat(k)];
lta(k,1) = (1/nl)* [(ltaz(k-1)*nl)-fildat(k-nl)-fildat(k)];
end
for l = nl+1: nt
sra(l) = sta(l)/lta(l);
end
itm = find(sra > thresh);
if ~isempty(itm)
itmax=itm(1);
end
tp =itmax*dt; % P-wave arriving time
fprintf(‘P-Wave detection time for threshold 4 = %f secondn’, tp);

itm1 = find(sra > thresh1);
if ~isempty(itm1)
itmax1 = itm1(1);
end
tp1 = itmax1*dt; % P-wave arriving time
fprintf(‘P-Wave detection time for threshold 3 = %f secondn’, tp1);
%% S-wave arrival time
pkHts = 0.72; % 10 percent
[pk2,t22] = findpeaks(Acc_NS_dlycompensated,Fs,’MinPeakHeight’,pkHts*max(Acc_ver_dlycompensated),’Npeaks’,1);
[pk3,t33] = findpeaks(Acc_EW_dlycompensated,Fs,’MinPeakHeight’,pkHts*max(Acc_ver_dlycompensated),’Npeaks’,1);

display(sprintf(‘S-wave found on EW component at %f seconds and on NS componet at %f seconds,’, t33,t22));

if(t22<t33)
display(‘S-wave detected first on North-South component’);
else
display(‘S-wave detected first on East-West component’);
end
ts = min(t22,t33);
line([ts,ts],[min(get(gca,’Ylim’))],’linestyle’,’–‘,’linewidth’,2,’color’,’red’);

%% Tauc , Pd and Magnitude calculations
vel_sq = vel.^2;
dis_sq = dis.^2;
r1 = trapz(vel_sq((itmax):(itmax+300)));
r2 = trapz(dis_sq((itmax):(itmax+300)));
r = r1/r2;
tauc = 2*pi/sqrt(r);
pd = max(dis((itmax):(itmax+300)));
mag_tauc = (log(tauc) + 3.45)/0.47 %Coefficients varies from region to region
mag_pd = (0.873*((log(pd)+6.3)/0.513))+4.74 %Coefficients varies from region to regiongiven below is the modified STA/LTA equation.

from my understanding i have formulated the code as follows.Kindly go through the code and correct me if my notion is wrong

display(sprintf(‘S-wave found on EW component at %f seconds and on NS componet at %f seconds,’, t33,t22));

from my understanding i have formulated the code as follows.Kindly go through the code and correct me if my notion is wrong

display(sprintf(‘S-wave found on EW component at %f seconds and on NS componet at %f seconds,’, t33,t22));

Cart

Cart