了解一个新提出的STA / LTA算法

45视图(30天)

显示旧的评论

PARVATHY NAIR 2023年1月23日

0
链接

这个问题直接联系

//www.tianjin-qmedu.com/matlabcentral/answers/1898725-understanding-a-newly-proposed-sta-lta-algorithm

评论道: PARVATHY NAIR2023年1月27日

下面是修改后的STA / LTA方程。

从我的理解我制定如下的代码。请通过代码和纠正我如果我的想法是错误的

                          clc
                         
                          清晰的所有
                         
                          % %数据输入
                         
                             错误使用importdata
无法打开文件。

                          Acc_EW = importdata (“ADIB.HHE.dat”)
                         
                          Acc_NS = importdata (“ADIB.HHN.dat”);
                         
                          Acc_ver = importdata (“ADIB.HHZ.dat”);
                         
                          Fs = 100;%采样频率
                         
                          % %信号预处理
                         
                          %滤波器设计
                         
                          digfilt = designfilt (“lowpassiir”,“PassbandFrequency”,20岁,“StopbandFrequency”25岁的“PassbandRipple”, 1“StopbandAttenuation”现年60岁的“SampleRate”,200);
                         
                          %过滤数据
                         
                          Acc_EW_filt =过滤器(digfilt Acc_EW);
                         
                          Acc_NS_filt =过滤器(digfilt Acc_NS);
                         
                          Acc_ver_filt =过滤器(digfilt Acc_ver);
                         
                          佛罗里达州公路巡警部门的= 0.8;%高通滤波器cutofff频率
                         
                          (b1, a1) =黄油(3,当天/ Fs,“高”);%三阶巴特沃斯高通滤波器
                         
                          fildat =过滤器(b1、a1、Acc_ver);%加速度数据过滤
                         
                          韦尔= cumtrapz (fildat)。/ Fs;%集成加速度数据的速度
                         
                          (b2 a2) =黄油(3,当天/ Fs,“高”);%三阶巴特沃斯高通滤波器
                         
                          fildat1 =过滤器(b2, a2,韦尔);%过滤速度数据
                         
                          说= cumtrapz (fildat1)。/ Fs;%整合速度位移数据
                         
                          peakToPeakRange = max (fildat) - min (fildat);
                         
                          dt = 1 / f;%采样时间
                         
                          元=长度(fildat);%输入信号的长度
                         
                          时间= (1:nt)。* dt;%的时间输入信号
                         
                          % % STA-LTA算法气油比纵波检测
                         
                          stw = 0.2;%短时间窗长度
                         
                          ltw = 70;%长时间窗口长度
                         
                          打= 3;%的阈值
                         
                          thresh1 = 4;
                         
                          % t = 1;
                         
                          问=修复(ltw / dt);%。长时间数据点的窗口
                         
                          ns =修复(stw / dt);%。短时间内数据点的窗口
                         
                          元=长度(fildat);
                         
                          sra = 0(1元);
                         
                          % %这地方我已经修改根据摘录论文——“框架
                         
                          % %为自动化的地震事件检测基于自适应去噪
                         
                          % %过滤器的
                         
                          为本地语言+ k = 1: nt
                         
                          staz (k, 1) = (1 / ns) * trapz (abs (fildat (k - n, k)));
                         
                          ltaz (k, 1) = (1 / nl) * trapz (abs (fildat (k-nl: k)));
                         
                          sta (k, 1) = (1 / ns) * ((staz (k - 1) * ns) -fildat (k - n) -fildat (k)];
                         
                          英国网球协会(k, 1) = (1 / nl) * ((ltaz (k - 1) * nl) -fildat (k-nl) -fildat (k)];
                         
                          结束
                         
                          为l =问+ 1:nt
                         
                          sra (l) = sta(左)/ lta(左);
                         
                          结束
                         
                          itm =找到(sra >打);
                         
                          如果~ isempty (itm)
                         
                          itmax = itm (1);
                         
                          结束
                         
                          tp = itmax * dt;% p波到达时间
                         
                          流(“纵波检测时间阈值4 = % f第二\ n '、tp);
                         
                          itm1 =找到(sra > thresh1);
                         
                          如果~ isempty (itm1)
                         
                          itmax1 = itm1 (1);
                         
                          结束
                         
                          tp1 = itmax1 * dt;% p波到达时间
                         
                          流(“纵波检测时间阈值3 = % f第二\ n '、tp1);
                         
                          % %横波到达时间
                         
                          pkHts = 0.72;% 10%
                         
                          [pk2, t22] = findpeaks (Acc_NS_dlycompensated Fs,“MinPeakHeight”,pkHts *马克斯(Acc_ver_dlycompensated),“Npeaks”1);
                         
                          [pk3, t33] = findpeaks (Acc_EW_dlycompensated Fs,“MinPeakHeight”,pkHts *马克斯(Acc_ver_dlycompensated),“Npeaks”1);
                         
                          显示器(sprintf (横波在电子战的组件在% f秒和NS % f秒组合,”、t33 t22));
                         
                          如果(t22 < t33)
                         
                          显示器(“南北组件横波检测到第一”);
                         
                          其他的
                         
                          显示器(“东西方组件横波检测到第一”);
                         
                          结束
                         
                          ts = min (t22 t33);
                         
                          线([ts, ts], [min (get (gca),“Ylim”))),“线型”,“——”,“线宽”2,“颜色”,“红色”);
                         
                          % % Tauc, Pd和大小的计算
                         
                          vel_sq =韦尔。^ 2;
                         
                          dis_sq =说。^ 2;
                         
                          r1 = trapz (vel_sq ((itmax): (itmax + 300)));
                         
                          r2 = trapz (dis_sq ((itmax): (itmax + 300)));
                         
                          r = r1 / r2;
                         
                          tauc = 2π/√(r);
                         
                          pd = max (dis ((itmax): (itmax + 300)));
                         
                          mag_tauc =(日志(tauc) + 3.45) / 0.47%系数因地区而异
                         
                          mag_pd =(0.873 *((日志(pd) + 6.3) / 0.513)) + 4.74%系数因地区而异

7评论
显示6年长的评论隐藏6年长的评论

PARVATHY NAIR 2023年1月23日

Framework_for_Automated_Earthquake_Event_Detection_Based_on_Denoising_by_Adaptive_Filter.pdf

谢谢 @Mathieu诺

同一篇论文。我再次附上相同的。和i will resend the code again

                                clc
                               
                                清晰的所有
                               
                                % %数据输入
                               
                                Acc_EW = importdata (“ADIB.HHE.dat”)
                               
                                Acc_NS = importdata (“ADIB.HHN.dat”);
                               
                                Acc_ver = importdata (“ADIB.HHZ.dat”);
                               
                                Fs = 100;%采样频率
                               
                                % %信号预处理
                               
                                %滤波器设计
                               
                                digfilt = designfilt (“lowpassiir”,“PassbandFrequency”,20岁,“StopbandFrequency”25岁的“PassbandRipple”, 1“StopbandAttenuation”现年60岁的“SampleRate”,200);
                               
                                %过滤数据
                               
                                Acc_EW_filt =过滤器(digfilt Acc_EW);
                               
                                Acc_NS_filt =过滤器(digfilt Acc_NS);
                               
                                Acc_ver_filt =过滤器(digfilt Acc_ver);
                               
                                佛罗里达州公路巡警部门的= 0.8;%高通滤波器cutofff频率
                               
                                (b1, a1) =黄油(3,当天/ Fs,“高”);%三阶巴特沃斯高通滤波器
                               
                                fildat =过滤器(b1、a1、Acc_ver);%加速度数据过滤
                               
                                韦尔= cumtrapz (fildat)。/ Fs;%集成加速度数据的速度
                               
                                (b2 a2) =黄油(3,当天/ Fs,“高”);%三阶巴特沃斯高通滤波器
                               
                                fildat1 =过滤器(b2, a2,韦尔);%过滤速度数据
                               
                                说= cumtrapz (fildat1)。/ Fs;%整合速度位移数据
                               
                                peakToPeakRange = max (fildat) - min (fildat);
                               
                                dt = 1 / f;%采样时间
                               
                                元=长度(fildat);%输入信号的长度
                               
                                时间= (1:nt)。* dt;%的时间输入信号
                               
                                % % lowpassiir滤波器消除延迟
                               
                                Acc_EW_dlycompensated = filtfilt (digfilt Acc_EW);
                               
                                Acc_NS_dlycompensated = filtfilt (digfilt Acc_NS);
                               
                                Acc_ver_dlycompensated = filtfilt (digfilt Acc_ver);
                               
                                % %的信噪比
                               
                                % noise_ratio =信噪比(fildat) %返回dBc真正的正弦输入信号的信噪比,x,采样速率fs。
                               
                                % %计算不包括电源包含在最低的n次谐波,包括最基本的。
                               
                                % % fs的默认值是1。n的默认值是6。
                               
                                % % fildat (. *) = 10 ^ (3/20);%增加x 3 dB的力量:
                               
                                % %策划频谱
                               
                                % fftsig = fft (fildat);%进行傅里叶变换
                               
                                % fftSig = fftshift (fftSig);% fftshift适用于我们习惯于把它放在形式
                               
                                % f = (- nt / 2:元/ 2 - 1)/ nt;%,计算频率轴,定义的采样率
                               
                                %图;
                               
                                %的阴谋(f, abs (fftSig));
                               
                                %的标题(“正弦的FFT级”);
                               
                                %包含(频率(赫兹));
                               
                                % ylabel(级);
                               
                                % % STA-LTA算法气油比纵波检测
                               
                                stw = 0.2;%短时间窗长度
                               
                                ltw = 70;%长时间窗口长度
                               
                                打= 3;%的阈值
                               
                                thresh1 = 4;
                               
                                % t = 1;
                               
                                问=修复(ltw / dt);%。长时间数据点的窗口
                               
                                ns =修复(stw / dt);%。短时间内数据点的窗口
                               
                                元=长度(fildat);
                               
                                sra = 0(1元);
                               
                                % %这地方我已经修改根据摘录论文——“框架
                               
                                % %为自动化的地震事件检测基于自适应去噪
                               
                                % %过滤器的
                               
                                为本地语言+ k = 1: nt
                               
                                staz (k, 1) = (1 / ns) * trapz (abs (fildat (k - n, k)));
                               
                                ltaz (k, 1) = (1 / nl) * trapz (abs (fildat (k-nl: k)));
                               
                                sta (k, 1) = (1 / ns) * ((staz (k - 1) * ns) -fildat (k - n) -fildat (k)];
                               
                                英国网球协会(k, 1) = (1 / nl) * ((ltaz (k - 1) * nl) -fildat (k-nl) -fildat (k)];
                               
                                结束
                               
                                为l =问+ 1:nt
                               
                                sra (l) = sta(左)/ lta(左);
                               
                                结束
                               
                                itm =找到(sra >打);
                               
                                如果~ isempty (itm)
                               
                                itmax = itm (1);
                               
                                结束
                               
                                tp = itmax * dt;% p波到达时间
                               
                                流(“纵波检测时间阈值4 = % f第二\ n '、tp);
                               
                                itm1 =找到(sra > thresh1);
                               
                                如果~ isempty (itm1)
                               
                                itmax1 = itm1 (1);
                               
                                结束
                               
                                tp1 = itmax1 * dt;% p波到达时间
                               
                                流(“纵波检测时间阈值3 = % f第二\ n '、tp1);
                               
                                % %横波到达时间
                               
                                pkHts = 0.72;% 10%
                               
                                [pk2, t22] = findpeaks (Acc_NS_dlycompensated Fs,“MinPeakHeight”,pkHts *马克斯(Acc_ver_dlycompensated),“Npeaks”1);
                               
                                [pk3, t33] = findpeaks (Acc_EW_dlycompensated Fs,“MinPeakHeight”,pkHts *马克斯(Acc_ver_dlycompensated),“Npeaks”1);
                               
                                显示器(sprintf (横波在电子战的组件在% f秒和NS % f秒组合,”、t33 t22));
                               
                                如果(t22 < t33)
                               
                                显示器(“南北组件横波检测到第一”);
                               
                                其他的
                               
                                显示器(“东西方组件横波检测到第一”);
                               
                                结束
                               
                                ts = min (t22 t33);
                               
                                线([ts, ts], [min (get (gca),“Ylim”))),“线型”,“——”,“线宽”2,“颜色”,“红色”);
                               
                                % % Tauc, Pd和大小的计算
                               
                                vel_sq =韦尔。^ 2;
                               
                                dis_sq =说。^ 2;
                               
                                r1 = trapz (vel_sq ((itmax): (itmax + 300)));
                               
                                r2 = trapz (dis_sq ((itmax): (itmax + 300)));
                               
                                r = r1 / r2;
                               
                                tauc = 2π/√(r);
                               
                                pd = max (dis ((itmax): (itmax + 300)));
                               
                                mag_tauc =(日志(tauc) + 3.45) / 0.47%系数因地区而异
                               
                                mag_pd =(0.873 *((日志(pd) + 6.3) / 0.513)) + 4.74%系数因地区而异
                               
                                % %的地震检波器的距离
                               
                                dist = (ts-tp) * 8;
                               
                                显示器(sprintf (“地震估计% f公里从地震仪的dist))
                               
                                % %加速度图
                               
                                图(1);
                               
                                次要情节(1,1)
                               
                                情节(时间、fildat tp tp, ylim,“r”,“线宽”,1)
                               
                                %的阴谋(时间、fildat)
                               
                                标题(“加速度数据”);
                               
                                包含(的时间(秒));
                               
                                ylabel (“加速度(厘米/秒^ 2)”);
                               
                                网格在
                               
                                网格小
                               
                                % %速度图
                               
                                次要情节(1、2)
                               
                                情节(时间,韦尔)
                               
                                标题(“速度数据”);
                               
                                包含(的时间(秒));
                               
                                ylabel (“速度(厘米/秒)”);
                               
                                网格在
                               
                                网格小
                               
                                % %位移图
                               
                                次要情节(3、1,3)
                               
                                情节(时间,说)
                               
                                标题(“位移数据”);
                               
                                包含(的时间(秒));
                               
                                ylabel (“位移(cm)”);
                               
                                网格在
                               
                                网格小
                               
                                % %策划Spectogram原始信号和检测p波初至
                               
                                图(2)
                               
                                盒子在
                               
                                持有在
                               
                                次要情节(1,1)
                               
                                情节(时间、fildat tp tp, ylim,“r”,“线宽”,2)
                               
                                持有在
                               
                                情节(时间、fildat [tp1 tp1], ylim,‘g’,“线宽”,2)
                               
                                %线([tp1 tp1],[0100],“颜色”、“绿色”、“线宽”,2);
                               
                                %的标题(“加速度数据”);
                               
                                包含(的时间(秒));
                               
                                ylabel (“加速度(厘米/秒^ 2)”);
                               
                                网格在
                               
                                网格小
                               
                                轴紧
                               
                                盒子在
                               
                                s =光谱图(abs (fildat), 256250256200年,“桠溪”);
                               
                                次要情节(1、2)
                               
                                %的标题(“加速度数据的光谱图”);
                               
                                光谱图(abs (fildat), 256250256200年,“桠溪”)
                               
                                tp_in_min = tp / 60;
                               
                                tp_in_min1 = tp1/60;
                               
                                %线([tp_in_min tp_in_min],[0100],“颜色”,“红”,“线宽”,2);
                               
                                线([tp_in_min1 tp_in_min1], [0100],“颜色”,“绿色”,“线宽”2);
                               
                                网格在
                               
                                网格小
                               
                                轴紧;
                               
                                盒子在
                               
                                次要情节(3、1,3)
                               
                                thresh_spec =光谱图(abs (fildat), 256250256200年,“MinThreshold”,-50,“桠溪”);
                               
                                thresh_spec1 = abs (thresh_spec);
                               
                                %标题(与-50分贝的加速度谱图数据阈值”);
                               
                                光谱图(abs (fildat), 256250256200年,“MinThreshold”,-50,“桠溪”)
                               
                                tp_in_min = tp / 60;
                               
                                tp_in_min1 = tp1/60;
                               
                                %线([tp_in_min tp_in_min],[0100],“颜色”,“红”,“线宽”,2);
                               
                                线([tp_in_min1 tp_in_min1], [0100],“颜色”,“绿色”,“线宽”2);
                               
                                网格在
                               
                                网格小
                               
                                轴紧;
                               
                                % print (‘-dpdf’,‘r600’,‘Temp.pdf’)
                               
                                % % Vizualizing地震仪的位移
                               
                                %计算速度与加速度数据使用一个积分器。
                               
                                %积分器与TF过滤器= 1 / (1 - z ^ 1)
                               
                                Nr = 1;
                               
                                博士= [1];
                               
                                %消除趋势加速数据删除漂移、集成和分裂与采样频率速度数据
                               
                                Ver_acc_dt =去趋势(Acc_ver_dlycompensated);
                               
                                NS_acc_dt =去趋势(Acc_NS_dlycompensated);
                               
                                EW_acc_dt =去趋势(Acc_EW_dlycompensated);
                               
                                vel_ver_nodrift =过滤器(Nr,博士,Ver_acc_dt) / Fs;
                               
                                vel_NS_nodrift =过滤器(Nr,博士,NS_acc_dt) / Fs;
                               
                                vel_EW_nodrift =过滤器(Nr,博士,EW_acc_dt) / Fs;
                               
                                %集成速度数据和分裂与采样频率位移数据
                               
                                dis_ver =过滤器(Nr,博士,vel_ver_nodrift) / Fs; A
                               
                                dis_NS =过滤器(Nr,博士,vel_NS_nodrift) / Fs;
                               
                                dis_EW =过滤器(Nr,博士,vel_EW_nodrift) / Fs;
                               
                                %情节位置使用3 d图来可视化地震检波器的动作
                               
                                图(3);
                               
                                plot3 (dis_NS dis_EW dis_ver);
                               
                                网格在;视图([-45,30]);
                               
                                包含(“厘米计算方向”);
                               
                                ylabel (“南方向厘米”);
                               
                                zlabel (在cm中垂直方向的);
                               
                                标题(“3 d地震检波器的位移”);
                               
                                集(gcf,“名字”,“地震检波器轨迹”);
                               
                                集(gcf,“单位”,“归一化”);

PARVATHY NAIR 2023年1月27日

把你的时间,请帮助我。

你可以避免硬件implementaation focud前2 - 3页。

很抱歉打扰你 @Mathieu诺 ,谢谢你的合作

登录置评。

在回答这个问题。

答案(0)

在回答这个问题。

类别

MATLAB 数学图和网络算法

找到更多的在图和网络算法在帮助中心和文件交换

社区寻宝

找到宝藏在MATLAB中央,发现社区如何帮助你!

开始狩猎!

了解一个新提出的STA / LTA算法

7评论
显示6年长的评论隐藏6年长的评论

答案(0)

另请参阅

类别

标签

社区寻宝

了解一个新提出的STA / LTA算法

7评论 显示6年长的评论隐藏6年长的评论

答案(0)

另请参阅

类别

标签

社区寻宝

7评论
显示6年长的评论隐藏6年长的评论