一种时频分析方法和系统
阅读说明：本技术 一种时频分析方法和系统 (Time-frequency analysis method and system ) 是由 轩建平 李锐 唐律 张青 于 2020-11-25 设计创作，主要内容包括：本发明公开了一种时频分析方法和系统。该方法包括步骤：对输入信号进行第一次信号分解；对第一次信号分解得到的信号进行第二次信号分解,获得三个长度相等的信号；对第二次信号分解得到的三个信号分别进行加窗,得到三个加窗信号；对三个加窗信号分别进行对称离散傅里叶变换,得到三个频域信号；分别对三个频域信号的虚部进行求和,分别对三个频域信号的实部进行求和；根据虚部求和结果计算信号瞬时频率；根据计算的瞬时频率及实部求和结果计算信号的直流分量、原点相位和幅值。本发明实现了瞬时频率、直流分量、瞬时幅值和瞬时相位的计算,适用于平稳信号和非平稳信号,计算精度高,对信号长度几乎无要求,抗直流干扰能力强,计算时间短。(The invention discloses a time-frequency analysis method and a time-frequency analysis system. The method comprises the following steps: performing first signal decomposition on an input signal; performing second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length; windowing three signals obtained by the second signal decomposition respectively to obtain three windowed signals; symmetric discrete Fourier transform is respectively carried out on the three windowed signals to obtain three frequency domain signals; summing the imaginary parts of the three frequency domain signals respectively, and summing the real parts of the three frequency domain signals respectively; calculating the instantaneous frequency of the signal according to the result of the summation of the imaginary parts; and calculating the direct current component, the origin phase and the amplitude of the signal according to the calculated instantaneous frequency and the real part summation result. The invention realizes the calculation of instantaneous frequency, direct current component, instantaneous amplitude and instantaneous phase, is suitable for stable signals and non-stable signals, has high calculation precision, almost no requirement on signal length, strong direct current interference resistance and short calculation time.)
技术领域
本发明属于信号分析技术领域，更具体地，涉及一种时频分析方法和系统。
背景技术
信号是信息的载体，人们认识世界万事万物及通讯均基于对信息的理解。获取信号中的信息依赖于信号分析。通常实际信号均为时变信号，科学家们多采用非平稳信号代指时变信号。时频分析被广泛应用于分析非平稳信号，如分析引力波，动物的语言，脑电波等信号。随着人类认识世界的深入，对时频分析方法的要求越来越高。
现有技术中存在多种时频分析方法。例如线性时频描述法(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的计算公式为：
频域信号Y2的计算公式为：
频域信号Y3的计算公式为：
其中，m的取值为{m∈Z|-(N-1)/2≤m≤(N-1)/2}。
优选的，所述对三个频域信号的虚部进行求和的计算公式为：
其中Im()表示取虚部运算；Im1为频域信号Y1的虚部求和结果，Im2为频域信号Y2的虚部求和结果，Im3为频域信号Y3的虚部求和结果；
所述对三个频域信号的实部进行求和的计算公式为：
Re1为频域信号Y1的实部求和结果，Re2为频域信号Y2的实部求和结果，Re3为频域信号Y3的实部求和结果。
优选的，所述计算信号相位差和瞬时频率是根据以下方程组得到：
为信号x1和信号x2的相位差，为信号sj的原点相位，A为信号sj的幅值，γ是一个是预设的参数，该参数具有平移不变性。
优选的，所述计算信号的直流分量、原点相位和幅值是根据以下方程组得到：
其中，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}；
现有离散傅里叶变换(DFT)的快速算法均为快速离散傅里叶变换(FFT)，该算法是常规离散傅里叶变换(ODFT)的快速算法，其并不适用于SDFT；采用SDFT的快速算法可以减少计算时间；理论上，SDFT的快速算法多种多样，不能一一列举；凡是采用的快速算法计算公式(1)应当等同于本步骤。
(5)频域积分求和步骤：
依次对Y1，Y2和Y3的虚部进行积分求和，得到Im1，Im2和Im3；依次对Y1，Y2和Y3的实部进行积分求和，得到Re1，Re2和Re3；以Y1为例，其虚部积分计算公式如公式(2)所示；其中Im()表示取虚部运算；
以Y1为例，其实部积分计算公式如公式(3)所示，
频域信号(Y1，Y2和Y3)的虚部奇对称，因而Y1，Y2和Y3的虚部和计算公式分别有三种，以计算Im1为例，另外两种计算公式如公式(4)和(5)所示；Im2和Im3有类似的计算公式，不一一列举；其中Im()表示取虚部运算；
(6)计算频率步骤：
根据SDFT频域虚部积分特性，简谐信号的参数有如下三个等式关系：
频域信号(Y1，Y2和Y3)的实部偶对称，因而Y1，Y2和Y3的实部和计算公式分别有三种，以计算Re1为例，另外两种计算公式如公式(9)和(10)所示；Re2和Re3有类似的计算公式，不一一列举；其中Re()表示取实部运算；
其中N为三等长信号的长度，为信号x1和信号x2的相位差，为信号x2的原点相位，A为信号x2的幅值，γ是一个是具备平移不变性的参数；本方法是一种相位差法，利用信号在某时间段内的相位差可以计算得到该信号的频率；信号x1和信号x2的时间间隔等于信号x2和信号x3的时间间隔，两信号之间的时间间隔T为k/fs，其中fs为采样频率；信号x1和信号x2的相位差计算公式如公式(11)所示；
根据相位差法的原理，信号的频率计算公式如公式(12)所示；
公式(6)-(8)可简化为一个非线性方程组，简化后有且只有三个未知参数，因而相位差有多种计算方式，采用公式(11)计算相位差只是其中的一种较为简单的形式；凡采用公式(6)-(8)构建等式关系获取相位差，即可认定为该算法与本发明算法等同；
(7)计算其他频率参数步骤：
根据SDFT频域实部积分特性，简谐信号的参数有如下三个等式关系：
其中N为三等长信号的长度，DC表示信号中的直流分量，将前述公式(7)计算得到的相位差带入公式(13)-(15)，可以得到信号sj的直流分量DC，原点相位和幅值A的计算公式如公式(16)-(18)所示；
采用公式(13)-(15)计算直流分量DC，原点相位和幅值A的时候，根据求解顺序的不同，直流分量DC，原点相位和幅值A的计算公式有多种形式；本实施例中先计算直流分量DC，实际上也可以先计算相位或者幅值；凡采用公式(13)-(15)构建等式关系获取信号参数，即可认定为该算法与本发明算法等同；
(8)相位解缠绕步骤：
相位解缠绕：以一个简谐信号为例，其中A为简谐信号的幅值，f为简谐信号的频率，π为圆周率，为信号的初始相位，为时刻t的相位，理论上相位的取值范围为(-∞，∞)；实际计算时相位被限制在其主值(不管是区间(-π,π]还是[0,2π))，这样的相位称为缠绕的相位；将缠绕的相位还原到区间(-∞，∞)的过程即为相位解缠绕。相位的离散偏导数被广泛应用于相位解缠绕。当简谐信号的相位被限于其主值时，以(-π,π]为例，相位呈周期性变化。由-π到π的过程是一个渐变的过程，再由π突变为-π，然后又渐变为π，反复出现，呈周期性；该变化轮廓明显，层次均匀，突变点为相位周期分界点。因而，理想情况下可以提取相位的离散偏导数，对偏导数进行积分，即可达到相位解缠的目的。
优选的，通过计算相位的离散偏导数来进行相位解缠绕，计算公式如公式(19)所示；
其中i的取值范围(1：K-N)；然后判断ΔΦ(i)的大小，如果ΔΦ(i)小于-π，则至增加2π；
(9)重新计算频率步骤：
通过微分瞬时相位计算瞬时频率(重新计算瞬时频率)的计算公式如公式(20)所示，
其中i的取值范围为(2：K-N)；
本发明实施例的一种时频分析系统，包括：
第一信号分解模块，用于获取输入信号，识别输入信号的信号长度，对输入信号进行第一次信号分解；
第二信号分解模块，用于对第一次信号分解得到的信号进行第二次信号分解，获得为三个长度相等的信号；
加窗模块，用于对第二次信号分解得到的三个信号分别进行加窗，得到三个加窗信号；
对称离散傅里叶变换模块，用于对三个加窗信号分别进行对称离散傅里叶变换，得到三个频域信号；
求和模块，用于分别对三个频域信号的虚部进行求和，分别对三个频域信号的实部进行求和；
第一计算模块，用于根据三个频域信号的三个虚部求和结果计算信号瞬时频率；
第二计算模块，用于根据三个频域信号的三个实部求和结果及计算得到的瞬时频率计算信号直流分量、原点相位和幅值。
系统的实现原理、技术效果与上述方法相同，此处不再赘述。
必须说明的是，上述任一实施例中，方法并不必然按照序号顺序依次执行，只要从执行逻辑中不能推定必然按某一顺序执行，则意味着可以以其他任何可能的顺序执行。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。