CN107992448B - 一种基于绝对值的直接反余弦瞬时频率求解方法 - Google Patents

一种基于绝对值的直接反余弦瞬时频率求解方法 Download PDF

Info

Publication number
CN107992448B
CN107992448B CN201711238456.1A CN201711238456A CN107992448B CN 107992448 B CN107992448 B CN 107992448B CN 201711238456 A CN201711238456 A CN 201711238456A CN 107992448 B CN107992448 B CN 107992448B
Authority
CN
China
Prior art keywords
instantaneous
frequency
instantaneous frequency
inverse cosine
method based
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.)
Active
Application number
CN201711238456.1A
Other languages
English (en)
Other versions
CN107992448A (zh
Inventor
张庆
贾林山
何小高
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201711238456.1A priority Critical patent/CN107992448B/zh
Publication of CN107992448A publication Critical patent/CN107992448A/zh
Application granted granted Critical
Publication of CN107992448B publication Critical patent/CN107992448B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Measuring Frequencies, Analyzing Spectra (AREA)

Abstract

一种基于绝对值的直接反余弦瞬时频率求解方法，首先将调频信号进行取反余弦操作，得到非单调的瞬时相位，然后直接对非单调的瞬时相位进行前向差分，得到包含有负瞬时频率的初始瞬时频率，对包含有负瞬时频率的初始瞬时频率取绝对值，得到全部为正值的瞬时频率，使用三阶中值滤波器对全部为正值的瞬时频率进行滤波，最终得到调频信号的瞬时频率；本发明相较于传统的基于相位展开的直接反余弦法，具有过程简单、运算效率高的优点。

Description

一种基于绝对值的直接反余弦瞬时频率求解方法
技术领域
本发明属于瞬时频率求解方法技术领域，具体涉及一种基于绝对值的直接反余弦瞬时频率求解方法。
背景技术
瞬时频率是调幅-调频(AM-FM)信号的重要调制参数，对于探究非平稳、非线性过程的详细机制具有重要意义。目前，瞬时频率的求解思路是先将信号分解成为一系列单分量信号，然后对这些单分量信号进行解调。比较常用的包括Hilbert-Huang变换(NordenE.Huang.etc.The empirical mode decomposition and the Hilbert spectrum fornonlinear and non-stationary time series analysis[J].Proceedings of the royalsociety A:Mathematical,Physical and Engineering Sciences.1998:pp903-995.)，经验AM-FM分解(Norden E.Huang et al,on instantaneous frequency,Advances inAdaptive Data Analysis 1(2),2009:pp177-229.)，局部均值分解(Jonathan S.Smith,The local mean decomposition and its application to EEG perception data,J.R.Soc.Interface 2(5),2005:pp443–454.)三种方式。Hilbert-Huang变换是首先使用经验模式分解将信号分解为一系列本征模态函数，然后使用希尔伯特变换求取本征模态函数的瞬时频率；经验AM-FM分解和Hilbert-Huang变换一样首先将信号分解为本征模式函数，然后将本征模式函数分解为调频部分的包络函数和调幅部分的纯调频信号，然后使用直接正交法对调频部分的纯调频信号求解瞬时频率；局部均值分解是将信号分解成为一系列的乘积函数，每一乘积函数是由调幅部分的包络函数和调频部分的纯调频信号组成，使用直接反余弦法求取乘积函数的纯调频信号的瞬时频率。
传统的直接反余弦法是首先对纯调频信号求反余弦得到非单调相位，然后对相位进行展开得到展开后的单调递增的瞬时频率，最后对展开后的相位求导得到瞬时频率。直接反余弦法是一种性能比较好的纯调频信号求解方法，相较于Hilbert变换、Teager能量法等其他求解方法，具有无端点效应、求解速度快、不会出现负频率的优点。但是，直接反余弦法也存在一些问题。在使用反余弦函数，由于反余弦函数是定义在区间[-1,1]上的单调函数，而纯调频信号则是非单调函数，那么对纯调频信号求反余弦后，得到的相位函数就是非单调的。根据瞬时频率的定义，瞬时频率应当是相位的导数，且瞬时频率应当都是正的，于是，将非单调的相位函数展开为单调递增的瞬时相位是获取瞬时频率的关键一步。针对这一问题，有文献(任达千,杨世锡,吴昭同,等.信号瞬时频率直接计算法与Hilbert变换及Teager能量法比较[J].机械工程学报,2013,49(9):42-48.)设定了一系列十分复杂的判定规则，将相位展开；也有文献(张亢,程军圣,杨宇,等.基于分段波形的信号瞬时频率计算方法[J].湖南大学学报(自科版),2011,38(11):54-59.)将纯调频信号分为“全波”段，然后对每一个全波段中的相位函数进行调整，使得每一个全波段的瞬时相位是单调递增的，最终获得完整的单调递增的瞬时相位。这些处理方法都十分繁杂，不易于理解且由于判别过程繁杂，影响了运算效率。
发明内容
为了克服上述现有技术的缺点，本发明的目的在于提供一种基于绝对值的直接反余弦瞬时频率求解方法，在保证求解结果正确的情况下，能有效降低算法复杂度，提高运算效率。
为了达到上述目的，本发明采取的技术方案是：
一种基于绝对值的直接反余弦瞬时频率求解方法，包括以下步骤：
1)设有调频信号s(t)，信号长度为N，满足-1≤s(t)≤1，先将调频信号s(t)取反余弦，得到非单调的瞬时相位
Figure BDA0001489449640000021
Figure BDA0001489449640000022
2)对非单调的瞬时相位
Figure BDA0001489449640000023
进行前向差分，设信号采样频率为fs，得到包含有负瞬时频率的初始瞬时频率f0(t)：
Figure BDA0001489449640000031
上式中diff表示前向差分运算；
3)对初始瞬时频率f0(t)取绝对值，得到全部为正值的瞬时频率f1(t)；
f1(t)＝abs(f0(t))
上式中，abs表示取绝对值操作；
4)使用三阶中值滤波器对全部为正值的瞬时频率f1(t)进行滤波，最终得到调频信号s(t)的瞬时频率f(t)，
f(t)＝med3(f1(t))
上式中，med3表示三阶中值滤波操作。
本发明的有益效果为：
(1)本发明采用直接反余弦法求解纯调频信号瞬时频率，探明了相位展开并不是获取正的瞬时频率的必要条件，绕开了繁琐的相位展开过程，降低了直接反余弦法的使用门槛。
(2)本发明与现有的基于相位展开的直接反余弦法相比，具有算法简单，运算效率高的优点，更适合在线运算。
附图说明
图1是本发明方法流程图。
图2是实施例调频信号s(t)的波形图。
图3是实施例调频信号s(t)的非单调的瞬时相位
Figure BDA0001489449640000032
图4是实施例调频信号s(t)的包含有负瞬时频率的初始瞬时频率f0(t)。
图5是实施例调频信号s(t)的全部为正值的瞬时频率f1(t)。
图6是使用三阶中值滤波器对f1(t)处理得到的无毛刺瞬时频率f(t)与真实瞬时频率对照图。
图7是使用传统的基于相位展开的直接反余弦法求得的s(t)的无毛刺瞬时频率fold(t)与真实瞬时频率对照图。
具体实施方式
以下结合附图和实施例作对本发明进一步的详细说明。
参照图1，一种基于绝对值的直接反余弦瞬时频率求解方法，包括以下步骤：
1)产生一个仿真信号作为调频信号s(t)，公式为：
Figure BDA0001489449640000041
设定采样频率为2Hz，t∈[1,300]，其波形如图2所示；将调频信号s(t)取反余弦，得到非单调的瞬时相位
Figure BDA0001489449640000042
如图3所示；
2)对非单调的瞬时相位
Figure BDA0001489449640000043
进行前向差分，得到包含有负瞬时频率的初始瞬时频率f0(t)，如图4所示；
3)对初始瞬时频率f0(t)取绝对值，得到全部为正值的瞬时频率f1(t)，如图5所示；
4)由于全部为正值的瞬时频率f1(t)尚存在许多毛刺，这些由直接反余弦法在极值点附近存在计算误差导致的，需要滤除，这里使用三阶中值滤波器对毛刺进行滤除，最终得到调频信号s(t)的调频信号s(t)的无毛刺瞬时频率f(t)，将f(t)与真实瞬时频率绘制在一幅图中，如图6所示；
排除绘图代码，上述运算过程共耗时0.000165秒。
对比例：为作对比，使用传统的基于相位展开的直接反余弦法(任达千,杨世锡,吴昭同,等.信号瞬时频率直接计算法与Hilbert变换及Teager能量法比较[J].机械工程学报,2013,49(9):42-48.)处理与本实施例相同的调频信号s(t)，该方法首先对纯调频信号应用反余弦函数求相位，这一步和本发明是一致的，即本发明中的非单调的瞬时相位
Figure BDA0001489449640000044
然后对非单调的瞬时相位
Figure BDA0001489449640000045
进行展开，求解得到单调递增的瞬时相位，然后对单调递增的瞬时相位进行前向差分最终可以得到瞬时频率调频信号s(t)的瞬时频率fold(t)，将fold(t)与真实瞬时频率绘制在一幅图中，如图7所示；
排除绘图代码，本过程共耗时0.001895秒。
对比图6和图7中可以看出，两种方法的计算结果都与真实结果非常接近。为了量化所求得的瞬时频率与真实值之间的误差，使用均方误差(MSE)来描述误差，其计算公式为：
Figure BDA0001489449640000051
上式中，f(t)表示求得的瞬时频率，在本实施例中分别为使用本发明提出的方法求得的调频信号s(t)的瞬时频率f(t)以及使用传统的基于相位展开的直接反余弦法求得的瞬时频率fold(t)，fr(t)为真实的瞬时频率，n为序列长度。MSE越接近于0说明求解结果越精确。
针对本实施例中的仿真信号25秒到275秒信号区段，求出本发明的瞬时频率的MSE为1.7518，传统的基于相位展开的直接反余弦法的瞬时频率的MSE为1.7513。对比可以发现两者的MSE在小数点后3位相等，说明精度相当。
通过以上对比，可以看出，本发明与传统的基于相位展开的直接反余弦瞬时频率求解方法相比，在保证精度的前提下，前者的运行效率是后者的10倍以上，这是由于本发明无需对逐一相位展开，减少了大量的运算，也即体现了本发明在计算效率上有明显优势。
需要说明的是，这两个过程和结果都是在同一台计算机环境下使用MATLAB2016B软件得到的，具体计算机环境如下所述：
CPU：I7-6700HQ 2.6GHz；
RAM：8GB；
ROM：128GB SSD+1TB HDD；
GPU：GTX960。

Claims (1)

1.一种基于绝对值的直接反余弦瞬时频率求解方法，其特征在于，包括以下步骤：
1)设有调频信号s(t)，信号长度为N，满足-1≤s(t)≤1，先将调频信号s(t)取反余弦，得到非单调的瞬时相位
Figure FDA0002263768900000011
Figure FDA0002263768900000012
2)对非单调的瞬时相位
Figure FDA0002263768900000013
进行前向差分，设信号采样频率为fs，得到包含有负瞬时频率的初始瞬时频率f0(t)：
Figure FDA0002263768900000014
上式中diff表示前向差分运算；
3)对初始瞬时频率f0(t)取绝对值，得到全部为正值的瞬时频率f1(t)：
f1(t)＝abs(f0(t))
上式中，abs表示取绝对值操作；
4)使用三阶中值滤波器对全部为正值的瞬时频率f1(t)进行滤波，最终得到调频信号s(t)的瞬时频率f(t)：
f(t)＝med3(f1(t))
上式中，med3表示三阶中值滤波操作。
CN201711238456.1A 2017-11-30 2017-11-30 一种基于绝对值的直接反余弦瞬时频率求解方法 Active CN107992448B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711238456.1A CN107992448B (zh) 2017-11-30 2017-11-30 一种基于绝对值的直接反余弦瞬时频率求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711238456.1A CN107992448B (zh) 2017-11-30 2017-11-30 一种基于绝对值的直接反余弦瞬时频率求解方法

Publications (2)

Publication Number Publication Date
CN107992448A CN107992448A (zh) 2018-05-04
CN107992448B true CN107992448B (zh) 2020-03-31

Family

ID=62034595

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711238456.1A Active CN107992448B (zh) 2017-11-30 2017-11-30 一种基于绝对值的直接反余弦瞬时频率求解方法

Country Status (1)

Country Link
CN (1) CN107992448B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108628801B (zh) * 2018-05-10 2020-07-24 西安交通大学 一种改进型直接正交瞬时频率求解方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0618672A1 (de) * 1993-03-31 1994-10-05 ANT Nachrichtentechnik GmbH Verfahren zur Demodulation von frequenzmodulierten Signalen
CN101887407A (zh) * 2010-07-16 2010-11-17 哈尔滨工业大学 一种基于希尔伯特-黄变换的设备或系统机内测试信号特征提取方法
CN103973411A (zh) * 2014-04-10 2014-08-06 沈阳理工大学 一种差分空时编码信号盲检测方法
CN105044456A (zh) * 2015-07-21 2015-11-11 电子科技大学 一种基于正交子带的电网瞬时频率测量与跟踪方法
CN107091988A (zh) * 2017-06-15 2017-08-25 西安交通大学 一种三相电机电流信号瞬时频率和瞬时包络提取方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0618672A1 (de) * 1993-03-31 1994-10-05 ANT Nachrichtentechnik GmbH Verfahren zur Demodulation von frequenzmodulierten Signalen
CN101887407A (zh) * 2010-07-16 2010-11-17 哈尔滨工业大学 一种基于希尔伯特-黄变换的设备或系统机内测试信号特征提取方法
CN103973411A (zh) * 2014-04-10 2014-08-06 沈阳理工大学 一种差分空时编码信号盲检测方法
CN105044456A (zh) * 2015-07-21 2015-11-11 电子科技大学 一种基于正交子带的电网瞬时频率测量与跟踪方法
CN107091988A (zh) * 2017-06-15 2017-08-25 西安交通大学 一种三相电机电流信号瞬时频率和瞬时包络提取方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
The local mean decomposition and its application to EEG perception data;Jonathan S.Smith;《Journal of the Royal Society interface》;20050728;第443-454页 *
一种新的估计瞬时频率的方法-经验包络法;郑近德 等;《振动与冲击》;20121231;第31卷(第17期);第86-90页 *
信号瞬时频率直接计算法与Hilbert变换及Teager能量法比较;任达千;《机械工程学报》;20130531;第49卷(第9期);第42-48页 *

Also Published As

Publication number Publication date
CN107992448A (zh) 2018-05-04

Similar Documents

Publication Publication Date Title
CN101919695B (zh) 一种基于小波变换的心电信号qrs波检测方法
CN106685478B (zh) 基于信号时频图像信息提取的跳频信号参数估计方法
CN105302940B (zh) 一种基于循环相关熵的载频估计方法
Qin et al. Research on iterated Hilbert transform and its application in mechanical fault diagnosis
Tang et al. Hilbert-Huang transform for ECG de-noising
CN102508206A (zh) 基于小波包去噪和功率谱熵的线性调频信号参数估计方法
CN110909480B (zh) 一种水轮机振动信号的去噪方法与装置
CN107992448B (zh) 一种基于绝对值的直接反余弦瞬时频率求解方法
CN109347452A (zh) 基于分段线性函数的双频功放数字预失真装置及方法
CN105721371B (zh) 一种基于循环谱相关的常用数字调制信号识别方法
CN103308829B (zh) 一种gis单次局放信号提取与触发时刻调整方法
CN102832942A (zh) 基于分数阶Fourier变换的三角线性调频连续波特征提取方法
CN108090270B (zh) 一种基于形态学滤波和盲源分离的暂态振荡参数识别方法
CN109586728B (zh) 基于稀疏贝叶斯的调制宽带转换器框架下信号盲重构方法
CN107342750B (zh) 适用于多奈奎斯特区的分数延迟优化方法及其实现结构
CN108363994A (zh) 基于经验模态分解改进的乘性噪声去除技术
CN105022917B (zh) 一种信号精确提取与处理方法
CN109617051B (zh) 一种新能源电力系统低频振荡参数辨识方法
CN104065597A (zh) 一种基于小波能量分布熵的bpsk/qpsk信号识别方法
CN113627398B (zh) 一种基于自适应重构滤波的信号特征检测方法
Boudraa Instantaneous frequency estimation of fm signals by ψb-energy operator
CN104331862A (zh) 一种并联型分数阶零相位滤波器及其滤波方法
CN108134604B (zh) 一种基于直接法求解瞬时频率的毛刺定位与消除方法
CN108628801B (zh) 一种改进型直接正交瞬时频率求解方法
CN117454155B (zh) 基于ssaf和emd的igbt声发射信号提取方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant