CN110514921B - 一种电力电子变换器非平稳信号中非线性现象的识别方法 - Google Patents

Abstract

本发明提供一种电力电子变换器非平稳信号中非线性现象的识别方法。该方法包括：对多频分量的非平稳信号进行EMD分解和基于互信息量的延拓，得到延拓时间序列。对延拓时间序列进行EMD分解，通过有效因子将无用的伪IMF分量滤除。将得到的有效IMF分量重构后得到重构信号。用复Morlet小波对重构信号进行连续小波变换，并截取原信号相应的时域部分，得到小波变换下的时‑频平面图。对重构信号进行VMD分解，将分解的信号进行希尔伯特变换，并截取原信号相应的时域部分，得到希尔伯特时‑频平面图。本发明所使用的方法能够观察到电力电子变换器中存在的非线性现象，并同时获取该现象发生大致的时间范围。

Description

一种电力电子变换器非平稳信号中非线性现象的识别方法

技术领域

本发明属于电力电子变换器稳定性分析领域，涉及一种电力电子变换器非平稳信号中非线性现象的识别方法。

背景技术

在电力电子变换器的稳定性分析中，频域分析方法经常作为一种辅助分析方法发挥着一定的作用。频域分析方法是常用的谱分析方法，如功率谱估计，傅里叶变换等等，其中最为主要的是傅里叶变换。FFT是对一个离散的时间序列进行频域的变换，这样做的好处是对于电力电子中一个混杂的信号，存在着各种频率的信号，如噪声、振荡、锁频、倍周期、准周期等等，通过时域无法观测到信号的特征，但是在频域范围内能够读取信号的频率特性。

在平稳信号的分析中，由于具有e^jwt这一算子，傅里叶变换有着其他谱分析方法所没有的快速性和普适性。但因为傅里叶谱是信号的整体统计特性，它是信号整个时域内的积分，没有局部化分析信号的能力，故多只适用于平稳信号的分析。

根据设计需求，对电力电子变换器进行控制的时候，往往由于参数调整不当而导致电力电子变换器出现非线性行为，如倍周期分岔、边界碰撞分岔、Hopf分岔、混沌现象等等。因此，为了准确地识别这些现象并及时地将其控制到正常状态，需要对上述的现象出现的时域和频域两方面的信息进行同时提取。

考虑到时域和频域的局部化矛盾的情形，1946年Gabor提出了STFT变换，1984年J.Morlet提出了小波变换的概念，时频分析这一分析方法越来越广泛地被人们引用。基于Heisenberg测不准原理，时域信号和其傅里叶变换的能量跨度是不可能同时趋于无穷小的，所以必定会丢失一些信息，只能在高频处缩小时间尺度以换来频域的高分辨率或者在低频区域增大时间尺度获取时域的高分辨率。小波变换是一种多尺度的时频分析方法，对不同频率用不同的尺度进行分析，但是存在两个应用上的问题，一个是小波函数的需要根据具体的应用场合进行选择，一个是小波变换是沿着时间方向平移，使得调制信号会随着包络方向传播，导致相位信息仅仅是局部的，失去了物理意义。为此，Pinnegar将S变换引入小波变换，在连续小波函数中乘上一个相位项，相当于STFT和小波变换的结合。这样利用S变换对相位进行了修正。20世纪60年代burg在分析地震信号时提出最大熵谱值估值理论，该理论克服了传统谱分析分辨率不高和频率漏泄等缺点，使得谱分析仅以一个新阶段，称之为现代谱分析阶段。

为了表征某一时间点的频率成分，需要将整体谱推广到局部谱中去。20世纪初Hilbert变换已经存在，但对于多频率的信号，计算获得的瞬时幅值和频率没有实际意义，因此没有得到广泛的应用；随着后来20世纪末的EMD方法的提出，1998年N.E.Huang等人将两者结合的方法Hilbert-Huang Transform(HHT)广泛地运用于信号的时频分析，该方法成功提取出了具有实际意义的频率分量；然而EMD分解会在低频区产生不合理的频率特征，即实际不存在的频率；且首先分解出的IMF分量具有的频率范围过宽，一定程度上降低了频率分辨率；此外，对于低能量的频率成分的分离效果较差。为了改进这一不足，Peng等利用小波包进行经验模态分解，保证了每个IMF处于窄频带范围。

本发明在电力电子变换器中引入了基于信号延拓、有效因子改进的时频分析法，通过使用其对电力电子变换器的非平稳信号进行分析，可以得到一定的识别效果。

发明内容

本发明的目的在于提供一种改进的时频分析方法观察整体时域内变换器的非线性现象的变化，涉及一种电力电子变换器系统参数阶跃变化时输出非平稳信号中非线性现象识别方法。克服了传统频域分析方法对含有直流分量的交流信号识别不精确的情况，并且采用信号延拓方法和筛选EMD分解伪分量的方法提高了时频分析的精确度，在得到多周期运行信息的同时获取倍周期分岔以及混沌发生的时间。

本发明通过如下技术方案实现。

步骤一、电力电子变换器中某一参数发生阶跃变化时，在相应时域内采集系统的状态变量的时间序列信号x(t)，信号长度为N；

步骤二、基于互信息的概念对信号x(t)的左端和右端分别进行信号延拓，，得到延拓后的信号X(t),相应的左延拓长度为N1，右延拓长度为N2，因此对应时间序列从[t_min,t_max]变为[t_min-N1×Δt,t_max+N2×Δt]，其中Δt为信号x(t)的采样间隔；对延拓时间序列X(t)进行集合经验模态分解(EMD,Empirical Mode Decomposition)分解，通过Pearson相关系数法将伪IMF(Intrinsic mode function，本征模态函数)分量滤除，得到IMF分量R-IMF，将R-IMF分量重构后，得到重构信号Y(t)；

步骤三、选择复Morlet小波(Complex Morlet,Cmor)对Y(t)进行连续小波变换，对求得的小波系数截取原时域部分，即截取[t_min,t_max]区间中对应的小波系数(在步骤二中进行了基于原信号x(t)与互信息量的延拓，为的是消除时频分析中存在的边界效应，因此需要在求得的小波系数中去掉该增加的部分)，得到小波变换下的时-频平面图；

步骤四、对重构信号Y(t)进行VMD分解(详细见权利要求6)，将分解的信号{u(t)}进行Hilbert变换，对求得的幅值矩阵(频率尺度×时间尺度的二维矩阵)原时域部分，即截取[t_min,t_max]区间(原因同步骤三)，得到Hilbert时-频平面图。

进一步的，在步骤二中，基于互信息量对电力电子变换器状态变量的时序序列进行延拓的过程包括:为消除或减弱EMD和Hilbert换变换普遍存在的较为严重的边界效应，考虑到端点延拓和镜像延拓对于非平稳信号效果不明显，这里引入一种基于互信息量或者条件熵的方法进行改进，对于信号X，Y，定义Y相对于X的信息量I为：

I(X,Y)＝H(X)-H(X|Y) (12)

H(X)表示信号X的信息熵，H(X|Y)表示Y与X的条件熵，条件熵值越大表示X与Y的相似性越小，差异度越大，信息量I用来表示两个信号的相似程度。

进一步的，信号延拓还包括以下过程：

1)计算待分析信号x(t)的极大值点，取横坐标为时间，纵坐标为信号幅值.

2)假设x(t)存在m个极大值点，对于左延拓而言，依据第一个极大值点的横坐标值与时间起始值确定数据左端截取长度N1，右端采用设置的默认长度值N2(与总极大值点数和总数据长度有关)；

对于右延拓而言，依据最后一个极大值点的横坐标值与时间结束值确定右端截取长度N2，左端采用设置的默认长度值N1(与总极大值点数和总数据长度有关)；

对于两种延拓情形，求取延拓数据段的步骤是类似的，区别在于N1和N2的选取。按照各极大值点所在数据段的相对位置不变一次进行数据左端N1个数据和右端N2个数据的截取，得到m组长度为N1+N2+1的数据段X₁,X₂,X₃,……,X_m；

3)式(2)分别计算了X₂,X₃,……,X_m与X₁的互信息量I₂,I₃,……,I_m，取互信息量最大值max(I₂,I₃,……,I_m)对应的下标最小(若最大值不唯一考虑最小下标，对应满足要求的数据段中离延拓的端点最近的一个)的互信息量I_i及其所对应的X_i作为参考数据段(具体作用见如下第4)点)。

式中，p(x,y)是当前X_i和X₁的联合概率密度函数，而p(x)和p(y)分别是X₁和X_i的边缘概率密度函数。

4)延拓数据段：当X₂与X₁的相似程度最大，则有两种获取延拓数据的方法：

i)读取原信号x(t)左端点的数值，记作A；在X₁右端点与X₂的极值点之间寻找与A点纵坐标相等的点，若能匹配到A点，则记为B₁,B₂,B₃,...,选择其中离X₂数据段极值点最近的一点作为延拓数据段的右端点，X₁的右端点作为延拓数据段的左端点，至此左延拓数据段的完成；

类似地，对于右延拓数据段的求取，首先读取原信号x(t)右端点的数值，记作A；在X₁左端点与X₂的极值点之间寻找与A点纵坐标相等的点，若能匹配到A点，则记为B₁,B₂,B₃,...,选择其中离X₂数据段极值点最近的一点作为延拓数据段的左端点，X₁的左端点作为延拓数据段的右端点，至此右延拓数据段的完成；

ii)若i)方法未能实施，即未能找到与匹配A点数据值匹配的点，计算X₁数据段的左端点与X₂数据段的左端点的差值，并将该差值叠加到X₁的右端点到X₂的左端点之间的数据中，至此左延拓完成。

类似地，对于右延拓数据段的求取，计算X₁数据段的右端端点与X₂数据段的右端点的差值，并将该差值叠加到X₁的左端点到X₂的右端点之间的数据中，至此右延拓完成。

进一步的，在步骤二中，重构信号的Y(t)获取具体过程如下：

对EMD分解得到的m个IMF分量，通过公式(4)所示的Pearson相关系数法进行筛选，得到有效因子Q_j：

式中，IMF(j,i)表示第j个IMF分量的第i个离散值Q_j表示第j个IMF分量与原信号的相关程度，包含的原信号频率成分越多该系数越大；x(i)对应长度为N的离散时间序列{x(i)|i＝1,...,N}的第i个值。

对Q_j设定阈值筛选出l个R-IMF分量c₁,c₂,…,c_l。于是重构信号

Y(t)＝c₁+c₂+…+c_l (4)

进一步的，在步骤三中，为分析包含电力电子变换器非线性现象的状态变量的时间序列，基于Cmor小波对重构的信号进行连续小波变换Continuous Wavelet Transform(CWT)的过程为：

小波分析的困难往往在于选择出合适的小波函数。Cmor是一种非正交小波，是复小波函数的一种，具有对称性，不存在尺度函数，但该小波函数在时域和频域都具有很好的局域性。在电力信号的分析中这类小波往往具有较好地表现特性，在之后的分析中，本方案采用带宽参数f_b为5，中心频率f_c为3.5，支撑长度L为(-8,8),网格数N为1000(默认2⁸)的Cmor小波函数进行应用。

该小波函数是基于时域上进行定义的，表达式为：

式中，x为函数的自变量，这里表示时间，i为虚数单位。

观察上式，可得Cmor小波实质上是一种复余弦调制的Gaussian波，其傅里叶变换为：

式中，f表示频率。

若将电感电流信号作为待分析的状态变量，将其在该小波基函数下展开，即连续小波变换：

式中，x_L(t)表示采集的电感电流信号，Ψ_a,τ(t)表示小波基函数,满足

Ψ^*表示是单位化后小波基函数Ψ的共轭函数。表示单位化后的小波基函数，a表示尺度大小,τ表示位移大小，i表示虚数单位。

根据得到的连续小波系数可以绘制出小波变换下的时频图。

进一步的，对于步骤四中的VMD分解过程以及Hilbert时频图的求解过程为：

1)VMD分解过程

VMD将本征模态函数(IMF)定义为一个调幅-调频信号，即u_k(t)＝A_k(t)cos[φ_k(t)](A_k(t)为幅值函数，φ_k(t)为相位函数)。易得瞬时频率ω_k(t)＝d[φ_k(t)]/dt,若多分量信号x(t)由K(需要预先判定)个有限带宽的IMF分量u_k组成，且各IMF的中心频率为ω_k，VMD方法建立的约束变分模型为：

式中，{u_k}表示原信号分解得到的K个IMF分量的集合，{ω_k}表示各分量的中心频率的集合(可经过平均化初值之后迭代)，δ(t)为脉冲函数。

为解决式(8)所示的约束变分问题，将通过ALM将其转变为非约束性问题，即引入增广拉格朗日函数：

式中，α为二次惩罚因子，与普通的拉格朗日函数不同的是，为了其能够更好运用交替方向乘子法(ADMM,Alternating Direction Method of Multipliers)求解方程(9),添加了一个惩罚项，这使得问题更加偏向了凸优化，便于采用对偶上升算法进行全局寻优。

最后得到一个K×(原数据长度+延拓数据段长度)的二维IMF分量幅值矩阵。

3)Hilbert变换

对1)中得到的幅值矩阵进行Hilbert变换

式中，H_j(t)表示第j个IMF分量u_j序列对应的希尔伯特变换，j＝1,2,...,K；τ表示积分变量。

将H_j(t)中的时间t离散化后可以得到{H_j(i)|i＝1,...,N_s}，N_s表示延拓后的数据长度，将频率区间划分为N_f个部分，将{H_j(i)|j＝1,2,...,K}分别对应到相应地频率部分中去，最后可得到大小为N_f×N_s的二维时频矩阵

本发明相对于现有技术具有如下的优点及效果：

电力电子变换器中常用的分析非线性现象的方法往往需要建模，过程繁琐，而且每次建模对应着单一的电路拓扑结构，单一的频率分析方法如快速傅里叶变换含有的信息不够充分，且只适用于平稳信号的分析。本文提出了优化后的时频分析方法，分别通过小波分析方法和Hilbert时频分析方法来观察信号的时频图。通过时频图，可以直观的看出信号稳定性出现变化的时间点，此外通过采用基于互信息量的边界延拓减弱了边界效应的影响，提高了时频分析的准确性。

附图说明

图1a为本发明实施方式中的实施步骤流程图。

图1b为电流型buck电路原理图。

图2a为本发明实施方式中的原信号参考数据段图。

图2b为本发明实施方式中的左延拓信号图。

图3a为本发明实施方式中使用的小波函数波形实部图。

图3b为本发明实施方式中使用的小波函数波形虚部图。

图4为对应于图2中伴随参考电流I_ref阶跃跳动电感电流的时域波形。

图5a为本发明实施方式中buck电路模型中电感电流信号中参考数据段图。

图5b为本发明实施方式中buck电路模型中电感电流信号双边延拓信号图。

图6为本发明实施方式中buck电路模型中电感电流EMD分解重构信号图。

图7a为本发明实施方式中t∈(0.44,0.52)时跟随参考电流I_ref阶跃跳动电感电流的小波变换二维时-频图。

图7b为本发明实施方式中t∈(0.496,0.55)时跟随参考电流I_ref阶跃跳动电感电流的小波变换二维时-频图。

图8a为本发明实施方式中t∈(0.44,0.51)时跟随参考电流I_ref阶跃跳动电感电流的Hilbert变换二维时-频图。

图8b为本发明实施方式中t∈(0.51,0.55)时跟随参考电流I_ref阶跃跳动电感电流的Hilbert变换二维时-频图。

具体实施方式

以下结合附图对本发明的具体实施作进一步描述，但本发明的实施和保护不限于此。

如图1a所示，本实施方式对本发明方法进行详细描述：

步骤一、本实例使用电流型buck电路(输入电压E＝20V,电感L＝3.3mH,电容C＝1000uF,负载R＝19Ω，参考电压V_ref＝11.3V,频率f＝2.5kHz,VD表示二极管，S表示开关管)，如图1b所示。采集参考电流I_ref阶跃下电感电流的时域信号x_L(t),取信号长度N为9000。根据已有论文的分析，该电力电子变换器可随着I_ref的递增通过倍周期分岔走向混沌，设置I_ref阶跃跳动时刻，I_ref的阶跃变化波形与相应的电感电流时域波形x_L如图4所示的非平稳信号时域波形(在t＝0.45s的时候I_ref从0.6A跳跃到0.9006A，在t＝0.49s的时候I_ref从0.9006A跳跃到1.25A，在t＝0.53s的时候I_ref从1.25A跳跃到1.28A)。

步骤二、基于互信息的概念对信号x_L(t)的左端和右端分别进行信号延拓，左延拓情形见图5a，得到延拓后的信号X_L(t)如图5b(其中左延拓数据段的N1＝17,N2＝20，右延拓数据段的N1＝39,N2＝39，计算得到左延拓数据段总长度L1＝33，右延拓长度为L2＝23)。

对延拓时间序列X_L(t)进行EMD分解，通过Pearson相关系数法将无用的伪IMF(Intrinsic mode function)分量滤除，求得前五个IMF分量的相关系数满足要求，分别为0.9433、0.3468、0.4689、0.2544、0.1261，因此得到五个有效的IMF分量(R-IMF)。将R-IMF分量重构后，得到重构信号Y_L(t)，如图6。

步骤三、选择复Morlet(Complex Morlet,Cmor)小波，运用公式(7)对Y_L(t)进行连续小波变换，并截取原信号相应的时域部分，得到小波变换下的时-频平面图，如图7a所示的t∈(0.44,0.52)时的Hilbert变换时频图和图7b的t∈(0.496,0.55)Hilbert变换时频图。

步骤四、对重构信号Y(t)进行变分模态分解分解，得到频率由低向高排序的分量{u_k(t)|k＝1,2,...,K}，对{u_k(t)|k＝1,2,...,K}进行希尔伯特变换，同样截取原信号对应时域部分，相应的时频图，如图8a所示t∈(0.44,0.51)时周期一，周期二，周期四的Hilbert变换时频图和图8bt∈(0.51,0.55)时周期四和混沌状态的Hilbert变换时频图。

Claims (3)

1.一种电力电子变换器非平稳信号中非线性现象的识别方法，其特征在于:所述方法包括以下步骤:

步骤一、电力电子变换器中某一参数发生阶跃变化时，在相应时域内采集电力电子变换器的状态变量的时间序列信号x(t)，信号长度为N；

步骤二、基于互信息的概念对信号x(t)的左端和右端分别进行信号延拓，得到延拓后的信号X(t),相应的左延拓长度为N1，右延拓长度为N2，因此对应时间序列从[t_min,t_max]变为[t_min-N1×Δt,t_max+N2×Δt]，其中Δt为信号x(t)的采样间隔；对延拓时间序列X(t)进行集合经验模态分解(EMD,Empirical Mode Decomposition)，通过Pearson相关系数法将伪IMF(Intrinsic mode function，本征模态函数)分量滤除，得到IMF分量R-IMF，将R-IMF分量重构后，得到重构信号Y(t)；基于互信息量对电力电子变换器状态变量的时序序列进行延拓的过程包括:

基于互信息量或者条件熵的方法进行改进，对于信号X，Y，定义Y相对于X的信息量I为：

I(X,Y)＝H(X)-H(X|Y) (1)

H(X)表示信号X的信息熵，H(X|Y)表示Y与X的条件熵，条件熵值越大表示X与Y的相似性越小，差异度越大，信息量I用来表示两个信号的相似程度；

信号延拓还包括以下过程：

1)计算待分析信号x(t)的极大值点，取横坐标为时间，纵坐标为信号幅值；

2)假设x(t)存在m个极大值点，对于左延拓而言，依据第一个极大值点的横坐标值与时间起始值确定数据左端截取长度N1，右端采用设置的默认长度值N2；

对于右延拓而言，依据最后一个极大值点的横坐标值与时间结束值确定右端截取长度N2，左端采用设置的默认长度值N1；

对于两种延拓情形，求取延拓数据段的步骤是类似的，区别在于N1和N2的选取；按照各极大值点所在数据段的相对位置不变一次进行数据左端N1个数据和右端N2个数据的截取，得到m组长度为N1+N2+1的数据段X₁,X₂,X₃,……,X_m；

3)式(2)分别计算了X₂,X₃,……,X_m与X₁的互信息量I₂,I₃,……,I_m，取互信息量最大值max(I₂,I₃,……,I_m)对应的下标最小；

式中，p(x,y)是当前X_i和X₁的联合概率密度函数，而p(x)和p(y)分别是X₁和X_i的边缘概率密度函数；

4)延拓数据段：当X₂与X₁的相似程度最大，则有两种获取延拓数据的方法：

i)读取原信号x(t)左端点的数值，记作A；在X₁右端点与X₂的极值点之间寻找与A点纵坐标相等的点，若能匹配到A点，则记为B₁,B₂,B₃,...,选择其中离X₂数据段极值点最近的一点作为延拓数据段的右端点，X₁的右端点作为延拓数据段的左端点，至此左延拓数据段的完成；

类似地，对于右延拓数据段的求取，首先读取原信号x(t)右端点的数值，记作A；在X₁左端点与X₂的极值点之间寻找与A点纵坐标相等的点，若能匹配到A点，则记为B₁,B₂,B₃,...,选择其中离X₂数据段极值点最近的一点作为延拓数据段的左端点，X₁的左端点作为延拓数据段的右端点，至此右延拓数据段的完成；

ii)若i)方法未能实施，即未能找到与匹配A点数据值匹配的点，计算X₁数据段的左端点与X₂数据段的左端点的差值，并将该差值叠加到X₁的右端点到X₂的左端点之间的数据中，至此左延拓完成；

类似地，对于右延拓数据段的求取，计算X₁数据段的右端端点与X₂数据段的右端点的差值，并将该差值叠加到X₁的左端点到X₂的右端点之间的数据中，至此右延拓完成；

步骤三、选择复Morlet小波(Complex Morlet,Cmor)对Y(t)进行连续小波变换，对求得的小波系数截取原时域部分，即截取[t_min,t_max]区间中对应的小波系数，得到小波变换下的时-频平面图；

步骤四、对重构信号Y(t)进行VMD分解，将分解的信号{u(t)}进行Hilbert变换，对求得的幅值矩阵即频率尺度×时间尺度的二维矩阵原时域部分，即截取[t_min,t_max]区间时域部分，得到Hilbert时-频平面图，重构信号的Y(t)获取具体过程如下：

对EMD分解得到的m个IMF分量，通过公式(3)所示的Pearson相关系数法进行筛选，得到有效因子Q_j：

式中，IMF(j,i)表示第j个IMF分量的第i个离散值，Q_j表示第j个IMF分量与原信号的相关程度，包含的原信号频率成分越多该系数Q_j越大；x(i)对应长度为N的离散时间序列{x(i)|i＝1,...,N}的第i个值；

对Q_j设定阈值筛选出l个R-IMF分量c₁,c₂,...,c_l，得到重构信号

Y(t)＝c₁+c₂+...+c_l (4)。

2.根据权利要求1所述的一种电力电子变换器非平稳信号中非线性现象的识别方法，其特征在于:在步骤三中，基于Cmor小波对重构的信号进行连续小波变换的过程为：

采用带宽参数f_b为5，中心频率f_c为3.5，支撑长度为L，网格数N为1000的Cmor小波函数进行应用，Cmor小波函数是基于时域上进行定义的，表达式为：

式中，x为函数的自变量，这里表示时间，i为虚数单位；

公式(5)的傅里叶变换为：

式中，f表示频率；

根据得到的连续小波系数绘制出小波变换下的时频图。

3.根据权利要求1所述的一种电力电子变换器非平稳信号中非线性现象的识别方法，其特征在于:在步骤四中，基于VMD分解和Hilbert变换求着幅值矩阵的过程为：

1)VMD分解过程

VMD将本征模态函数(IMF)定义为一个调幅-调频信号，即u_k(t)＝A_k(t)cos[φ_k(t)]，其中A_k(t)为幅值函数，φ_k(t)为相位函数；易得瞬时频率ω_k(t)＝d[φ_k(t)]/dt,若多分量信号x(t)由K个有限带宽的IMF分量u_k组成，且各IMF的中心频率为ω_k，VMD方法建立的约束变分模型为：

式中，{u_k}表示原信号分解得到的K个IMF分量的集合，{ω_k}表示各分量的中心频率的集合，δ(t)为脉冲函数；

为解决式(8)所示的约束变分问题，通过ALM将其转变为非约束性问题，即引入增广拉格朗日函数；

与普通的拉格朗日函数不同的是，为了其能够更好运用交替方向乘子法(ADMM,Alternating Direction Method of Multipliers)求解,添加了一个惩罚项，这使得问题更加偏向了凸优化，便于采用对偶上升算法进行全局寻优；

最后得到一个K×(原数据长度+延拓数据段长度)的二维IMF分量幅值矩阵；

2)Hilbert变换

对1)中得到的幅值矩阵进行Hilbert变换

式中，H_j(t)表示第j个IMF分量c_j序列对应的希尔伯特变换，j＝1,2,...,K；τ表示积分变量；

将H_j(t)中的时间t离散化后可以得到{H_j(i)|i＝1,...,N_s}，N_s表示延拓后的数据长度，将频率区间划分为N_f个部分，将{H_j(i)|j＝1,2,...,K}分别对应到相应地频率部分中去，最后可得到大小为N_f×N_s的二维时频矩阵。

Images