CN111639541A - 基于频率变化率的自适应同步压缩时频分析方法 - Google Patents

基于频率变化率的自适应同步压缩时频分析方法 Download PDF

Info

Publication number
CN111639541A
CN111639541A CN202010368021.4A CN202010368021A CN111639541A CN 111639541 A CN111639541 A CN 111639541A CN 202010368021 A CN202010368021 A CN 202010368021A CN 111639541 A CN111639541 A CN 111639541A
Authority
CN
China
Prior art keywords
frequency
time
adaptive
signal
short
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.)
Pending
Application number
CN202010368021.4A
Other languages
English (en)
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.)
Nanjing University of Science and Technology
Original Assignee
Nanjing University of Science and Technology
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 Nanjing University of Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN202010368021.4A priority Critical patent/CN111639541A/zh
Publication of CN111639541A publication Critical patent/CN111639541A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G06F2218/10—Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
    • 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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Complex Calculations (AREA)
  • Probability & Statistics with Applications (AREA)
  • Operations Research (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Signal Processing (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Discrete Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)

Abstract

本发明公开了一种基于频率变化率的自适应同步压缩时频分析方法，包括：对输入信号进行时域采样，获得信号的离散序列；对信号进行短时傅里叶变换；计算每一个时频点的瞬时频率和瞬时频率变化率的估计值；依据峰值搜索的办法估计信号个数，对信号时频脊线位置的频率变化率进行提取；根据信号的频率变化率，构造新的窗函数；以新的窗函数对信号进行短时傅里叶变换，得到自适应短时傅里叶系数；根据自适应同步压缩公式，得到频率精确估计值；对自适应短时傅里叶系数进行重排，得到时频分布结果。本发明突破了海森伯格不确定原理的限制，使得窗宽可以随着信号的频率变化率自适应调整，使得在整个时频谱中都有很高的分辨率。

Description

基于频率变化率的自适应同步压缩时频分析方法
技术领域
本发明属于信号处理领域，特别是一种基于频率变化率调节窗长的自适应同步压缩时频分析方法。
背景技术
微多普勒信号可以准确反映目标的运动特征和物理结构。多普勒信号大都为非平稳信号，传统的频域分析方法如傅里叶变换无法获得信号时频域信息，因此时频分析是分析该类信号的重要工具。传统的时频分析方法包括短时傅里叶变换(STFT)、魏格纳-威尔分布(WVD)、S变换(ST)、小波变换(WT)等。
Dennis Gabor于1946年提出了短时傅里叶变换，其基本思想是通过加窗实现信号的分段傅里叶变换，从而得到信号的时变特性。但STFT仅仅使用一个确定的窗函数，即STFT不具备自动调节能力。小波变换继承了STFT局部化分析的思想，克服了其窗口固定的缺点，能提供一个随频率变化的窗，但是小波基函数一旦确定后，其特性就固定了，而且小波基设计难度较大，还有容许性条件的约束，同时存在时频分辨率不足、尺度频率转换复杂等缺陷。S变换，引入了可变高斯窗函数，且该函数的时窗宽度与频率成反比。该方法根据频率调节窗长，得到的时频谱在低频部分频率分辨率高，在高频部分频率分辨率低，但是这种反比关系使得窗函数在局部出现窗长过宽和过窄的问题，进而导致在低频处由于时窗宽度过宽引起的时间定位不准确问题，以及在高频处由于频窗宽度过宽引起的频率定位不准确问题。
上述传统的时频分析方法受到海森伯格不确定原理的限制，时频分辨率无法满足信号高分辨率的处理需求，促使一些新兴时频方法的提出。Gaurav Thakur等人于2011年提出同步压缩短时傅里叶变换，其基本思想是把信号的短时傅里叶变换结果在一定频率范围内的时频能量仅沿频率范围“挤压”到信号的瞬时频率附近，从而使信号的时频分辨率显著提升。同步压缩短时傅里叶变换作为时频分析的一种后处理工具，突破了海森伯格不确定原理的限制，但是由于短时傅里叶变换窗口固定，分析频率变换较快的信号存在时频聚集性不高的问题。
专利申请号为CN201510600036.8，发明名称为“一种基于同步压缩分数阶的小波变换方法”的中国专利，首先通过分数阶小波变换得出在最优分数阶上的小波系数,再同步挤压技术对时变信号处理,该方法有较好的时频分析结果，但局限性在于不同信号难以适应固定小波基函数的特性，从而导致时频分辨率变差。
由上可知，现有的时频分析方法还存在不足，需进一步改进来实现高精度的时频分析。
发明内容
本发明的目的在于提供一种基于频率变化率的自适应同步压缩时频分析方法。
实现本发明目的的技术解决方案为：一种基于频率变化率的自适应同步压缩时频分析方法，包括以下步骤：
步骤1、对输入信号进行时域采样，获得信号的离散序列；
步骤2、对信号进行短时傅里叶变换，得到短时傅里叶系数；
步骤3、根据短时傅里叶变换分析结果进一步计算出每一个时频点的瞬时频率和瞬时频率变化率的估计值；
步骤4、依据峰值搜索的办法估计信号个数N,并对信号时频脊线位置的频率变化率进行提取；
步骤5、根据信号的频率变化率，构造自适应窗函数；
步骤6、以自适应窗函数对信号进行短时傅里叶变换，得到自适应短时傅里叶系数；
步骤7、根据自适应同步压缩公式，得到频率精确估计值；
步骤8、根据得到的频率估计对自适应短时傅里叶系数进行重排，得到时频分布结果；
本发明与现有技术相比，其显著优点为：1)本发明可以根据信号频率变化率对窗函数窗宽进行自适应调整，以解决同步压缩在分析频率变换较快的信号存在时频聚集性不高的问题；2)本发明针对不同的输入信号和不同的分辨率要求，能够通过调整控制参数的取值，实现高分辨率时频分析，具有很强的灵活性。
下面结合附图对本发明作进一步详细描述。
附图说明
图1是本发明基于频率变化率的自适应同步压缩时频分析方法的流程图。
图2是本发明实施例1的时频分析结果图。
图3是同步压缩实施例1的时频分析结果图。
图4是本发明实施例2的时频分析结果图。
图5是同步压缩实施例2的时频分析结果图。
具体实施方式
结合图1，一种基于频率变化率的自适应同步压缩时频分析方法，包括以下步骤：
步骤1、对输入信号进行时域采样，获得信号的离散序列；
步骤2、对信号进行短时傅里叶变换，得到短时傅里叶系数；
步骤3、根据短时傅里叶变换分析结果进一步计算出每一个时频点的瞬时频率和瞬时频率变化率的估计值；
步骤4、依据峰值搜索的办法估计信号个数N,并对信号时频脊线位置的频率变化率进行提取；
步骤5、根据信号的频率变化率，构造自适应窗函数；
步骤6、以自适应窗函数对信号进行短时傅里叶变换，得到自适应短时傅里叶系数；
步骤7、根据自适应同步压缩公式，得到频率精确估计值；
步骤8、根据得到的频率估计对自适应短时傅里叶系数进行重排，得到时频分布结果；
进一步地，步骤1对输入信号进行时域采样，获得信号的离散序列；
进一步地，步骤2对信号进行短时傅里叶变换，得到短时傅里叶系数G(t,η)，具体为：
Figure BDA0002477219870000031
其中τ积分变量。
进一步地，步骤3根据短时傅里叶变换分析结果进一步计算出每一个时频点的瞬时频率
Figure BDA0002477219870000032
和瞬时频率变化率
Figure BDA0002477219870000033
的估计值，具体为：
Figure BDA0002477219870000034
Figure BDA0002477219870000035
步骤4、依据峰值搜索的办法估计信号个数N,并对信号时频脊线位置的频率变化率进行提取；
步骤4-1、对于G(t,η)时间轴上任意时间点tk，信号个数N即为峰值的个数，同步压缩系数峰值即满足：
G(tk,η)＞γ，
γ为常数，γ＝max(G(tk,η))/2；
步骤4-2、
(1)对于G(t,η)时间轴上第k个时间点tk，从η＝0开始搜索峰值，并记录位置
(tk,η1)...(tk,ηN)，
(2)
Figure BDA0002477219870000041
(3)重复步骤4-2(1)(2)直至所有时间点完成,得到频率变化率q1(t),q2(t)...qN(t)；
步骤5、根据信号的频率变化率，构造自适应窗函数；
步骤5-1、确定窗宽变化的参数λ；
(1)确定时窗取值范围[σ2,σ1]，其中σ1为最大时窗长度，σ2为最小时窗长度，之后通过不等式确定参数λ的取值范围,不等式为：
Figure BDA0002477219870000042
在(0,1)取值范围中随机选取λ的值；
步骤5-2、将步骤5-1确定的参数代入窗宽变化公式，窗宽变化的公式为：
其中σ1为最大时窗长度，σ2为最小时窗长度，qi(t)为个信号分量的频率变化率，λi为对应参数,N为信号个数；
步骤5-3、将σ(t)代入高斯窗函数构造自适应窗函数gσ(t)为：
Figure BDA0002477219870000043
步骤6、以自适应窗函数对信号进行短时傅里叶变换，得到自适应短时傅里叶系数G1(t,η)；
Figure BDA0002477219870000051
其中τ为积分变量；
步骤7、根据自适应同步压缩公式，得到频率精确估计值
Figure BDA0002477219870000052
Figure BDA0002477219870000053
G1(t,η)自适应短时傅里叶系数，
Figure BDA0002477219870000054
表示窗函数为g2＝(τ/σ2(t))g′(τ/σ(t))的STFT系数；
步骤8、根据得到的频率估计对自适应短时傅里叶系数进行重排，得到时频分布结果；
步骤8-1、对于瞬时频率轴上第l个频率点，计算Δωl,Δωl＝ωl-ωl1，其中l∈[1,n]，n为信号采样总点数，确定重排区间[ωl-Δωl,ωl+Δωl]；
步骤8-2、进行能量重排，得到ASST(t,ω)，计算公式为：
Figure BDA0002477219870000055
其中，ωl为同步压缩后的频率，ηk为自适应短时傅里叶变换谱上的离散化频率点；
步骤8-3、重复步骤8-1、步骤8-2直至所有瞬时频率点完成计算，得到时频结果。
下面以两种信号为例验证本发明的时频分析方法。
实施例1
仿真信号为一频率呈正弦变化的非线性调频信号，解析式为：
s(t)1＝ej2π(6cos(10πt)+240t)t∈(0,1)
信号采样频率fs＝1024Hz，图2为采用自适应同步压缩时频分析方法得到的时频谱。该信号频率随时间呈正弦变化，频率变化较快，取σ1＝64，λ＝0.005，σ2＝30的时频结果中可以清楚地看到频率随时间的变化轨迹，图3为原同步压缩的分析结果。比较图2、图3可以看出自适应同步压缩具有更好的时频分析结果。
实施例2
仿真信号为两个非线性调频信号x1(t)，x2(t)的叠加，解析式为：
Figure BDA0002477219870000061
Figure BDA0002477219870000062
s2(t)＝x1(t)+x2(t)t∈(0,1)
信号采样频率fs＝1024Hz，图4为非线性调频信号经过本实施例方法处理后得到的时频图，取σ1＝90,σ2＝40,λ1＝0.003,λ2＝0.001的时频结果中可以清楚地看到频率随时间的变化轨迹，图5为原同步压缩的分析结果。比较图4、图5可以看出自适应同步压缩具有更好的时频分析结果。
本发明的方法使得窗宽可以随着信号的频率变化率自适应调整，使得在整个时频谱中都有很高的分辨率。

Claims (6)

1.一种基于频率变化率的自适应同步压缩时频分析方法，其特征在于，包括以下步骤：
步骤1、对输入信号进行时域采样，获得信号的离散序列；
步骤2、对信号进行短时傅里叶变换，得到短时傅里叶系数；
步骤3、根据短时傅里叶变换分析结果进一步计算出每一个时频点的瞬时频率和瞬时频率变化率的估计值；
步骤4、依据峰值搜索的办法估计信号个数,并对信号时频脊线位置的频率变化率进行提取；
步骤5、根据信号的频率变化率，构造自适应窗函数；
步骤6、以自适应窗函数对信号进行短时傅里叶变换，得到自适应短时傅里叶系数；
步骤7、根据自适应同步压缩公式，得到频率精确估计值；
步骤8、根据得到的频率估计对自适应短时傅里叶系数进行重排，得到时频分布结果。
2.根据权利要求1所述的基于频率变化率的自适应同步压缩时频分析方法，其特征在于，步骤4依据峰值搜索的办法估计信号个数N，并对信号时频脊线位置的频率变化率进行提取，具体为：
步骤4-1：对于G(t,η)时间轴上任意时间点tk，信号个数N即为峰值的个数，同步压缩系数峰值即满足：
G(tk,η)＞γ
γ为常数，γ＝max(G(tk,η))/2；
步骤4-2：对于G(t,η)时间轴上第k个时间点tk，从η＝0开始搜索峰值，并记录位置(tk,η1)...(tk,ηN)；
步骤4-3：
Figure FDA0002477219860000011
步骤4-4：重复步骤4-2、步骤4-3直至所有时间点完成,得到频率变化率q1(t),q2(t)...qN(t)。
3.根据权利要求2所述的基于频率变化率的自适应同步压缩时频分析方法，其特征在于，步骤5所述根据信号的频率变化率，构造自适应窗函数，具体为：
步骤5-1、确定窗宽变化的参数λ；
(1)确定时窗取值范围[σ2,σ1]，其中σ1为最大时窗长度，σ2为最小时窗长度，之后通过不等式确定参数λ的取值范围,不等式为：
Figure FDA0002477219860000021
在(0,1)取值范围中随机选取λ的值；
步骤5-2、将步骤5-1确定的参数代入窗宽变化公式，窗宽变化的公式为：
Figure FDA0002477219860000022
其中σ1为最大时窗长度，σ2为最小时窗长度，qi(t)为个信号分量的频率变化率，λi为对应参数，N为信号个数；
步骤5-3、将σ(t)代入高斯窗函数构造自适应窗函数gσ(t)为：
Figure FDA0002477219860000023
4.根据权利要求3所述的基于频率变化率的自适应同步压缩时频分析方法，其特征在于，步骤6所述以自适应窗函数对信号进行短时傅里叶变换，得到自适应短时傅里叶系数，具体为：
Figure FDA0002477219860000024
其中τ为积分变量。
5.根据权利要求4所述的基于频率变化率的自适应同步压缩时频分析方法，其特征在于，步骤7所述根据自适应同步压缩公式，得到频率精确估计值，具体为：
Figure FDA0002477219860000025
G1(t,η)自适应短时傅里叶系数，
Figure FDA0002477219860000026
表示窗函数为g2＝(τ/σ2(t))g′(τ/σ(t))的STFT系数。
6.根据权利要求5所述的基于频率变化率的自适应同步压缩时频分析方法，其特征在于，步骤8所述根据得到的频率估计对自适应短时傅里叶系数进行重排，得到时频分布结果，具体为：
步骤8-1、对于瞬时频率轴上第l个频率点，计算Δωl,Δωl＝ωl-ωl-1，其中l∈[1,n]，n为信号采样总点数，确定重排区间[ωl-Δωl,ωl+Δωl]；
步骤8-2、进行能量重排，得到ASST(t,ω)，计算公式为：
Figure FDA0002477219860000031
其中，ωl为同步压缩后的频率，ηk为自适应短时傅里叶变换谱上的离散化频率点；
步骤8-3、重复步骤8-1、步骤8-2直至所有瞬时频率点完成计算，得到时频结果。
CN202010368021.4A 2020-04-30 2020-04-30 基于频率变化率的自适应同步压缩时频分析方法 Pending CN111639541A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010368021.4A CN111639541A (zh) 2020-04-30 2020-04-30 基于频率变化率的自适应同步压缩时频分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010368021.4A CN111639541A (zh) 2020-04-30 2020-04-30 基于频率变化率的自适应同步压缩时频分析方法

Publications (1)

Publication Number Publication Date
CN111639541A true CN111639541A (zh) 2020-09-08

Family

ID=72329961

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010368021.4A Pending CN111639541A (zh) 2020-04-30 2020-04-30 基于频率变化率的自适应同步压缩时频分析方法

Country Status (1)

Country Link
CN (1) CN111639541A (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112834993A (zh) * 2020-12-31 2021-05-25 中国地质大学（武汉） 一种用于ipix雷达信号目标检测的lmsct时频分析方法
CN113197583A (zh) * 2021-05-11 2021-08-03 广元市中心医院 一种基于时频分析和循环神经网络的心电图波形分割方法
CN113280963A (zh) * 2021-05-26 2021-08-20 南通河海大学海洋与近海工程研究院 一种基于改进s变换的实时索力识别方法
CN113822565A (zh) * 2021-09-16 2021-12-21 中国电建集团华东勘测设计研究院有限公司 一种风机监测数据时频特征分级细化分析的方法
CN114417933A (zh) * 2022-01-24 2022-04-29 福州大学 一种基于能量重心法的扫频卷积变换时频分析方法
CN114942419A (zh) * 2022-07-26 2022-08-26 中国石油大学(华东) 长积累时间下的船只散射点三自由度微动特征提取方法
CN116055262A (zh) * 2023-01-18 2023-05-02 中国人民解放军国防科技大学 基于同步挤压小波变换的通信信号载频盲估计方法、系统及介质
CN116055262B (zh) * 2023-01-18 2024-05-28 中国人民解放军国防科技大学 基于同步挤压小波变换的通信信号载频盲估计方法、系统及介质

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112834993A (zh) * 2020-12-31 2021-05-25 中国地质大学（武汉） 一种用于ipix雷达信号目标检测的lmsct时频分析方法
CN113197583A (zh) * 2021-05-11 2021-08-03 广元市中心医院 一种基于时频分析和循环神经网络的心电图波形分割方法
CN113280963A (zh) * 2021-05-26 2021-08-20 南通河海大学海洋与近海工程研究院 一种基于改进s变换的实时索力识别方法
CN113822565A (zh) * 2021-09-16 2021-12-21 中国电建集团华东勘测设计研究院有限公司 一种风机监测数据时频特征分级细化分析的方法
CN113822565B (zh) * 2021-09-16 2023-02-14 中国电建集团华东勘测设计研究院有限公司 一种风机监测数据时频特征分级细化分析的方法
CN114417933A (zh) * 2022-01-24 2022-04-29 福州大学 一种基于能量重心法的扫频卷积变换时频分析方法
CN114942419A (zh) * 2022-07-26 2022-08-26 中国石油大学(华东) 长积累时间下的船只散射点三自由度微动特征提取方法
CN116055262A (zh) * 2023-01-18 2023-05-02 中国人民解放军国防科技大学 基于同步挤压小波变换的通信信号载频盲估计方法、系统及介质
CN116055262B (zh) * 2023-01-18 2024-05-28 中国人民解放军国防科技大学 基于同步挤压小波变换的通信信号载频盲估计方法、系统及介质

Similar Documents

Publication Publication Date Title
CN111639541A (zh) 基于频率变化率的自适应同步压缩时频分析方法
CN109343020B (zh) 一种基于改进窗函数的s变换时频分析方法
US7054792B2 (en) Method, computer program, and system for intrinsic timescale decomposition, filtering, and automated analysis of signals of arbitrary origin or timescale
CN108009347B (zh) 基于同步压缩联合改进广义s变换的时频分析方法
CN110514921B (zh) 一种电力电子变换器非平稳信号中非线性现象的识别方法
CN107085140B (zh) 基于改进的SmartDFT算法的非平衡系统频率估计方法
CN106597408B (zh) 基于时频分析和瞬时频率曲线拟合的高阶pps信号参数估计方法
KR101294681B1 (ko) 기상 신호 처리장치 및 그 처리방법
CN107561497B (zh) Fsk及多种非线性调频信号识别及参数估算方法
CN111912521A (zh) 一种非平稳信号的频率检测方法和存储介质
CN111427018A (zh) 一种雷达干扰装备干扰效果评估方法
US20070063887A1 (en) Method of determining the velocity field of an air mass by high resolution doppler analysis
CN115954017A (zh) 一种基于hht的发动机小样本声音异常故障识别方法及系统
WO2003102810B1 (en) Waveform analysis
CN107966687B (zh) 基于部分自相关谱的mimo雷达信号调制类型识别方法
CN110426681B (zh) 一种基于同步提取s变换的lfm信号参数估计方法
CN107748387A (zh) 一种高分辨的薄互储层含气性检测方法
CN112328956A (zh) 一种强频变信号时频分析方法
CN109120562B (zh) 一种基于频谱累加匹配的mfsk信号频率估计方法
CN113552543B (zh) 基于set-stiaa的空间微动目标时频分析方法
CN115630271A (zh) 信号频率估计方法、装置、设备及存储介质
CN112883787B (zh) 一种基于频谱匹配的短样本低频正弦信号参数估计方法
CN114674410A (zh) 一种分量数时变的水声信号瞬时频率估计方法
CN113553901A (zh) 一种基于自适应窗长的改进同步提取时频分析方法
Khurram et al. Detection of oscillatory modes in power systems using empirical wavelet transform

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20200908

RJ01 Rejection of invention patent application after publication