understanding a newly proposed STA/LTA algorithm
36 views (last 30 days)
Show older comments
given below is the modified STA/LTA equation.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1271470/image.png)
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 second\n', 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 second\n', 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 region
9 Comments
Jorge Luis Paredes Estacio
on 30 Jul 2023
Hi, @PARVATHY NAIRParbathy. I am reading that paper about the STA/LTA alborithm. I was wondering of you managed to finish that code as I am working on detecting P- and S-wave from a great bunch of records. I've got another algorithm but I would like to compare it with the one you may finish. My email is zs19084@bristol.ac.uk or jlparedese@gmail.com. Thank you.
Answers (0)
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!