CN112505413B - 一种时频分析方法和系统 - Google Patents

一种时频分析方法和系统 Download PDF

Info

Publication number
CN112505413B
CN112505413B CN202011337810.8A CN202011337810A CN112505413B CN 112505413 B CN112505413 B CN 112505413B CN 202011337810 A CN202011337810 A CN 202011337810A CN 112505413 B CN112505413 B CN 112505413B
Authority
CN
China
Prior art keywords
signal
signals
frequency domain
frequency
decomposition
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
CN202011337810.8A
Other languages
English (en)
Other versions
CN112505413A (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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong 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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN202011337810.8A priority Critical patent/CN112505413B/zh
Publication of CN112505413A publication Critical patent/CN112505413A/zh
Application granted granted Critical
Publication of CN112505413B publication Critical patent/CN112505413B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
    • G01R23/16—Spectrum analysis; Fourier analysis

Landscapes

  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Resistance Or Impedance (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种时频分析方法和系统。该方法包括步骤：对输入信号进行第一次信号分解；对第一次信号分解得到的信号进行第二次信号分解，获得三个长度相等的信号；对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；根据虚部求和结果计算信号瞬时频率；根据计算的瞬时频率及实部求和结果计算信号的直流分量、原点相位和幅值。本发明实现了瞬时频率、直流分量、瞬时幅值和瞬时相位的计算，适用于平稳信号和非平稳信号，计算精度高，对信号长度几乎无要求，抗直流干扰能力强，计算时间短。

Description

一种时频分析方法和系统
技术领域
本发明属于信号分析技术领域，更具体地，涉及一种时频分析方法和系统。
背景技术
信号是信息的载体，人们认识世界万事万物及通讯均基于对信息的理解。获取信号中的信息依赖于信号分析。通常实际信号均为时变信号，科学家们多采用非平稳信号代指时变信号。时频分析被广泛应用于分析非平稳信号，如分析引力波，动物的语言，脑电波等信号。随着人类认识世界的深入，对时频分析方法的要求越来越高。
现有技术中存在多种时频分析方法。例如线性时频描述法(Linear time-frequency representation)、双线性时频分布法(Bilinear time-frequencydistribution)、时变高阶光谱法(Time-varying higher order spectra)、自适应参数时频分析法(Adaptive parametric time-frequency analysis)、时频ARMA模型法 (Time-frequency ARMA models)和自适应非参数时频分析法(Adaptive non-parametric time-frequency analysis)等方法。
现有时频分析方法存在问题一：对于非平稳信号频率参数计算存在较大误差。受海森堡不确定性准则(时域和频域的分辨率不能同时达到最高水平)的约束，现有方法不能采用少数几个样本得到准确的局部频率参数(幅值，频率和相位)。当信号为平稳信号时，现有时频分析方法能够给出较为准确的结果。当信号为非平稳信号时，现有时频分析方法给出的分析结果存在误差。一般来说，信号的频率参数变化越快，误差就越大。非平稳信号的时频分析的理论基础是任何非平稳信号可以看成短时间平稳信号的首尾相连，选择较多的样本动摇了这一理论基础，进而造成非平稳信号频率参数估计存在较大的误差。
现有时频分析方法还存在问题二：丢失了幅值和相位变化信息。幅值反映了信号的能量强度，相位反映了信号的位置。以短时傅里叶变换(STFT)为例，该时频图只能近似反映时间和频率的关系，丢失了幅值和相位信息。
发明内容
针对现有技术的至少一个缺陷或改进需求，本发明提供了一种时频分析方法和系统，采用对称离散傅里叶变换进行时频分析，利用了对称离散傅里叶变换的积分特性实现了频率参数计算，还可计算瞬时幅值和瞬时相位，适用于平稳信号和非平稳信号，计算精度高，对信号长度无要求，抗直流干扰能力强，计算时间短。
为实现上述目的，按照本发明的第一方面，提供了一种时频分析方法，包括步骤：
获取输入信号，识别输入信号的信号长度，根据信号长度对输入信号进行第一次信号分解；
对第一次信号分解得到的信号进行第二次信号分解，获得为三个长度相等的信号；
对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；
对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；
分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；
根据三个频域信号的三个虚部求和结果计算信号瞬时频率；
根据三个频域信号的三个实部求和结果及计算得到的瞬时频率计算信号的直流分量、原点相位和幅值。
优选的，所述第一次信号分解包括步骤：
将第一次信号分解前的输入信号记为S(i)，将输入信号S(i)的长度记为K，将信号S(i)分解为K-L+1个长度为L的短信号sj作为第一次信号分解的输出信号，其中j的取值为1,2,3，…，K-L+1；其具体步骤为：选择S(1)至S(L)作为第1个短信号s1，选择选择S(2)至S(L+1)作为第2个短信号s2，选择S(3)至S(L+2) 作为第3个短信号s3，如此依次进行下去，直至获得第K-L+1个短信号sK-L+1；其中L为奇数。
优选的，所述第二次信号分解包括步骤：
将第一次信号分解得到的信号记为sj，信号sj的长度记为L，将信号sj分解为三个长度相等的信号，分别记为x1，x2和x3，其中信号x1为输入序列的第 1个到第L-2k个样本，其中信号x2为输入序列的第1+k个到第L-k个样本，其中信号x3为输入序列的第1+2k个到第L个样本，其中k为正整数，且N为奇数。
优选的，所述离散傅里叶变换包括步骤：
将三个加窗信号分别记为y1，y2和y3，将三个频域信号分别记为Y1，Y2和Y3，加窗信号的长度记为N，频域信号Y1的计算公式为：
Figure BDA0002797772450000031
频域信号Y2的计算公式为：
Figure BDA0002797772450000032
频域信号Y3的计算公式为：
Figure BDA0002797772450000033
其中，m的取值为{m∈Z|-(N-1)/2≤m≤(N-1)/2}。
优选的，所述对三个频域信号的虚部进行求和的计算公式为：
Figure BDA0002797772450000034
Figure BDA0002797772450000035
Figure BDA0002797772450000041
其中Im()表示取虚部运算；Im1为频域信号Y1的虚部求和结果，Im2为频域信号Y2的虚部求和结果，Im3为频域信号Y3的虚部求和结果；
所述对三个频域信号的实部进行求和的计算公式为：
Figure BDA0002797772450000042
Figure BDA0002797772450000043
Figure BDA0002797772450000044
Re1为频域信号Y1的实部求和结果，Re2为频域信号Y2的实部求和结果， Re3为频域信号Y3的实部求和结果。
优选的，所述计算信号相位差和瞬时频率是根据以下方程组得到：
Figure BDA0002797772450000045
Figure BDA0002797772450000046
Figure BDA0002797772450000047
Figure BDA0002797772450000048
为信号x1和信号x2的相位差，
Figure BDA0002797772450000049
为信号sj的原点相位，A为信号sj的幅值，γ是一个是预设的参数，该参数具有平移不变性。
优选的，所述计算信号的直流分量、原点相位和幅值是根据以下方程组得到：
Figure BDA00027977724500000410
Figure BDA00027977724500000411
Figure BDA00027977724500000412
其中，DC表示信号中的直流分量。
优选的，时频分析方法还包括步骤：
在计算原点相位后，对原点相位进行相位解缠绕，根据相位解缠绕后的相位重新计算信号瞬时频率。
优选的，所述相位解缠绕通过计算相位的离散偏导数得到。
按照本发明的第二方面，提供了一种时频分析系统，包括：
第一信号分解模块，用于获取输入信号，识别输入信号的信号长度，对输入信号进行第一次信号分解；
第二信号分解模块，用于对第一次信号分解得到的信号进行第二次信号分解，获得为三个长度相等的信号；
加窗模块，用于对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；
对称离散傅里叶变换模块，用于对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；
求和模块，用于分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；
第一计算模块，用于根据三个频域信号的三个虚部求和结果计算信号瞬时频率；
第二计算模块，用于根据三个频域信号的三个实部求和结果及计算得到的瞬时频率计算信号直流分量、原点相位和幅值。
总体而言，本发明与现有技术相比，具有有益效果：
(1)本发明为一种解析算法，当信号为平稳信号时，可以获得精度极高的解析解，频率的误差小于1*10-12Hz；
(2)当信号为非平稳信号时，虽然直接计算得到的瞬时频率具有较大的误差，但是其瞬时相位具有极高的精度，通过对瞬时相位进行微分重新计算瞬时频率，即间接计算瞬时频率，结果表明重新计算的瞬时频率具有较高的精度。
(3)本发明在计算瞬时频率的同时，还可以计算瞬时幅值和瞬时相位，优于传统的基于时频图的时频信号分析方法。
(4)本发明对信号的长度几乎没有要求，采用五个样本就可以计算得到局部频率值，还可以得到局部幅值、局部相位和直流分量。
(5)本发明为一种解析算法，其利用SDFT变换的虚部不含有直流分量信息计算信号的频率，可以消除直流分量对频率估计的影响，因而本算法具备极强的抗直流干扰能力。
(6)本发明计算时间很短，500次独立仿真试验总共花费19.6秒。
(7)本发明可以用于传感器的矫正，虽然采用减去信号均值的方法可以矫正传感器的零点漂移问题，但是矫正后仍然存在误差。本发明提供了一种计算信号中的直流分量解析解的方法，该方法可以完全消除传感器零点漂移的问题。假设简谐信号的幅值为1个单位，直流分量的计算误差小于1*10-10个单位。
附图说明
图1是本发明实施例的时频分析方法流程图；
图2是本发明实施例的信号第二次分解示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明实施例提供了一种时频分析方法，用于解决现有时频分析方法不能采用少数几个样本计算瞬时频率参数的问题和信号中的直流分量的计算等问题。
参考图1，本发明实施例提供了一种时频分析方法，包括信号第一次分解步骤、信号第二次分解步骤、信号加窗步骤、对称离散傅里叶变换步骤、频域积分求和步骤、计算频率步骤、计算其他频率参数步骤、相位解缠绕步骤、重新计算频率步骤，共九个步骤。后两个步骤并非必须的，在进行平稳信号分析时，采用前七个步骤即可以计算得到信号的频率，相位和幅值，采用本算法得到的频域参数具有极高的精度。在进行非平稳信号分析时，采用前七个步骤可以计算得到信号的瞬时频率，瞬时相位和瞬时幅值，采用第七和第八步骤，可以重新计算瞬时频率，提高瞬时频率的计算精度。各个步骤的优选实现方式如下：
(1)信号第一次分解步骤：
获取输入信号，假设输入信号S(i)的长度为K，将信号S(i)分解为K-L+1 个长度为L的短信号sj；其具体步骤为：选择S(1)至S(L)作为第1个短信号s1，选择选择S(2)至S(L+1)作为第2个短信号s2，选择S(3)至S(L+2)作为第3个短信号s3，如此依次进行下去，直至获得第K-L+1个短信号sK-L+1；本发明算法对信号sj的长度(L)几乎没有要求，仅要求L为大于等于5的奇数。
(2)信号第二次分解步骤：
假设输入信号为sj,信号sj的长度为L；信号第二次分解步骤将输入信号sj分解为三个长度相等的信号(x1，x2和x3)，如图2所示；其中信号x1为输入序列的第1个到第L-2k个样本，其中信号x2为输入序列的第1+k个到第L-k 个样本，其中信号x3为输入序列的第1+2k个到第L个样本；其中k为大于等于1的整数，实际计算时，优选k值取为1。
(3)信号加窗步骤：
窗函数：在信号处理中，窗函数是一种在给定区间之外全部为0的实值函数；给定区间的长度称为窗函数的长度，通常窗函数的长度与信号的长度相等；离散信号分析中，信号的长度为其样本数；通过查询信号分析教材，可以获得多种形式的窗函数；常见的窗函数有矩形窗、汉宁窗、高斯窗、海明窗和三角窗等窗函数。
将上述分解得到的三个信号(x1，x2和x3)依次加窗，依次得到y1，y2和 y3，所谓信号加窗即将信号与窗函数进行Hadamard积，所谓Hadamard积是一个向量乘法运算，即将两长度相等的信号的对应元素相乘，得到一个长度与之前两信号长度相等的信号；窗函数有很多种，本方法不局限于某一种窗函数，优选矩形窗作为本发明实施例算法的默认窗函数；
(4)对称离散傅里叶变换步骤：
连续傅里叶变换(Continuous Fourier transform，CFT)是一个数学变换，其无穷积分导致其现实适用性很差。在数字时代，各行各业中使用的傅里叶变换均为离散傅里叶变换(Discrete Fourier transform，DFT)。当前DFT有两种主要形式，第一种是常规离散傅里叶变换(Ordinary Discrete Fourier Transform， ODFT)，第二种是对称离散傅里叶变换(Symmetric Discrete Fourier Transform， SDFT)。ODFT以其快速傅里叶变换(FastFourier transform，FFT)而被广泛应用，SDFT则应用较少，现有的时频分析方法均基本都是基于ODFT。与ODFT 相比，SDFT更加接近CFT，其具有与CFT类似的对称性和积分特性，而ODFT 不具备这些性质。
本步骤将上述三个加窗信号(y1，y2和y3)依次进行对称离散傅里叶变换 (SDFT)，得到Y1，Y2和Y3；以信号y1为例，信号y1的对称离散傅里叶变换如公式(1)所示，其中m的取值为{m∈Z|-(N-1)/2≤m≤(N-1)/2}；
Figure BDA0002797772450000081
现有离散傅里叶变换(DFT)的快速算法均为快速离散傅里叶变换(FFT)，该算法是常规离散傅里叶变换(ODFT)的快速算法，其并不适用于SDFT；采用SDFT的快速算法可以减少计算时间；理论上，SDFT的快速算法多种多样，不能一一列举；凡是采用的快速算法计算公式(1)应当等同于本步骤。
(5)频域积分求和步骤：
依次对Y1，Y2和Y3的虚部进行积分求和，得到Im1，Im2和Im3；依次对 Y1，Y2和Y3的实部进行积分求和，得到Re1，Re2和Re3；以Y1为例，其虚部积分计算公式如公式(2)所示；其中Im()表示取虚部运算；
Figure BDA0002797772450000082
以Y1为例，其实部积分计算公式如公式(3)所示，
Figure BDA0002797772450000083
频域信号(Y1，Y2和Y3)的虚部奇对称，因而Y1，Y2和Y3的虚部和计算公式分别有三种，以计算Im1为例，另外两种计算公式如公式(4)和(5) 所示；Im2和Im3有类似的计算公式，不一一列举；其中Im()表示取虚部运算；
Figure BDA0002797772450000091
Figure BDA0002797772450000092
(6)计算频率步骤：
根据SDFT频域虚部积分特性，简谐信号的参数有如下三个等式关系：
Figure BDA0002797772450000093
Figure BDA0002797772450000094
Figure BDA0002797772450000095
频域信号(Y1，Y2和Y3)的实部偶对称，因而Y1，Y2和Y3的实部和计算公式分别有三种，以计算Re1为例，另外两种计算公式如公式(9)和(10) 所示；Re2和Re3有类似的计算公式，不一一列举；其中Re()表示取实部运算；
Figure BDA0002797772450000096
Figure BDA0002797772450000097
其中N为三等长信号的长度，
Figure BDA00027977724500000910
为信号x1和信号x2的相位差，
Figure BDA00027977724500000911
为信号 x2的原点相位，A为信号x2的幅值，γ是一个是具备平移不变性的参数；本方法是一种相位差法，利用信号在某时间段内的相位差可以计算得到该信号的频率；信号x1和信号x2的时间间隔等于信号x2和信号x3的时间间隔，两信号之间的时间间隔T为k/fs，其中fs为采样频率；信号x1和信号x2的相位差
Figure BDA00027977724500000912
计算公式如公式(11)所示；
Figure BDA0002797772450000098
根据相位差法的原理，信号的频率计算公式如公式(12)所示；
Figure BDA0002797772450000099
公式(6)-(8)可简化为一个非线性方程组，简化后有且只有三个未知参数，因而相位差
Figure BDA00027977724500000913
有多种计算方式，采用公式(11)计算相位差
Figure BDA00027977724500000914
只是其中的一种较为简单的形式；凡采用公式(6)-(8)构建等式关系获取相位差，即可认定为该算法与本发明算法等同；
(7)计算其他频率参数步骤：
根据SDFT频域实部积分特性，简谐信号的参数有如下三个等式关系：
Figure BDA0002797772450000101
Figure BDA0002797772450000102
Figure BDA0002797772450000103
其中N为三等长信号的长度，DC表示信号中的直流分量，将前述公式(7) 计算得到的相位差
Figure BDA0002797772450000104
带入公式(13)-(15)，可以得到信号sj的直流分量DC，原点相位
Figure BDA0002797772450000105
和幅值A的计算公式如公式(16)-(18)所示；
Figure BDA0002797772450000106
Figure BDA0002797772450000107
Figure BDA0002797772450000108
采用公式(13)-(15)计算直流分量DC，原点相位
Figure BDA00027977724500001013
和幅值A的时候，根据求解顺序的不同，直流分量DC，原点相位
Figure BDA0002797772450000109
和幅值A的计算公式有多种形式；本实施例中先计算直流分量DC，实际上也可以先计算相位或者幅值；凡采用公式(13)-(15)构建等式关系获取信号参数，即可认定为该算法与本发明算法等同；
(8)相位解缠绕步骤：
相位解缠绕：以一个简谐信号
Figure BDA00027977724500001010
为例，其中A为简谐信号的幅值，f为简谐信号的频率，π为圆周率，
Figure BDA00027977724500001011
为信号的初始相位，
Figure BDA00027977724500001012
为时刻t 的相位，理论上相位的取值范围为(-∞，∞)；实际计算时相位被限制在其主值(不管是区间(-π,π]还是[0,2π))，这样的相位称为缠绕的相位；将缠绕的相位还原到区间(-∞，∞)的过程即为相位解缠绕。相位的离散偏导数被广泛应用于相位解缠绕。当简谐信号的相位被限于其主值时，以(-π,π]为例，相位呈周期性变化。由-π到π的过程是一个渐变的过程，再由π突变为-π，然后又渐变为π，反复出现，呈周期性；该变化轮廓明显，层次均匀，突变点为相位周期分界点。因而，理想情况下可以提取相位的离散偏导数，对偏导数进行积分，即可达到相位解缠的目的。
优选的，通过计算相位的离散偏导数来进行相位解缠绕，计算公式如公式 (19)所示；
Figure BDA0002797772450000111
其中i的取值范围(1：K-N)；然后判断ΔΦ(i)的大小，如果ΔΦ(i)小于-π，则
Figure BDA0002797772450000112
至
Figure BDA0002797772450000113
增加2π；
(9)重新计算频率步骤：
通过微分瞬时相位计算瞬时频率(重新计算瞬时频率)的计算公式如公式 (20)所示，
Figure BDA0002797772450000114
其中i的取值范围为(2：K-N)；
本发明实施例的一种时频分析系统，包括：
第一信号分解模块，用于获取输入信号，识别输入信号的信号长度，对输入信号进行第一次信号分解；
第二信号分解模块，用于对第一次信号分解得到的信号进行第二次信号分解，获得为三个长度相等的信号；
加窗模块，用于对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；
对称离散傅里叶变换模块，用于对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；
求和模块，用于分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；
第一计算模块，用于根据三个频域信号的三个虚部求和结果计算信号瞬时频率；
第二计算模块，用于根据三个频域信号的三个实部求和结果及计算得到的瞬时频率计算信号直流分量、原点相位和幅值。
系统的实现原理、技术效果与上述方法相同，此处不再赘述。
必须说明的是，上述任一实施例中，方法并不必然按照序号顺序依次执行，只要从执行逻辑中不能推定必然按某一顺序执行，则意味着可以以其他任何可能的顺序执行。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种时频分析方法，其特征在于，包括步骤：
获取输入信号，识别输入信号的信号长度，对输入信号进行第一次信号分解；
对第一次信号分解得到的信号进行第二次信号分解，获得三个长度相等的信号；
对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；
对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；
分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；
根据三个频域信号的三个虚部求和结果计算信号瞬时频率，所述计算信号瞬时频率是根据以下方程组得到：
Figure FDA0003197052430000011
Figure FDA0003197052430000012
Figure FDA0003197052430000013
Im1、Im2、Im3分别为三个频域信号的虚部求和结果，N为第二次信号分解获得的单个信号的长度，
Figure FDA0003197052430000014
为第二次信号分解获得的第一个信号和第二个信号的相位差，
Figure FDA0003197052430000015
为第一次信号分解得到的信号的原点相位，A为第一次信号分解得到的信号的幅值，γ是一个是预设的参数，该参数具备平移不变性；
根据计算的瞬时频率及三个频域信号的实部和计算信号的直流分量、原点相位和幅值，所述计算信号的直流分量、原点相位和幅值是根据以下方程组得到：
Figure FDA0003197052430000016
Figure FDA0003197052430000017
Figure FDA0003197052430000021
其中，Re1、Re2、Re3分别为三个频域信号的实部求和结果，DC表示第一次信号分解得到的信号的直流分量。
2.如权利要求1所述的一种时频分析方法，其特征在于，所述第一次信号分解包括步骤：
将第一次信号分解前的输入信号记为S(i)，将输入信号S(i)的长度记为K，将信号S(i)分解为K-L+1个长度为L的短信号sj作为第一次信号分解的输出信号，其中j的取值为1,2,3，…，K-L+1；其具体步骤为：选择S(1)至S(L)作为第1个短信号s1，选择S(2)至S(L+1)作为第2个短信号s2，选择S(3)至S(L+2)作为第3个短信号s3，如此依次进行下去，直至获得第K-L+1个短信号sK-L+1；其中L为奇数。
3.如权利要求1所述的一种时频分析方法，其特征在于，所述第二次信号分解包括步骤：
将第一次信号分解得到的信号记为sj，信号sj的长度记为L，将信号sj分解为三个长度为N的信号，分别记为x1，x2和x3，其中信号x1为输入序列的第1个到第L-2k个样本，其中信号x2为输入序列的第1+k个到第L-k个样本，其中信号x3为输入序列的第1+2k个到第L个样本，其中k为正整数，且N为奇数。
4.如权利要求1所述的一种时频分析方法，其特征在于，所述对称离散傅里叶变换包括步骤：
将三个加窗信号分别记为y1，y2和y3，将三个频域信号分别记为Y1，Y2和Y3，加窗信号的长度记为N，频域信号Y1的计算公式为：
Figure FDA0003197052430000022
频域信号Y2的计算公式为：
Figure FDA0003197052430000031
频域信号Y3的计算公式为：
Figure FDA0003197052430000032
其中，m的取值为{m∈Z|-(N-1)/2≤m≤(N-1)/2}。
5.如权利要求4所述的一种时频分析方法，其特征在于，所述对三个频域信号的虚部进行求和的计算公式为：
Figure FDA0003197052430000033
Figure FDA0003197052430000034
Figure FDA0003197052430000035
其中Im()表示取虚部运算；Im1为频域信号Y1的虚部求和结果，Im2为频域信号Y2的虚部求和结果，Im3为频域信号Y3的虚部求和结果；
所述对三个频域信号的实部进行求和的计算公式为：
Figure FDA0003197052430000036
Figure FDA0003197052430000037
Figure FDA0003197052430000038
Re1为频域信号Y1的实部求和结果，Re2为频域信号Y2的实部求和结果，Re3为频域信号Y3的实部求和结果。
6.如权利要求1所述的一种时频分析方法，其特征在于，还包括步骤：
在计算原点相位后，对原点相位进行相位解缠绕，根据相位解缠绕后的相位重新计算信号瞬时频率。
7.如权利要求6所述的一种时频分析方法，其特征在于，所述相位解缠绕通过计算相位的离散偏导数得到。
8.一种时频分析系统，其特征在于，包括：
第一信号分解模块，用于获取输入信号，识别输入信号的信号长度，对输入信号进行第一次信号分解；
第二信号分解模块，用于对第一次信号分解得到的信号进行第二次信号分解，获得为三个长度相等的信号；
加窗模块，用于对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；
对称离散傅里叶变换模块，用于对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；
求和模块，用于分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；
第一计算模块，用于根据三个频域信号的三个虚部求和结果计算信号瞬时频率，所述计算信号瞬时频率是根据以下方程组得到：
Figure FDA0003197052430000041
Figure FDA0003197052430000042
Figure FDA0003197052430000043
Im1、Im2、Im3分别为三个频域信号的虚部求和结果，N为第二次信号分解获得的单个信号的长度，
Figure FDA0003197052430000044
为第二次信号分解获得的第一个信号和第二个信号的相位差，
Figure FDA0003197052430000045
为第一次信号分解得到的信号的原点相位，A为第一次信号分解得到的信号的幅值，γ是一个是预设的参数，该参数具备平移不变性；
第二计算模块，用于根据三个频域信号的三个实部求和结果及计算得到的瞬时频率计算信号直流分量、原点相位和幅值，所述计算信号的直流分量、原点相位和幅值是根据以下方程组得到：
Figure FDA0003197052430000051
Figure FDA0003197052430000052
Figure FDA0003197052430000053
其中，Re1、Re2、Re3分别为三个频域信号的实部求和结果，DC表示第一次信号分解得到的信号的直流分量。
CN202011337810.8A 2020-11-25 2020-11-25 一种时频分析方法和系统 Active CN112505413B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011337810.8A CN112505413B (zh) 2020-11-25 2020-11-25 一种时频分析方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011337810.8A CN112505413B (zh) 2020-11-25 2020-11-25 一种时频分析方法和系统

Publications (2)

Publication Number Publication Date
CN112505413A CN112505413A (zh) 2021-03-16
CN112505413B true CN112505413B (zh) 2021-09-14

Family

ID=74958639

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011337810.8A Active CN112505413B (zh) 2020-11-25 2020-11-25 一种时频分析方法和系统

Country Status (1)

Country Link
CN (1) CN112505413B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113281566B (zh) * 2021-05-11 2023-11-14 重庆矩子兴智能科技有限公司 一种基于组合复信号相位差的频率估计方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5164670A (en) * 1990-09-17 1992-11-17 Syracuse University Multidimensional magnetic resonance system using selective discrete Fourier transformation (SDFT)
US8559568B1 (en) * 2012-01-04 2013-10-15 Audience, Inc. Sliding DFT windowing techniques for monotonically decreasing spectral leakage
CN105991257B (zh) * 2015-01-23 2020-10-23 北京三星通信技术研究有限公司 基于滤波器组的信号生成、发送和接收方法及其装置
CN106053936B (zh) * 2016-06-17 2019-05-14 海南大学 一种获取电学信号瞬时频率的方法及系统
CN109444515B (zh) * 2018-09-26 2021-01-26 徐文涛 一种基于sdft算法的无功、不平衡与谐波检测方法
CN109142867A (zh) * 2018-10-25 2019-01-04 闽南理工学院 基于改进滑窗离散傅里叶变换的谐波检测方法、设备
CN110008434B (zh) * 2019-03-20 2020-11-17 华中科技大学 一种高精度的简谐信号参数估计方法
CN110674456B (zh) * 2019-09-26 2022-11-22 电子科技大学 一种信号采集系统的时频转换方法

Also Published As

Publication number Publication date
CN112505413A (zh) 2021-03-16

Similar Documents

Publication Publication Date Title
Agrez Weighted multipoint interpolated DFT to improve amplitude estimation of multifrequency signal
Duda DFT interpolation algorithm for Kaiser–Bessel and Dolph–Chebyshev windows
Kawahara et al. Tandem-STRAIGHT: A temporally stable power spectral representation for periodic signals and applications to interference-free spectrum, F0, and aperiodicity estimation
Eichstädt et al. Deconvolution filters for the analysis of dynamic measurement processes: a tutorial
CN112505413B (zh) 一种时频分析方法和系统
Carini et al. Nonlinear system identification using Wiener basis functions and multiple-variance perfect sequences
Petri Frequency-domain testing of waveform digitizers
CN112162152A (zh) 基于相位直线拟合的正弦波相参脉冲串信号频率估计方法
EP3396670B1 (en) Speech signal processing
Dragotti et al. Exact sampling results for signals with finite rate of innovation using Strang-Fix conditions and local kernels
Rajamani et al. A novel method for designing allpass digital filters
JP2812184B2 (ja) 音声の複素ケプストラム分析装置
Orović et al. A class of highly concentrated time-frequency distributions based on the ambiguity domain representation and complex-lag moment
CN112652321B (zh) 一种基于深度学习相位更加友好的语音降噪系统及方法
CN111289800B (zh) 一种基于广义回归神经网络的小电阻振动监测方法
Petrović et al. New procedure for harmonics estimation based on Hilbert transformation
Funaki Sparse Time-Varying Complex AR (TV-CAR) speech analysis based on Adaptive LASSO
Todoran et al. Discrete hilbert transform. numeric algorithms
Shavelis et al. Signal sampling according to time-varying bandwidth
Zhang et al. Determining the length of sliding window by using frequency decomposition
Burtea et al. Estimating the frequencies of vibration signals using a machine learning algorithm with explained predictions
Li Half-infinite sampling and its FT
Wolf et al. Amplitude and frequency estimator for aperiodic multi-frequency noisy vibration signals of a tram gearbox
Serov et al. Estimation of the Signal Parameters Measurement Error for the Case of ADC Nonlinearity Approximation by Chebyshev Polynomial
Goos et al. Continuous time frequency domain LPV state space identification via periodic time-varying input-output modeling

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