understanding a newly proposed STA/LTA algorithm

Question

PARVATHY NAIR 2023-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1898725-understanding-a-newly-proposed-sta-lta-algorithm

评论： Kardes ASLAN 2025-2-24

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 
Error using importdata
Unable to open file.
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

12 个评论
显示 10更早的评论隐藏 10更早的评论

PARVATHY NAIR 2023-1-23

在 MATLAB Online 中打开

Framework_for_Automated_Earthquake_Event_Detection_Based_on_Denoising_by_Adaptive_Filter.pdf

thankyou @Mathieu NOE

its the same paper.I will attach the same again.and i will resend the code again

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
%% Removing the delay introduced by lowpassiir filter
Acc_EW_dlycompensated = filtfilt(digfilt,Acc_EW);
Acc_NS_dlycompensated = filtfilt(digfilt,Acc_NS);
Acc_ver_dlycompensated = filtfilt(digfilt,Acc_ver);
 %% Finding the SNR
%  noise_ratio = snr(fildat) %returns the SNR in dBc of a real sinusoidal input signal, x, sampled at a rate fs. 
%  % The computation excludes the power contained in the lowest n harmonics, including the fundamental. 
%  % The default value of fs is 1. The default value of n is 6.
%  %fildat( .*) = 10^(3/20); %To increase the power of x by 3 dB:
%% plotting frequency spectrum
% fftsig = fft(fildat); %Taking fourier transform
% fftSig = fftshift(fftsig); %apply fftshift to put it in the form we are used to
% f = [-nt/2:nt/2-1]/nt; %Next, calculate the frequency axis, which is defined by the sampling rate
% figure;
% plot(f, abs(fftSig));
% title('magnitude FFT of sine');
% xlabel('Frequency (Hz)');
% ylabel('magnitude');
%% 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
%% Distance of earthquake from the seismometer
dist = (ts-tp)*8;
display(sprintf('Earthquake is estimated to be %f kilometers from the seismometer',dist))
%% Acceleration Plot
figure(1);
subplot(3,1,1)
plot(time,fildat,[tp tp],ylim,'r','LineWidth',1)
%plot(time,fildat)
title('Acceleration Data');
xlabel('Time (Sec)');
ylabel('Acceleration (cm/sec^2)');
grid on 
grid minor
%% Velocity Plot
subplot(3,1,2)
plot(time,vel)
title('Velocity Data');
xlabel('Time (Sec)');
ylabel('Velocity (cm/sec)');
grid on 
grid minor
%% Displacement Plot
subplot(3,1,3)
plot(time,dis)
title('Displacement Data');
xlabel('Time (Sec)');
ylabel('Displacement (cm)');
grid on 
grid minor
%% Plotting Spectogram of Original Signal and detecting the P-wave first arrival
figure(2)
box on
hold on
subplot(3,1,1)
plot(time,fildat,[tp tp],ylim,'r','LineWidth',2)
hold on
plot(time,fildat,[tp1 tp1],ylim,'g','LineWidth',2)
%line([tp1 tp1],[0,100],'Color','green','LineWidth',2);
%title('Acceleration Data');
xlabel('Time (Sec)');
ylabel('Acceleration (cm/sec^2)');
grid on 
grid minor
axis tight
box on
s = spectrogram(abs(fildat),256,250,256,200,'yaxis');
subplot(3,1,2)
%title('Spectrogram of Acceleration Data');
spectrogram(abs(fildat),256,250,256,200,'yaxis')
tp_in_min = tp/60;
tp_in_min1 = tp1/60;
%line([tp_in_min tp_in_min],[0,100],'Color','red','LineWidth',2);
line([tp_in_min1 tp_in_min1],[0,100],'Color','green','LineWidth',2);
grid on 
grid minor
axis tight;
box on
subplot(3,1,3)
thresh_spec = spectrogram(abs(fildat),256,250,256,200,'MinThreshold',-50,'yaxis');
thresh_spec1 = abs(thresh_spec);
%title('Spectrogram of Acceleration Data with -50dB Threshold');
spectrogram(abs(fildat),256,250,256,200,'MinThreshold',-50,'yaxis')
tp_in_min = tp/60;
tp_in_min1 = tp1/60;
%line([tp_in_min tp_in_min],[0,100],'Color','red','LineWidth',2);
line([tp_in_min1 tp_in_min1],[0,100],'Color','green','LineWidth',2);
grid on 
grid minor
axis tight;
%print('-dpdf','-r600','Temp.pdf')
%% Vizualizing the displacement of the seismometer
% Using an integrator to compute velocity from acceleration data.
% Integrator is filter with TF = 1/(1-Z^-1)
Nr = 1;
Dr = [1,-1];
% Detrending acceleration data it to remove drift,Integrating and dividing with sampling frequency to get velocity data
Ver_acc_dt = detrend(Acc_ver_dlycompensated);
NS_acc_dt = detrend(Acc_NS_dlycompensated);
EW_acc_dt = detrend(Acc_EW_dlycompensated);
vel_ver_nodrift = filter(Nr,Dr,Ver_acc_dt)/Fs; 
vel_NS_nodrift = filter(Nr,Dr,NS_acc_dt)/Fs;
vel_EW_nodrift = filter(Nr,Dr,EW_acc_dt)/Fs;
% Integrating velocity data and dividing with sampling frequency to get displacement data
dis_ver = filter(Nr,Dr,vel_ver_nodrift)/Fs;A
dis_NS = filter(Nr,Dr,vel_NS_nodrift)/Fs;
dis_EW = filter(Nr,Dr,vel_EW_nodrift)/Fs;
% Plot the positions using 3D plot to visualize seismometer movements
figure(3);
plot3(dis_NS,dis_EW,dis_ver);
grid on; view([-45,30]);
xlabel('N-S Direction in cm');
ylabel('E-W Direction in cm');
zlabel('Vertical Direction in cm');
title('Displacement of the seismometer in 3D');
set(gcf,'Name','Seismometer Trajectory');
set(gcf,'Units','Normalized');

汉羽张 2025-1-17

Please share me too. Email: 824284715@qq.com

Kardes ASLAN 2025-2-24

Dear Nair

I observed your matlab code. In my opinion, there are obvious errors in your code. First of all, equation 41 and equation 43 in the paper "Framework for Automated Earthquake Event Detection Based On Denoising by Adaptive Filter" have nothing to do with each other. These equations are 2 different types of calculating STA. However, you have related these equations to each other in the code you wrote. In the code you wrote, sta depends on staz. But it shouldn't be like that.

Similar situation is valid equations 42 and 44.

Furthermore, arrival time of p wave (tp) emerges greater than arrival time of s wave (ts).

Best regards...

请先登录，再进行评论。

请先登录，再回答此问题。

understanding a newly proposed STA/LTA algorithm

12 个评论
显示 10更早的评论隐藏 10更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

understanding a newly proposed STA/LTA algorithm

12 个评论 显示 10更早的评论隐藏 10更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

12 个评论
显示 10更早的评论隐藏 10更早的评论