CN104950335B - Enpemf信号归一化stft‑wvd时频分析方法 - Google Patents
Enpemf信号归一化stft‑wvd时频分析方法 Download PDFInfo
- Publication number
- CN104950335B CN104950335B CN201510209896.9A CN201510209896A CN104950335B CN 104950335 B CN104950335 B CN 104950335B CN 201510209896 A CN201510209896 A CN 201510209896A CN 104950335 B CN104950335 B CN 104950335B
- Authority
- CN
- China
- Prior art keywords
- stft
- wvd
- arrays
- array
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提供了一种ENPEMF信号归一化STFT‑WVD时频分析方法，对地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner‑Ville分布，得到STFT数组和WVD数组，然后选取STFT数组中的最大值归一化得到数组STFT_1，记录1值位置及最小值，将其中0值的数用最小值替换，选取WVD数组中相同位置数归一化得到临时数组A，点除以STFT‑1得到临时数组B，将其中大于设置值x的位置的数与WVD数组中对应位置的数全部置为0，输出临时数组B和WVD数组。本发明较好地消除交叉项干扰，沿袭WVD较高的时频分辨率，克服输入信号改变从而需要重新调节阀值与幂调节系数的缺点，结果理想，使用灵活。
Description
技术领域
本发明涉及一种ENPEMF信号归一化STFT-WVD时频分析方法，属于地震电磁前兆研究与地震预测领域。
背景技术
地震给人类的生活带来了巨大的灾难，据统计，全球的自然灾害之中，地震造成的死亡人数占全部自然灾害死亡人数的54％，堪称自然灾害之最。如何预测地震一直以来都是一个热门敏感的话题。然而，因地震预测有着地球内部的“不可入性”、大地震的“非频发性”、“地震物理过程的复杂性”三大困难，地震预测成为了公认的世界性难题，对于地震前兆预测对于人类生命安全和社会财产安全的具有很大的意义。
在现有的孕震信息研究过程中，STFT与WVD时频分析方法是常用的用来对采集到的大量地球天然脉冲电磁场信号进行分析的方法。
传统的STFT-WVD时频分析方法会用到以下公式对地球天然脉冲电磁场信号进行处理：
SWx(t,f)＝Wx(t,f){|SSTFTx(t,f)|>c}………………(1)
SWx(t,f)＝SSTFTx a(t,f)Wx b(t,f)…………………(2)
传统的STFT-WVD时频分析方法需要设置阀值c以及幂调节系数a和b。其中c为STFT谱的交叉项消除阀值。当STFT谱的值小于该阀值时返回0，如果大于该阀值则返回1。WVD中有交叉项对应STFT谱部分的数值肯定小于该阀值，用数字0与WVD相乘以消除交叉项；其中a和b式为幂调节系数，作用是增强两变换数值较大部分而消弱有交叉项部分。式(1)所示方法灵活性差，输入信号幅值或能量大小直接影响c值的选择，目前没有人能提出一种行之有效的自适应阀值选择方法，并且在实际信号中有用分量幅值、能量大小往往各不相同甚至差别较大，因此采用设置阀值消交叉项容易也将信息项消除；式(2)所示方法有些许改进，但同样幂调节系数的确定没有理论基础，如何根据待分析信号的特征确定幂调节系数有待进一步研究，对使用造成不便，同时存在难于避免的交叉项干扰，往往需要进一步利用滤波等方法对得到的信号进行进一步处理，才能得到便于解读的时频图和谱图。
发明内容
为了解决现有技术的不足，本发明提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，通过结合短时Fourier变换和Wigner-Ville分布各自特点较好地消除了交叉项的干扰，同时沿袭了WVD较高的时频分辨率。
本发明为解决其技术问题所采用的技术方案是：提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，包括以下步骤：
(1)将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布化，分别得到STFT数组和WVD数组；
(2)选取STFT数组的最大值max_st，将STFT数组中的各个数除以max_st以对STFT数组进行归一化，得到归一化后的数组STFT_1；
(3)记录数组STFT_1数值为1的数所在的位置(i，j)；记录数组STFT_1中非0值的最小值min_1；
(4)将数组STFT_1中的值为0的数全部用min_1替换；
(5)选取WVD数组中位置为(i，j)的数max_wvd，将WVD数组中的各个数除以max_wvd以对WVD数组进行归一化，得到临时数组A；
(6)临时数组A点除以数组STFT_1得到临时数组B，设置矩阵倍数比值上限值x，x的范围为1和2之间；选取临时数组B中大于x的数并记录它们的位置，将临时数组B中大于x的数全部置为0，将WVD数组中与临时数组B中大于x的数的相同位置的数全部置为0；
(7)输出经步骤(6)置零后的临时数组B以及WVD数组。
步骤(6)中矩阵倍数比值上限值x根据地球天然脉冲电磁场信号的幅值强弱进行设置，x随地球天然脉冲电磁场信号的幅值变强而增大。
步骤(7)输出的WVD数组进行二维时频谱显示。
本发明基于其技术方案所具有的有益效果在于：
本发明的ENPEMF信号归一化STFT-WVD时频分析方法结合短时Fourier变换和Wigner-Ville分布各自特点较好地消除了交叉项的干扰，同时沿袭了WVD较高的时频分辨率，与传统的STFT-WVD算法相比，本发明的ENPEMF信号归一化STFT-WVD时频分析方法不需要设置阀值或者调节幂调节系数，克服了输入信号改变从而需要重新调节阀值与幂调节系数的缺点，结果更加理想，使用更加灵活。
附图说明
图1是本发明的流程示意图。
图2是线性调频信号Wigner-Ville分布二维时频能量分布示意图。
图3是线性调频信号短时Fourier变换二维时频谱示意图。
图4是线性调频信号归一化STFT_WVD二维时频能量分布示意图。
图5是第18日通道3NH数据短时Fourier变换二维时频谱示意图。
图6是第18日通道3NH数据Wigner-Ville分布二维时频能量分布示意图。
图7是第18日通道3NH数据归一化STFT_WVD二维时频能量分布。
具体实施方式
下面结合附图和实施例对本发明作进一步说明。
本发明提供了一种ENPEMF信号归一化STFT-WVD时频分析方法，结合图1，包括以下步骤：
(1)利用matlab，将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布化，所述两种变换的函数分别为：
[ST,Ts,F]＝stft(x,Nw,nstep,h,dt)
[tfr,t,f]＝wignerVille(x,fs,fre_bins)
分别得到STFT数组和WVD数组，其中ST和tfr分别表示STFT数组和WVD数组；
(2)选取STFT数组的最大值max_st：
max_st＝max(max(ST))
将STFT数组中的各个数除以max_st以对STFT数组进行归一化，得到归一化后的数组STFT_1：
STFT_1＝ST/max_st
(3)记录数组STFT_1数值为1的数所在的位置(i，j)：
[i,j,v]＝find(STFT_1＝＝1)
记录数组STFT_1中非0值的最小值min_1：
[a,b,c]＝find(STFT_1)
min_1＝min(c)
(4)将数组STFT_1中的值为0的数全部用min_1替换：
STFT_1(STFT_1＝＝0)＝min_1
(5)选取WVD数组中位置为(i，j)的数max_wvd：
max_wvd＝tfr(i,j)
将WVD数组中的各个数除以max_wvd以对WVD数组进行归一化，得到临时数组A：
A＝tfr/max_wvd
(6)临时数组A点除以数组STFT_1得到临时数组B：
B＝A./STFT_1
根据地球天然脉冲电磁场信号的幅值强弱设置阵倍数比值上限值x，x随地球天然脉冲电磁场信号的幅值变强而增大，x可以设置为1.8，选取临时数组B中大于x的数并记录它们的位置，将临时数组B中大于x的数全部置为0：
tfr(B>1.8)＝0
将WVD数组中与临时数组B中大于x的数的相同位置的数全部置为0：
B(B>1.8)＝0
(7)输出经步骤(6)置零后的临时数组B以及WVD数组。
步骤(7)输出的WVD数组进行二维时频谱显示。
对于任何一种时频分布方法，有一个公认的观点：如果该时频分布方法对线性调频信号不能提供好的时频聚集性，那么它便不适合用作非平稳信号时频分析的工具。图2所示为仅经过Wigner-Ville分布化的线性调频信号的二维时频能量分布，图3所示为仅经过短时Fourier变换后的线性调频信号二维时频谱示意图，图4是利用本发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图。与图3相比利用本发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图克服了短时Fourier变换较差的时频聚集特性的缺点，和图2相比本发明的ENPEMF信号归一化STFT-WVD时频分析方法解决了Wigner-Ville分布固有交叉项的干扰的问题，而且在对信号处理前后不需要对参数进行设置、调整，因此本发明的ENPEMF信号归一化STFT-WVD时频分析方法是一种行之有效的时频分析工具。
18日芦山Ms7.0地震前通道3的NH数据(每秒接收的地球天然脉冲电磁场信号脉冲个数)短时Fourier变换的二维时频谱如图5所示，经过Wigner-Ville分布的二维时频能量分布示意图如图6所示，采用本发明的ENPEMF信号归一化STFT-WVD时频分析方法得到的二维时频能量分布示意图如图7所示。可以看到与使用短时Fourier变换得到的二维时频谱相比采用归一化STFT_WVD方法得到的二维时频能量分布图更加清晰的描述了信号在不同时间和频率上的能量、强度和信号能量集中的时间和频段；与使用Wigner-Ville分布得到的二维时频能量分布图相比采用归一化STFT_WVD方法得到的二维时频能量分布图排除了交叉项信号叠加对真实信息项在不同时间和频率上的能量的反映的干扰，对信号有用信息分量的提取和判断起重要作用。
Claims (3)
1.一种ENPEMF信号归一化STFT-WVD时频分析方法，其特征在于包括以下步骤：
(1)将地球天然脉冲电磁场信号分别作短时Fourier变换和Wigner-Ville分布化，分别得到STFT数组和WVD数组；
(2)选取STFT数组的最大值max_st，将STFT数组中的各个数除以max_st以对STFT数组进行归一化，得到归一化后的数组STFT_1；
(3)记录数组STFT_1数值为1的数所在的位置(i，j)；记录数组STFT_1中非0值的最小值min_1；
(4)将数组STFT_1中的值为0的数全部用min_1替换；
(5)选取WVD数组中位置为(i，j)的数max_wvd，将WVD数组中的各个数除以max_wvd以对WVD数组进行归一化，得到临时数组A；
(6)临时数组A点除以数组STFT_1得到临时数组B，设置矩阵倍数比值上限值x，x的范围为1和2之间；选取临时数组B中大于x的数并记录它们的位置，将临时数组B中大于x的数全部置为0，将WVD数组中与临时数组B中大于x的数的相同位置的数全部置为0；
(7)输出经步骤(6)置零后的临时数组B以及WVD数组。
2.根据权利要求1所述的ENPEMF信号归一化STFT-WVD时频分析方法，其特征在于：步骤(6)中矩阵倍数比值上限值x根据地球天然脉冲电磁场信号的幅值强弱进行设置，x随地球天然脉冲电磁场信号的幅值变强而增大。
3.根据权利要求1所述的ENPEMF信号归一化STFT-WVD时频分析方法，其特征在于：步骤(7)输出的WVD数组进行二维时频谱显示。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510209896.9A CN104950335B (zh) | 2015-04-28 | 2015-04-28 | Enpemf信号归一化stft‑wvd时频分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510209896.9A CN104950335B (zh) | 2015-04-28 | 2015-04-28 | Enpemf信号归一化stft‑wvd时频分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104950335A CN104950335A (zh) | 2015-09-30 |
| CN104950335B true CN104950335B (zh) | 2017-05-31 |
Family
ID=54165130
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510209896.9A Expired - Fee Related CN104950335B (zh) | 2015-04-28 | 2015-04-28 | Enpemf信号归一化stft‑wvd时频分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104950335B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107144874A (zh) * | 2017-07-12 | 2017-09-08 | 中国地质大学（武汉） | 一种基于对enpemf信号进行bswt‑ddtfa时频分析的方法及系统 |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106483555B (zh) * | 2016-09-28 | 2018-05-22 | 中国地质大学（武汉） | Enpemf数据的格林函数-spwvd时频分析方法 |
| CN107037486B (zh) * | 2017-03-31 | 2019-01-08 | 中国地质大学（武汉） | 地球天然脉冲电磁场数据处理的时频谱分析方法及系统 |
| CN107895141B (zh) * | 2017-10-23 | 2020-02-14 | 中国地质大学（武汉） | 一种用于含噪enpemf信号的镜像累加nnmp-sst时频分析方法 |
| CN107831549A (zh) * | 2017-11-20 | 2018-03-23 | 中国地质大学（武汉） | 一种enpemf信号的nmp倒谱sst时频方法 |
| CN109374968A (zh) * | 2018-12-14 | 2019-02-22 | 国网山东省电力公司电力科学研究院 | 一种基于stft-wvd变换的vfto频谱分析方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6512788B1 (en) * | 1998-11-02 | 2003-01-28 | Agilent Technologies, Inc. | RF output spectrum measurement analyzer and method |
| CN101944235B (zh) * | 2009-09-18 | 2012-02-01 | 哈尔滨工程大学 | 基于分数傅立叶变换的图像压缩方法 |
| CN102466819B (zh) * | 2010-11-03 | 2014-04-16 | 中国石油天然气集团公司 | 地震信号的频谱分析方法及装置 |
-
2015
- 2015-04-28 CN CN201510209896.9A patent/CN104950335B/zh not_active Expired - Fee Related
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107144874A (zh) * | 2017-07-12 | 2017-09-08 | 中国地质大学（武汉） | 一种基于对enpemf信号进行bswt‑ddtfa时频分析的方法及系统 |
| CN107144874B (zh) * | 2017-07-12 | 2019-02-05 | 中国地质大学（武汉） | 一种基于对enpemf信号进行bswt-ddtfa时频分析的方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104950335A (zh) | 2015-09-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104950335B (zh) | Enpemf信号归一化stft‑wvd时频分析方法 | |
| CN112904414B (zh) | 地声事件定位及失稳灾害预警方法、感知仪、监测系统 | |
| CN107991706B (zh) | 基于小波包多重阈值和改进经验模态分解的煤层水力压裂微震信号联合降噪方法 | |
| CN108345033B (zh) | 一种微地震信号时频域初至检测方法 | |
| CN101872616B (zh) | 端点检测方法以及使用该方法的系统 | |
| CN107316653B (zh) | 一种基于改进的经验小波变换的基频检测方法 | |
| CN107144879B (zh) | 一种基于自适应滤波与小波变换结合的地震波降噪方法 | |
| CN107300971B (zh) | 基于骨传导振动信号传播的智能输入方法及系统 | |
| CN106886044A (zh) | 一种基于剪切波与Akaike信息准则的微地震初至拾取方法 | |
| CN107688553A (zh) | 基于小波变换和逻辑回归算法检测心电波形特征的方法 | |
| CN104887263A (zh) | 一种基于心音多维特征提取的身份识别算法及其系统 | |
| CN105445801B (zh) | 一种消除二维地震资料随机噪音的处理方法 | |
| CN102799892A (zh) | 一种mfcc水下目标特征提取和识别方法 | |
| CN103116405A (zh) | 牙齿动作状态的脑肌电的实时检测控制装置及方法 | |
| CN107037486A (zh) | 地球天然脉冲电磁场数据处理的时频谱分析方法及系统 | |
| CN108520759B (zh) | 用于帕金森病语音检测的时频特征图像提取方法 | |
| CN105259579A (zh) | 一种基于地震数据瞬时属性的强振幅屏蔽层剔除方法 | |
| CN109377982B (zh) | 一种有效语音获取方法 | |
| CN109186752A (zh) | 基于图形处理器的水下声学信号采集、传输和检测系统 | |
| CN103578466A (zh) | 基于分数阶傅里叶变换的语音非语音检测方法 | |
| CN103308829A (zh) | 一种gis单次局放信号提取与触发时刻调整方法 | |
| CN110146922A (zh) | 高速铁路地震预警系统单双地震计干扰识别方法 | |
| CN106094033A (zh) | 奇异值分解的定向地震波束形成方法 | |
| CN113095113B (zh) | 一种用于水下目标识别的小波线谱特征提取方法及系统 | |
| CN103778914B (zh) | 基于信噪比加权模板特征匹配的抗噪语音识别方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170531 Termination date: 20190428 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |