CN107608935B - 基于时间重排压缩变换的冲击类信号时频分析与重构方法 - Google Patents

Abstract

本发明公开了一种基于时间重排压缩变换的冲击类信号时频分析与重构方法，包括：1)计算待分析离散信号的短时傅里叶变换，得到相应的时频复矩阵S_x[n,k]；2)对时间变量与信号的乘积作短时傅里叶变换，得到时频复矩阵S_tx[n,k]；3)计算群延时估计算子；4)对步骤1中得到的时频复矩阵S_x[n,k]仅沿时间方向重新排列，得到时间重排压缩变换后的时频矩阵V_x[m,k]；5)将时间重排压缩变换得到的时频矩阵V_x[m,k]的每一行元素相加，得到一维列向量，接着再对该列向量除以短时傅里叶变换所用窗函数的均值得到重构信号的频谱。本发明涉及的时间重排压缩变换对于分析冲击类信号，得到的时频图的聚集性比同步压缩变换更高，具有良好的抗噪性能。与传统时频重排相比，具有可重构，计算速度快的优点。

Description

基于时间重排压缩变换的冲击类信号时频分析与重构方法

技术领域

本发明属于非平稳信号的时频分析领域，特别涉及一种基于时间重排压缩变换的冲击类信号时频分析与重构方法。

背景技术

时频分析方法是将一维的时域信号变换为二维的联合时间-频率域(时频域)分布的方法，即综合考虑了时间和频率两个变量，能够揭示信号中各成分随时间或频率的变化规律。由于时频分析的时-频联合表征能力，使其广泛应用于包括振动信号分析、语音信号处理、图像处理等在内的多个领域。时频分析方法的研究目标主要有两个，一是提高可读性，二是提高可重构性。高可读性指的是时频图具有较高的时频聚集性，且没有交叉项的干扰，能够容易地从时频图中得到所需信息；可重构性指的是对时频图能够进行逆变换，将二维分布重构回到一维信号。所以衡量一个时频分析方法的优劣时，可以从可读性和可重构性这两个方面进行分析。

一般来讲，时频分析方法可以分为线性时频分析方法和二次时频分析方法两大类。短时傅里叶变换(STFT)、连续小波变换(CWT)属于线性时频分析方法。由于Heisenberg测不准原理的限制，线性时频分析方法无法得到很高的时频聚集性，可读性不高；二次时频分析方法有Wigner-Ville分布(WVD)、Cohen类分布等等，它们具有比线性时频分析方法更高的时频聚集性。但是二次时频分析方法会受到交叉项的干扰，同样会影响时频图的可读性。并且二次时频分析方法无法进行线性叠加，重构时域信号也较为困难。1995年学者Auger与Flandrin提出了一种称为时频重排的时频后处理方法，通过把时频分布系数从原先计算的位置，沿着时间轴和频率轴两个方向重新移动至能量分布重心位置，提高时频变换结果的聚集性。然而这种时频重排方法无法将二维时频分布逆变换为一维信号，即无法重构信号，而重构性质恰恰是某些应用领域对时频分析方法所提出的关键要求。2011年Daubechies等学者提出了另一种时频后处理方法，称为同步压缩变换。它通过对时频变换的结果仅沿频率方向重新排列，将每个时频点的能量叠加至能量中心，比传统时频分析(如短时傅里叶变换，小波变换)具有更高的时频聚集性。另外，同步压缩变换在两个方面优于传统时频重排，一是同步压缩变换可以对时频图进行逆变换，得到时域信号，具有良好的可重构性；二是同步压缩变换是一维积分，运算速度比二维积分的传统时频重排快。

然而，现有的同步压缩变换方法大多适合提取谐波类信号，其在时频图上大致以“横线”的形式出现，如图1(a)所示，而对于图1(b)所示的以“竖线”形式出现的冲击类信号，提取效果不明显。原因是对于冲击类信号，其瞬时频率变化率很大，在时频图上表现为时频脊线的斜率很大。而理论分析表明，对于这种情况，如果按照现有同步压缩变换的方法进行时频分析，会存在较大的瞬时频率估计误差，导致时频聚集性很差，并且重构误差也较大。因此，为了提高时频分析与重构的精度，需要对现有的同步压缩变换进行改进，使其能够准确提取冲击类信号特征。

发明内容

本发明的目的是提供一种基于时间重排压缩变换的冲击类信号时频分析与重构方法，以解决上述问题。

为达到以上目的，本发明采取如下方案：

基于时间重排压缩变换的冲击类信号时频分析与重构方法，包括以下步骤：

(1)计算待分析离散信号x[n]的短时傅里叶变换，得到相应的时频复矩阵S_x[n,k]；

(2)对时间变量与待分析信号的乘积作短时傅里叶变换，得到时频复矩阵S_tx[n,k]；

(3)将步骤(2)中得到的时频复矩阵S_tx[n,k]的各元素除以步骤(1)中得到的时频复矩阵S_x[n,k]的对应元素，再对结果取实部，得到群延时估计算子矩阵；

(4)利用步骤(3)中得到的群延时估计算子矩阵，对步骤(1)中得到的待分析信号的时频复矩阵S_x[n,k]仅沿时间方向重新排列，将每个时频点系数值叠加至群延时估计算子矩阵中的元素所指向的时间中心，得到时间重排压缩变换后的时频矩阵V_x[m,k]；

(5)从时频域重构回频域：将时间重排压缩变换得到的时频矩阵V_x[m,k]的每一行元素相加，得到一维列向量，接着再对该列向量除以短时傅里叶变换所用窗函数的均值得到重构信号的频谱。

进一步的，还包括以下步骤：(6)对重构信号的频谱进行傅里叶反变换，计算重构后的时域信号。

进一步的，步骤(2)中，计算离散后的时间变量t[n]与离散后的信号x[n]的乘积t[n]·x[n]，将其作为一个新的信号，对该信号作离散短时傅里叶变换

得到时频复矩阵S_tx[n,k]；

进一步的，步骤(3)中，按照公式

将步骤(2)中得到的时频复矩阵S_tx[n,k]的各元素除以步骤(1)中得到的时频复矩阵S_x[n,k]的对应元素，再对结果取实部，得到群延时估计算子矩阵

进一步的，步骤(4)中，利用步骤(3)得到的群延时估计算子

对原始信号的短时傅里叶变换时频矩阵S_x[n,k]沿时间方向进行重排；将每个时频点[n,k]处的系数值S_x[n,k]仅沿时间方向移动到与群延时估计算子值

距离最近的时间元素t[m]处，即按照公式

重新排列相加；当遍历完S_x[n,k]所有的时频点，完成重排运算，实现对时频图S_x[n,k]时间方向的压缩。

进一步的，步骤(5)中，将经过时间重排压缩变换得到的时频矩阵V_x[m,k]的每一行元素相加，得到一维列向量

接着再对该列向量除以短时傅里叶变换所用窗函数的均值

得到重构信号的频谱

横坐标是频率序列，纵坐标是频谱复值，包含幅值和相位两个信息。

进一步的，步骤(6)利用步骤(5)中计算的重构信号的频谱

进行离散傅里叶反变换

得到重构后的时域信号x_r[n]。

进一步的，步骤(1)待分析离散信号x[n]的长度为N点，n＝0,1,…,N-1,采样时间间隔为T，时间元素对应的坐标为t[n]＝nT，频率元素对应的坐标为f[k]＝k/NT,k＝0,1,…,N-1；短时傅里叶变换的时频复矩阵S_x[n,k]通过下面的公式

计算，其中i表示虚数单位，下同。

再进一步，步骤(3)中，考虑到实际信号总是会被噪声干扰，此时只有足够大的短时傅里叶变换模系数，才不会引起因复数相除而导致群延时估计算子计算不稳定的问题。因此，实际计算时仅利用模大于阈值γ的那些时频系数来计算群延时估计算子

Ξ_[n,k]＝{n,k∈N:|S_x[n,k]|>γ}

式中，阈值γ>0用于克服S_x[n,k]某些时频系数模太小时导致的群延时估计算子计算不稳定问题，阈值γ常取为10^-6。

相对于现有技术，本发明具有以下有益效果：

1、与短时傅里叶变换和基于短时傅里叶变换的同步压缩变换相比，本发明涉及的时间重排压缩变换针对分析冲击类信号，得到的时频聚集性更高，时频图具有更强的可读性。另外，对于冲击类信号的时频表示，时间重排压缩变换具有更强的抗噪性能。

2、与传统时频重排相比，时间重排压缩变换的逆变换可以将二维时频分布重构到一维频域，进而可以通过傅里叶反变换得到重构后的时域信号，因此具有传统时频重排不具有的重构特性。并且由于时间重排压缩变换是一维积分，计算时间复杂度比传统时频重排低。

本发明一种基于时间重排压缩变换的冲击类信号时频分析与重构方法，仅利用群延时估计算子，用它作为时间方向的指针，对时频变换的结果仅在时间方向上进行重排，实现冲击类信号特征提取。并且由于时间重排压缩变换是时间方向上的一维积分操作，因此时间重排压缩变换的逆变换可以将二维时频分布重构回一维频域，具有传统时频重排不具有的重构特性，并且运算速度比传统时频重排快。

附图说明

图1为背景技术中的两类信号时频形态；其中，图1(a)是谐波类信号的时频特征形态，图1(b)是冲击类信号的时频特征形态。图1(a)中横坐标表示时间，单位为s；纵坐标表示角频率，单位为rad/s；图1(b)中横坐标表示时间，单位为s；纵坐标表示角频率，单位为rad/s。

图2为本发明基于时间重排压缩变换的冲击类信号时频分析与重构方法的流程图。

图3为含噪声的冲击类仿真信号时域波形；图中横坐标表示时间，单位为s；纵坐标表示信号幅值，单位为mm。

图4是仿真信号的短时傅里叶变换时频图；其中横坐标表示时间，单位为s；纵坐标表示频率，单位为Hz。

图5是仿真信号的时间重排压缩变换时频图与同步压缩变换时频图的对比；其中，图5(a)是仿真信号的时间重排压缩变换时频图；图5(b)是仿真信号的同步压缩变换时频图；图5(a)中横坐标表示时间，单位为s；纵坐标表示频率，单位为Hz；图5(b)中横坐标表示时间，单位为s；纵坐标表示频率，单位为Hz。

图6是时间重排压缩变换的重构特性说明示意图；图6(a)对应原始的仿真信号；而图6(b)对应经过时间重排压缩变换并重构后得到的时域信号；图6(a)中横坐标表示时间，单位为s；纵坐标表示信号幅值，单位为mm；图6(b)中横坐标表示时间，单位为s；纵坐标表示信号幅值，单位为mm。

具体实施方式

下面结合一个实例来验证本发明在应用中的有效性，但本实例并不用于限制本发明。

为了说明时间重排压缩变换相对于传统短时傅里叶变换以及同步压缩变换的优越性，选取一个含噪声的仿真信号，其时域波形图如图3所示，可以看到仿真信号在时域上持续时间较短，表现为一定的冲击特性。仿真信号的表达式为：

式中相位

f∈[0,1MHz]，频率间隔为1KHz，即f＝0:1000:1×10⁶Hz。F^-1(·)表示傅里叶反变换；noise表示添加高斯白噪声。信号采样总时长为1ms，采样频率为2MHz，采样点数N＝2000，添加高斯白噪声使得信噪比为5dB。从表达式中看出仿真信号属于一种非平稳信号，由两个分量组成，两部分的群延时都随着频率的变化而变化，根据公式

可计算该信号两个分量的理论群延时分别为

请参阅图2所示，本发明一种基于时间重排压缩变换的冲击类信号时频分析与重构方法，包括下述步骤：

(1)用采样频率2MHz离散化仿真信号x(t)，得到离散信号x[n]。计算x[n]的离散短时傅里叶变换，窗函数取高斯窗，宽度参数σ＝0.01ms，窗函数沿时间方向每次移动一个点，得到大小为1001×2000的时频复矩阵S_x[n,k]，并将该矩阵中模小于1×10^-6的元素置为NAN。将短时傅里叶变换时频矩阵结合时间轴t和频率轴f，绘制成时频图，如图4所示。从图中可以看到仿真信号的短时傅里叶变换时频图能够大致反映群延时随频率的变化规律，然而它的时频聚集性和信噪比都不高，因此需要进一步的后处理操作。

(2)计算时间变量与离散信号的乘积t[n]·x[n]，并将其作为一个新的信号，作短时傅里叶变换

得到与步骤(1)中同样尺寸的时频复矩阵S_tx[n,k]。

(3)按照公式

将步骤(2)中得到的时频复矩阵S_tx[n,k]的各元素除以步骤(1)中得到的时频复矩阵S_x[n,k]的对应元素，再对结果取实部，得到群延时估计算子矩阵

(4)利用步骤(3)中得到的群延时估计算子矩阵

对步骤(1)得到的原始信号的短时傅里叶变换时频矩阵S_x[n,k]沿时间方向进行重排。将每个时频点[n,k]处的系数值S_x[n,k]仅沿时间方向移动到与对应的群延时估计算子矩阵元素值

距离最近的时间元素t[m]处，按照公式

重新排列，进行复数相加。当遍历完S_x[n,k]所有的时频点，则完成了重排运算，实现对时频图S_x[n,k]时间方向的压缩，得到如图5(a)所示的时间重排压缩变换时频图。

将得到图5(a)所示的时间重排压缩变换时频图的与图5(b)所示的同步压缩变换时频图进行比较，可以发现时间重排压缩变换对于分析冲击类信号，其时频聚集性优于同步压缩变换。并且由于时间重排压缩变换仅在时间方向重排，是一维积分操作，运算速度理论上应与同步压缩变换相当。经过实际测算，本发明涉及的时间重排压缩变换运算时间为1.4s，同步压缩变换运算时间为1.6s，两者相当。然而相比起来，时间和频率方向均进行重排的传统时频重排运算时间达到了3.2s，运算复杂度较大。另外，通过比较图5(a)和图5(b)，可以看到在分析含噪声冲击类信号的情况下，时间重排压缩变换时频图中噪声分布区域较少，与同步压缩变换相比具有更好的抗噪性能。

(5)将时间重排压缩变换得到的时频矩阵V_x[m,k]的每一行元素相加，得到一维列向量

接着再对该列向量除以短时傅里叶变换所用窗函数的均值

即可得到重构信号的频谱

其横坐标是频率序列，纵坐标是频谱复值，包含幅值和相位两个信息。

(6)可选择的，利用步骤(5)中计算的重构信号的频谱

进行离散傅里叶反变换

得到重构后的时域信号x_r[n]。如图6所示。其中图6(a)是原始信号时域波形，图6(b)是重构后的时域信号波形，可以看到两者几乎完全相同。为了定量评价时间重排压缩变换的重构性能，引入原信号x[n]与重构信号x_r[n]之间的均方根误差

作为指标，计算得到RMSE为8×10^-4，远小于信号幅值的数量级，表明时间重排压缩变换能够准确地重构冲击类信号。

Claims (6)

1.基于时间重排压缩变换的冲击类信号时频分析与重构方法，其特征在于，包括以下步骤：

(1)计算待分析离散信号x[n]的短时傅里叶变换，得到相应的时频复矩阵S_x[n,k]；

(2)对时间变量与待分析信号的乘积作短时傅里叶变换，得到时频复矩阵S_tx[n,k]；

(3)计算群延时估计算子矩阵：将步骤(2)中得到的时频复矩阵S_tx[n,k]的各元素除以步骤(1)中得到的时频复矩阵S_x[n,k]的对应元素，再对结果取实部，得到群延时估计算子矩阵；

(4)利用步骤(3)中得到的群延时估计算子矩阵，对步骤(1)中得到的待分析信号的时频复矩阵S_x[n,k]仅沿时间方向重新排列，将每个时频点系数值叠加至群延时估计算子矩阵中的元素所指向的时间中心，得到时间重排压缩变换后的时频矩阵V_x[m,k]；

(5)从时频域重构回频域：将时间重排压缩变换得到的时频矩阵V_x[m,k]的每一行元素相加，得到一维列向量，再对该列向量除以短时傅里叶变换所用窗函数的均值得到重构信号的频谱；

步骤(1)待分析离散信号x[n]的长度为N点，n＝0,1,…,N-1,采样时间间隔为T，时间元素对应的坐标为t[n]＝nT，频率元素对应的坐标为f[k]＝k/NT,k＝0,1,…,N-1；短时傅里叶变换的时频复矩阵S_x[n,k]通过公式

计算，其中i表示虚数单位；

步骤(3)中，按照公式

将步骤(1)和步骤(2)中分别得到的S_x[n,k]和S_tx[n,k]复矩阵中对应元素相除，并对结果取实部，得到群延时估计算子矩阵

步骤(4)中，利用步骤(3)得到的群延时估计算子

对原始信号的短时傅里叶变换时频矩阵S_x[n,k]沿时间方向进行重排；将每个时频点[n,k]处的系数值S_x[n,k]仅沿时间方向移动到与群延时估计算子值

距离最近的时间元素t[m]处，即按照公式

重新排列相加；当遍历完S_x[n,k]所有的时频点，完成重排运算，实现时间方向的压缩。

2.根据权利要求1所述的基于时间重排压缩变换的冲击类信号时频分析与重构方法，其特征在于，还包括以下步骤：

(6)对重构信号的频谱进行傅里叶反变换，计算重构后的时域信号。

3.根据权利要求1所述的基于时间重排压缩变换的冲击类信号时频分析与重构方法，其特征在于，步骤(2)中，计算离散后的时间变量t[n]与离散后的信号x[n]的乘积t[n]·x[n]，将其作为一个新的信号，对该信号作离散短时傅里叶变换

得到时频复矩阵S_tx[n,k]。

4.根据权利要求1所述的基于时间重排压缩变换的冲击类信号时频分析与重构方法，其特征在于，步骤(5)中，将经过时间重排压缩变换得到的时频矩阵V_x[m,k]的每一行元素相加，得到一维列向量

接着再对该列向量除以短时傅里叶变换所用窗函数的均值

得到重构信号的频谱

横坐标是频率序列，纵坐标是频谱复值，包含幅值和相位两个信息。

5.根据权利要求2所述的基于时间重排压缩变换的冲击类信号时频分析与重构方法，其特征在于，步骤(6)利用步骤(5)中计算的重构信号的频谱

进行离散傅里叶反变换

得到重构后的时域信号x_r[n]。

6.根据权利要求1所述的基于时间重排压缩变换的冲击类信号时频分析与重构方法，其特征在于，步骤(3)中，为保证数值计算的稳定性，计算时仅利用模大于阈值γ的时频系数来计算群延时估计算子矩阵

Ξ_[n,k]＝{n,k∈N:|S_x[n,k]|>γ}

式中，阈值γ取为10^-6。

Images