CN112883326B - 一种基于流式算法的自适应时频变换方法 - Google Patents

Abstract

本发明适用于信号处理领域，提供了一种基于流式算法的自适应时频变换方法，其中，所述自适应时频变换方法包括：对实数离散信号进行滑动时窗处理，得分段信号；对所述分段信号进行傅里叶级数展开并引入时变傅里叶系数，以将所述分段信号表示为包含时变傅里叶系数的形式；根据所述包含时变傅里叶系数的形式，利用流式算法自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱。本发明解析计算正则化条件下的数学欠定问题，避免了局部时频变换中的迭代算法，有效的节约了计算过程中的内存成本；同时，本发明参数选取灵活，可以自主选取频率采样间隔和频率范围，只保留有效的频率信息，节省了数据存储空间。

Description

一种基于流式算法的自适应时频变换方法

技术领域

本发明属于信号处理领域，尤其涉及一种基于流式算法的自适应时频变换方法。

背景技术

许多技术领域的信号大都是非平稳信号，如地球物理勘探、医学和通信等。传统的傅里叶变换只能表征信号的全局频谱特征，不适于处理和分析非平稳信号。时频分析研究了信号的频率信息随时间的变化规律，能更全面、细致地表征信号的特征，对非平稳信号的处理与解释具有重要意义。

20世纪40年代，Gabor提出了短时傅里叶变换(STFT)，该变换通过对信号进行加窗处理实现分段傅里叶变换，可以表征非平稳信号的时频特征，但是固定的窗函数会限制STFT的时频分辨率。小波变换通过伸缩平移基本小波对信号进行多尺度分析，实现了信号的多分辨率分析，但是小波变换存在尺度频率转换和小波基选取等问题。Stockwell等于1996年提出了S变换，它采用时窗长度与频率成反比的高斯函数，具有可变的时频分辨率。Liu等提出局部时频变换(LTFT)方法，该方法定义时变傅里叶系数来表征时频谱，且利用整形正则化条件求解关于时变傅里叶系数的最小二乘解的欠定问题。LTFT参数选取较灵活，通过平滑半径参数约束时频分辨率，且能够自主选取频率采样间隔和频率范围，但是LTFT中的迭代算法会导致运算速度慢、内存成本高等问题。

发明内容

本发明实施例的目的在于提供一种基于流式算法的自适应时频变换方法，旨在解决现有的局部时频变换方法存在运算复杂、内存成本高等问题。

本发明实施例是这样实现的，一种基于流式算法的自适应时频变换方法，包括：

对实数离散信号进行滑动时窗处理，得分段信号；

对所述分段信号进行傅里叶级数展开并引入时变傅里叶系数，以将所述分段信号表示为包含时变傅里叶系数的形式；

根据所述包含时变傅里叶系数的形式，利用流式算法自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱。

本发明实施例提供的基于流式算法的自适应时频变换方法，通过对实数离散信号进行滑动时窗处理，利用傅里叶级数匹配目标信号并引入时变傅里叶系数，将傅里叶分析与非平稳回归分析相结合，通过流式计算理论使时变傅里叶系数随着数据的改变而自适应更新，可以有效的表征非平稳信号的时频特征。本发明解析计算正则化条件下的数学欠定问题，避免了局部时频变换中的迭代算法，有效的节约了计算过程中的内存成本；同时，本发明参数选取灵活，可以自主选取频率采样间隔和频率范围，只保留有效的频率信息，节省了数据存储空间。

附图说明

图1为本发明实施例提供的一种基于流式算法的自适应时频变换方法的实现流程图；

图2为本发明实施例提供的另一种基于流式算法的自适应时频变换方法的实现流程图；

图3为本发明实施例提供的具有两个交叉线性调频信号的合成信号；

图4为本发明实施例提供的线性调频信号的理论频率；

图5为本发明实施例提供的经过本发明时频变换方法处理后的时频谱图；

图6为本发明实施例提供的经过本发明时频变换方法反变换后的信号图；

图7为本发明实施例提供的反变换后的信号与原始合成信号的差值图；

图8为本发明实施例提供的单道实际地震数据图；

图9为本发明实施例提供的经过本发明时频变换方法处理所得时频谱图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

可以理解，本申请所使用的术语“第一”、“第二”等可在本文中用于描述各种元件，但除非特别说明，这些元件不受这些术语限制。这些术语仅用于将第一个元件与另一个元件区分。举例来说，在不脱离本申请的范围的情况下，可以将第一xx脚本称为第二xx脚本，且类似地，可将第二xx脚本称为第一xx脚本。

本发明实施例为了解决现有的局部时频变换方法存在运算复杂、内存成本高等问题，在短时傅里叶变换和局部时频变换的基础上，引入流式计算的概念，提出了基于流式算法的自适应时频变换方法，通过对实数离散信号进行滑动时窗处理，利用傅里叶级数匹配目标信号并引入时变傅里叶系数，将傅里叶分析与非平稳回归分析相结合，通过流式计算理论使时变傅里叶系数随着数据的改变而自适应更新，可以有效的表征非平稳信号的时频特征。本发明解析计算正则化条件下的数学欠定问题，避免了局部时频变换中的迭代算法，有效的节约了计算过程中的内存成本；同时，本发明参数选取灵活，可以自主选取频率采样间隔和频率范围，只保留有效的频率信息，节省了数据存储空间。

为了进一步阐述本发明为实现预定发明目的所采取的技术手段及功效，以下结合附图及较佳实施例，对依据本发明的具体实施方式、结构、特征及其功效，详细说明如下。

图1示出了本发明实施例提供的一种基于流式算法的自适应时频变换方法的实现流程，为了便于说明，仅示出与本发明实施例相关的部分，详述如下：

步骤S101中，对实数离散信号进行滑动时窗处理，得分段信号。

在本发明实施例中，设x(n)(n＝0,1,…,N-1)为实数离散信号，n为信号的采样点数，N为信号的总长度。为了得到信号的局部时频谱，首先通过滑动时窗处理将信号分成短时间片段，分段信号可表示为：

x_l(m)＝x(m+l-h)w(m) (1)

其中，m＝0,1,…,M-1为分段信号的采样点数，M为时窗的长度。若M为偶数，h＝M/2；若M为奇数，h＝(M-1)/2。x_l(m)(l＝0,1,…,N-1)代表第l个分段信号。w(m)代表窗函数，若直接利用原始分段信号(即采用矩形窗函数)进行时频分析，会造成较大的频谱泄露，所以为了减小频谱泄露选用汉宁窗对信号进行处理，汉宁窗函数表达式为：

窗函数每次向前移动一个样点数，所以可以得到N个分段信号。

步骤S102中，对所述分段信号进行傅里叶级数展开并引入时变傅里叶系数，以将所述分段信号表示为包含时变傅里叶系数的形式。

在本发明实施例中，将分段信号x_l(m)用傅里叶级数展开式表示为：

其中，t_l,m＝(m+l-h)T为信号x_l(m)对应的时刻，T为信号的时间采样间隔，a_k和b_k为傅里叶系数，F为频率采样间隔。在离散傅里叶变换中，频率是有限的，令最大频率f_max＝KF，则k的取值范围变为[0,K]。令

为傅里叶系数(其中b₀＝0)，

代表傅里叶基函数，上式可简写为：

上述傅里叶系数

只能反映全局信号的频率特征，不能表征每个时刻及其邻域对应的频率特征，而非平稳信号的频率特征随着时间发生变化，因此假设系数

随时间变化，此时公式(4)可改写为：

其中，

为时变傅里叶系数，可以表征非平稳信号的时频特征。

步骤S103中，根据所述包含时变傅里叶系数的形式，利用流式算法自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱。

在本发明实施例中，如图2所示，所述步骤S103，包括：

步骤S201中，对所述包含时变傅里叶系数的形式进行矩阵变换处理，求解所述时变傅里叶系数的最小二乘解。

在本发明实施例中，根据公式(5)可以得到第l个分段信号的第m个样点数据x_l(m)的矩阵表达式：

可以看出式(6)为数学欠定问题，其未知数个数大于方程个数，方程的解存在不唯一性。

为了解决这一数学欠定问题，本发明基于流式计算框架，使时变系数随着数据点的变化而连续更新，具体分为下面几个步骤：

假设相邻两个数据点的系数相似，此时数据x_l(m)可用超定方程组表示：

上式可以用简化的块矩阵表示为：

其中，

为数据x_l(m)对应的基函数，

代表数据x_l(m)对应的时变傅里叶系数，

为前一个样点的数据x_l(m-1)对应的时变傅里叶系数，ε是控制相邻两个采样点的数据对应的傅里叶系数相似程度的常数因子，I为单位矩阵。此时系数C_l(m)可由下面最小二乘问题求得：

上式的最小二乘解为：

C_l(m)＝(Ψ_l(m)^TΨ_l(m)+ε²I)^-1(x_l(m)Ψ_l(m)^T+ε²C_l(m-1)) (10)

步骤S202中，根据所述时变傅里叶系数的最小二乘解以及谢尔曼-莫里森公式，确定流式计算框架下的自适应时变傅里叶系数。

在本发明实施例中，利用线性代数中的谢尔曼-莫里森(Sherman-Morrison)公式直接将公式(10)中的逆矩阵变换为如下形式：

其中，Ψ_l(m)为行向量，Ψ_l(m)^T代表Ψ_l(m)的转置(列向量)，所以分子中的Ψ_l(m)Ψ_l(m)^T为常量，分母中的Ψ_l(m)^TΨ_l(m)是一个矩阵。这种逆矩阵的变换形式避免了矩阵逆的直接求解或迭代运算，提高了计算速度。

将公式(11)带入公式(10)并进行简化可以得到流式计算框架下的自适应时变傅里叶系数：

从上式(12)可以看出，利用流式理论计算的系数通过加上一定尺度缩放的基函数进行自适应更新，且更新系数只需要初等代数运算(向量点积)，不需要进行迭代。

步骤S203中，按照所述流式计算框架自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱。

在本发明实施例中，设置每个分段信号的初始时变傅里叶系数为零值，重复步骤S201和S202使系数C_l(m)按照上述流式计算框架自适应更新：

初始值→C_l(0)→C_l(1)→…→C_l(M)，m＝0,1,…,M (13)

将该分段信号最终更新的系数

定义为第l个分段信号(即第l个时间样点)对应的傅里叶系数C(l)＝[C₀(l),C₁(l),…,C_K(l)]^T。其中，C_k(l)＝[a_k(l),b_k(l)]^T，则第l个分段信号(即第l个时间样点)的时频谱可以定义为：

其中，l和k分别代表时间和频率。基于流式算法的自适应时频变换可以通过Nyquist频率或者由使用者规定该时频变换的频率范围和频率采样间隔。

令l＝0,1,…,N-1，重复步骤S101、S102和S103，得到每个分段信号的傅里叶系数。将每个分段信号的傅里叶系数组合在一起，可以得到整体信号的时频谱。

根据下式(15)计算该时频变换的反变换：

本发明实施例提供的基于流式算法的自适应时频变换方法，通过对实数离散信号进行滑动时窗处理，利用傅里叶级数匹配目标信号并引入时变傅里叶系数，将傅里叶分析与非平稳回归分析相结合，通过流式计算理论使时变傅里叶系数随着数据的改变而自适应更新，可以有效的表征非平稳信号的时频特征。本发明解析计算正则化条件下的数学欠定问题，利用谢尔曼-莫里森公式直接将最小二乘解中的逆矩阵转化为初等代数运算，避免了局部时频变换中的复杂迭代算法，有效的节约了计算过程中的内存成本；同时，相对于短时傅里叶变换，本发明参数选取灵活，可以自由调节频率采样间隔和频率范围，在保证时频分辨率的情况下去掉不必要的频率信息，极大的节省了数据存储空间。

仿真实例1

图3为一个简单的合成信号，该信号包含两个交叉的线性调频信号，时间采样间隔T为2ms，最大传播时间为1s。图4为该线性调频信号的理论频率。图5为经过基于流式算法的自适应时频变换处理后的时频谱图。在进行时频变换时，设置窗口长度M为80，参数ε为0.001，频率采样间隔和频率范围通过Nyquist频率来决定。从图5可以看出，本发明方法可以在时频谱中直观的将两个线性调频信号分解开来，准确的刻画频率随时间的变化关系，有效的反映信号的时频特征。图6为利用公式(15)反变换重构得到的信号，图7为反变换后的信号与原始合成信号的差值，可以看出本发明方法具有较高的重构精度。

仿真实例2

图8为某地区的实际地震数据，时间长度为4s，时间采样间隔T为2ms，时间样点总数N为2001。短时傅里叶变换通常基于分段快速傅里叶变换，最高采样频率通过Nyquist频率来决定，即

所以频率范围为[0,250]Hz，频率采样间隔

(N_fft＝2048为2的整数次幂)，频率总采样点数为1025，则短时傅里叶变换时频谱需要存储2001×2005个数据。当最高采样频率大于实际数据的峰值频率时，时频谱中高于数据峰值频率的范围不包含有用信息，这浪费了存储空间。本发明方法可以根据实际情况选取频率范围和频率采样间隔，图9为经过本发明时频变换方法处理后的时频谱图，选取的窗口长度M为60，参数ε为0.001，频率采样间隔为F＝0.244141Hz，频率范围为[0,150]Hz，频率总采样点数为615，则本发明时频谱需要存储2001×615个数据。可以看出，在不影响时频分辨率的前提下，本发明方法极大的节省了存储成本。

综上，本发明提供了一种基于流式算法的自适应时频变换，能够准确刻画信号的频率随时间的变化关系，表征非平稳信号的时频特征；一方面，本发明利用流式算法自适应的更新时变傅里叶系数，通过谢尔曼-莫里森公式直接将最小二乘解中的逆矩阵转化为初等代数运算，避免了迭代运算，降低了计算成本；另一方面，本发明利用傅里叶级数匹配目标信号，可以根据实际情况自主选取时频变换的频率范围和频率采样间隔，只保留有效的频率范围，极大的节省了数据存储空间。

以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于流式算法的自适应时频变换方法，其特征在于，包括：

对实数离散信号进行滑动时窗处理，得分段信号；

对所述分段信号进行傅里叶级数展开并引入时变傅里叶系数，以将所述分段信号表示为包含时变傅里叶系数的形式；

根据所述包含时变傅里叶系数的形式，利用流式算法自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱；

所述根据所述包含时变傅里叶系数的形式，利用流式算法自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱的步骤，包括：

对所述包含时变傅里叶系数的形式进行矩阵变换处理，求解所述时变傅里叶系数的最小二乘解；

根据所述时变傅里叶系数的最小二乘解以及谢尔曼-莫里森公式，确定流式计算框架下的自适应时变傅里叶系数；

按照所述流式计算框架自适应更新所述时变傅里叶系数，确定所述实数离散信号的局部时频谱；

所述流式计算框架下的自适应时变傅里叶系数为：

其中，C_l(m-1)为分段信号前一个采样点对应的时变傅里叶系数，ε为控制同一分段信号相邻两个采样点对应的傅里叶系数相似程度的常数因子。

2.根据权利要求1所述的基于流式算法的自适应时频变换方法，其特征在于，所述对实数离散信号进行滑动时窗处理，得分段信号的步骤中，所述分段信号的具体计算公式为：

x_l(m)＝x(m+l-h)w(m)

其中，x(n)(n＝0,1,…,N-1)为实数离散信号；m＝0,1,…,M-1为分段信号的采样点数；M为时窗的长度；若M为偶数，h＝M/2；若M为奇数，h＝(M-1)/2；x_l(m)(l＝0,1,…,N-1)代表第l个分段信号；w(m)代表窗函数。

3.根据权利要求2所述的基于流式算法的自适应时频变换方法，其特征在于，所述窗函数为汉宁窗，所述汉宁窗的表达式为：

。

4.根据权利要求1所述的基于流式算法的自适应时频变换方法，其特征在于，所述将所述分段信号表示为包含时变傅里叶系数和傅里叶基函数的形式的步骤中，具体表示形式为：

其中，

为时变傅里叶系数；

代表傅里叶基函数；t_l,m＝(m+l-h)T为分段信号x_l(m)对应的时刻；F为频率采样间隔。

Images