CN105572473B - 高分辨率线性时频分析方法 - Google Patents

Abstract

本发明公开了一种高分辨率线性时频分析方法，其采用迭代自适应方法获取高分辨率瞬时频率估计结果，相比线性时频分析方法，时间‑频率二维分辨率得到大幅提高，可显著提升多个信号在时间‑频率域内的二维区分度；且采用线性变换对信号进行分析处理，相比非线性时频分析方法，有效消除了交叉项的影响，可进一步提升对多个时频混叠信号的时频分析处理效果；本发明通过少量迭代即可获得收敛的时频分析结果，并对线性调频信号和非线性调频信号均具备良好的分析处理性能，相比其它改进类时频分析方法，计算量大幅降低且适用性更强。

Description

高分辨率线性时频分析方法

技术领域

本发明属于无线电信号处理技术领域，特别是一种高分辨率线性时频分析方法。

背景技术

现代战争，特别是20世纪90年代发生的几起局部战争表明：无线电电子侦察技术已经从早期的辅助保障作用发展成为涉及整个作战行动的重要影响因素。电子对抗是电子战(EW，Electronic Warfare)中的重要分支，电子战微波接收机是电子侦察系统以及情报保障体系中的重要组成部分，其主要功能是接收复杂信号环境下的非合作信号，侦察接收的性能直接影响干扰的效果。它的主要侦察对象为雷达信号和通信信号。现代电子侦察信号环境中，信号形式越来越复杂，信号数量越来越密集，所处频段也在不断扩展。因此，在同一时间内接收到多个时域、频域同时混叠且具有不同调制特性信号的概率较高。因此，电子战和电子侦察接收机除了要求具有宽的接收带宽、高灵敏度、大动态范围、高分辨率和实时或准实时信号处理能力以外，还必须能够适应现代高密集度的信号环境，具有对多个时域、频域混叠信号进行分析处理的能力。

传统的信号分析一般在时域或者频域展开，分别构成了信号的时域分析或频域分析方法，使用的主要数学工具是傅里叶变换(Fourier Transform，FT)，只适用于统计量不随时间变化的平稳信号。虽然Kalman滤波、RLS算法等自适应滤波方法也可用于非平稳信号的分析，但是一般只限于对慢时变信号的跟踪，并不能得到一般时变信号的统计量等结果(如信号功率谱等)。换言之，这些信号处理方法不能满足多个时域、频域混叠信号对应的非平稳信号分析的特殊要求。

对非平稳信号的时频分析可以分为线性分析和非线性分析两大类，分别对应两种经典的时频分析方法：短时傅里叶变换法(Short Time Fourier Transform，STFT)和魏格纳-威尔变换法(Wigner-Ville Transform，WVT)。由于WVT方法为基于二次型的分析方法，在分析多信号时会受到交叉项的影响，限制了其在多分量信号分析中的应用效果。通过对传统STFT时频分析方法的研究可知，在对包含多个分量的时域混叠信号进行分析时，为了获取足够的瞬时频率分辨率，需要尽可能增大STFT中的分析窗函数宽度；然而，为了在时频分析中获取信号的局部特性(时域分辨率)，分析窗函数宽度又不能超过非平稳信号的局部平稳长度。因此，STFT中的滑窗操作需要在频域分辨率和时域分辨率之间寻求平衡，这也成为制约STFT瞬时频率分辨率的根本原因。

除了上述两种经典时频分析方法以外，国内外学者还提出了多种改进方法，如L.B.Almeida将分数阶Fourier变换引入信号处理领域，利用其线性变换性质消除交叉项影响，但是该方法仅适用于对线性调频信号的分析，应用范围有限；Barbarossa.S等人在1995年的IEEE Trans.on SP上发表的文章《Analysis of multicomponent LFM signals by acombined Wigner-Hough transform》中，提出利用魏格纳-霍夫变换(Wigner-HoughTransform，WHT)对混合线性调频信号进行分析处理，但是由于WHT方法计算量较大，影响了其工程应用价值；齐林等人在2003年的《中国科学》上发表的文章《基于分数阶Fourier变换的多分量LFM信号的检测和参数估计》中，基于分数阶傅里叶变换提出一种多分量线性调频信号分析技术，并利用牛顿迭代法降低计算复杂度；戴琼海等人在1997年的《电子学报》上发表的文章《基于Randon-STFT变换的含噪LFM信号子空间分解》中，将STFT与Randon变换相结合，提出了基于Randon-STFT变换进行线性调频信号检测分析；但是上述两中方法均局限于对线性调频信号的分析，其对于复杂调制信号的处理性能有待进一步考察。

发明内容

本发明的目的在于提供一种能够突破非平稳信号局部平稳长度对STFT线性时频分析方法瞬时频率分辨率的限制，从而获取不含交叉项的高分辨时频分布结果的高分辨率线性时频分析方法。

实现本发明目的的技术解决方案为：

一种高分辨率线性时频分析方法，其特征在于：包括以下步骤：

步骤1：对待分析信号进行截取，获取子段观测数据；

步骤2：构造频域导向矢量和导向矩阵；

2a)构造如下频域导向矢量：

其中，e^j·表示复指数，f_i表示分析频点，i＝1,…,N；

2b)构造如下导向矩阵：

A(f)＝[a(f₁),a(f₂),…,a(f_N)]

其中，f＝[f₁,f₂,…,f_N]；

步骤3：采用迭代自适应方法对子段观测数据x_l进行高分辨率瞬时频率估计；

3a)初始化子段观测数据x_l的相关矩阵：

其中，[·]^H表示共轭转置操作；

3b)初始化功率矩阵：设功率矩阵为N阶对角阵则对角阵的对角元素为：

3c)据功率矩阵对相关矩阵进行更新，得到更新后的相关矩阵

3d)根据更新后的相关矩阵重新计算功率矩阵重新计算后，对角阵的对角元素为：

3e)重复进行3c)、3d)步操作直至功率矩阵对角线元素估计值的平均变化量小于其均值的5％，可认为此时的功率矩阵估计结果已趋于稳定，其对角线元素构成的数据矢量即为相应的高分辨率瞬时频谱估计结果，其中diag(·)表示取矩阵对角线元素构成数据矢量的操作；

步骤4：对滑动截取得到的L段观测数据向量x_l(l＝1,…,L)依次采用步骤3所述的迭代自适应方法进行高分辨率瞬时频率估计，得到相应的高分辨率瞬时频谱估计结果按顺序排列成矩阵，即为最终的高分辨率线性时频分布结果：

本发明与现有技术相比，其显著优点：

(1)本发明采用迭代自适应方法获取高分辨率瞬时频率估计结果，相比线性时频分析方法，时间-频率二维分辨率得到大幅提高，可显著提升多个信号在时间-频率域内的二维区分度。

(2)本发明采用线性变换对信号进行分析处理，相比非线性时频分析方法，有效消除了交叉项的影响，可进一步提升对多个时频混叠信号的时频分析处理效果；

(3)本发明通过少量迭代即可获得收敛的时频分析结果，并对线性调频信号和非线性调频信号(复杂调制信号)均具备良好的分析处理性能，相比其它改进类时频分析方法，计算量大幅降低且适用性更强。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明高分辨率线性时频分析方法的流程图。

图2是本发明采用不同方法仿真的时间-频率域二维分布结果图；图2(a)是多分量信号的WVD方法时频分析结果，图2(b)是多分量信号的STFT方法时频分析结果，图2(c)是采用本发明方法的多分量信号时频分析结果。

图3是本发明采用不同方法仿真的时频分布在某一时刻的瞬时频谱图；图3(a)是WVD方法时频分析结果对应的时间切片，图3(b)是STFT方法时频分析结果对应的时间切片，图3(c)是采用本发明方法所得时频分析结果对应的时间切片。

具体实施方式

结合图1：

本发明一种高分辨率线性时频分析方法，包括以下步骤：

步骤1：对待分析信号进行截取，获取子段观测数据；

步骤2：构造频域导向矢量和导向矩阵；

2a)构造如下频域导向矢量：

其中，e^j·表示复指数，f_i表示分析频点，i＝1,…,N；

2b)构造如下导向矩阵：

A(f)＝[a(f₁),a(f₂),…,a(f_N)]

其中，f＝[f₁,f₂,…,f_N]；

步骤3：采用迭代自适应方法对子段观测数据x_l进行高分辨率瞬时频率估计；

3a)初始化子段观测数据x_l的相关矩阵：

其中，[·]^H表示共轭转置操作；

3b)初始化功率矩阵：设功率矩阵为N阶对角阵则对角阵的对角元素为：

3c)据功率矩阵对相关矩阵进行更新，得到更新后的相关矩阵

3d)根据更新后的相关矩阵重新计算功率矩阵重新计算后，对角阵的对角元素为：

3e)重复进行3c)、3d)步操作直至功率矩阵对角线元素估计值的平均变化量小于其均值的5％，可认为此时的功率矩阵估计结果已趋于稳定，其对角线元素构成的数据矢量即为相应的高分辨率瞬时频谱估计结果，其中diag(·)表示取矩阵对角线元素构成数据矢量的操作；

步骤4：对滑动截取得到的L段观测数据向量x_l(l＝1,…,L)依次采用步骤3所述的迭代自适应方法进行高分辨率瞬时频率估计，得到相应的高分辨率瞬时频谱估计结果按顺序排列成矩阵，即为最终的高分辨率线性时频分布结果：

步骤1的具体过程如下：

1a)设Δω为期望达到的时频分析瞬时频率分辨率，由测不准原理可知，短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt满足下列关系式：

ΔtΔω≥0.5

1b)取分析窗函数宽度为短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt的10％，设待分析信号的采样周期为T_s，则本发明中分析窗函数宽度对应的采样点数为：

为保证时频分析中获取瞬时频谱时傅立叶变换的计算效率，将以采样点数表示的分析窗函数宽度取为：

其中，表示向上取整操作；

1c)设以采样点数表示的待分析信号长度为N，整个截取过程中共截取得到L段观测数据，令数据向量x_l(l＝1,…,L)为分析窗函数在l时刻截取的子段观测数据，则x_l可表示为

x_l＝[x₁,x₂,…,x_M]^T

其中[·]^T表示转置操作且x_k(k＝1,…,M)为子段观测数据x_l中的第k个采样点。

步骤1a)中Δt的选取范围如下：为保证时频分析结果的时间分辨率，短时傅里叶变换线性时频分析方法中的分析窗函数宽度应取最短，即取为：

步骤1c)中的截取方式采用滑动截取：整个滑动过程中共截取得到L段观测数据为：

L＝N-M。

步骤1c)中的截取方式采用分块截取：整个分块截取过程中共截取得到L段观测数据为：

其中，表示向下取整操作。

实施例1：

采用滑动截取的一种高分辨率线性时频分析方法实现步骤如下：

步骤1：对待分析信号进行滑动截取，获取子段观测数据。

1a)设Δω为期望达到的时频分析瞬时频率分辨率，则由测不准原理可知，短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt满足下列关系式：

ΔtΔω≥0.5

出自《现代信号处理》，清华大学出版社，2002，第362页

为了在时频分析结果中获取最高的时间分辨率，短时傅里叶变换线性时频分析方法中的分析窗函数宽度应该取最小值，即

按照上式计算时，时频分析过程的计算速度最低；为了提升时频分析过程中计算速度，可以在范围内增大分析窗函数宽度。

1b)取分析窗函数宽度为短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt的10％，设待分析信号的采样周期为T_s，则本发明中分析窗函数宽度对应的采样点数为：

为保证时频分析中获取瞬时频谱时傅立叶变换的计算效率，将以采样点数表示的分析窗函数宽度取为：

其中，表示向上取整操作；

1c)设待分析信号长度(采样点数)为N，截取方式采用滑动截取，整个滑动截取过程中共截取得到L段观测数据为：

L＝N-M

令数据向量x_l(l＝1,…,L)为本发明分析窗函数在l时刻截取的子段观测数据，则x_l可表示为：

x_l＝[x₁,x₂,…,x_M]^T

其中[·]^T表示转置操作且x_k(k＝1,…,M)为子段观测数据x_l中的第k个采样点；

步骤2：构造频域导向矢量和导向矩阵。

2a)构造如下频域导向矢量：

其中，i＝1,…,N；

2b)构造如下导向矩阵：

A(f)＝[a(f₁),a(f₂),…,a(f_N)]

其中，f＝[f₁,f₂,…,f_N]；

步骤3：采用迭代自适应方法对子段观测数据x_l进行高分辨率瞬时频率估计。

3a)初始化子段观测数据x_l的相关矩阵：

其中，[·]^H表示共轭转置操作；

3c)据功率矩阵对相关矩阵进行更新：

3d)根据更新后的相关矩阵重新计算功率矩阵：

3e)重复进行3c)、3d)步操作直至功率矩阵对角线元素估计值的平均变化量小于其均值的5％，可认为此时的功率矩阵估计结果已趋于稳定，其对角线元素构成的数据矢量即为相应的高分辨率瞬时频谱估计结果，其中diag(·)表示取矩阵对角线元素构成数据矢量的操作。

步骤4：对滑动截取得到的L段观测数据向量x_l(l＝1,…,L)依次采用步骤3所述的迭代自适应方法进行高分辨率瞬时频率估计，得到相应的高分辨率瞬时频谱估计结果按顺序排列成矩阵，即为最终的高分辨率线性时频分布结果：

本发明的效果可以通过以下仿真结果进一步说明。

1、仿真数据：

设仿真数据由两个不同频率的短脉冲信号、一个平方调频信号和一个正弦调频信号组合而成。其中，两个短脉冲信号在时频域内较为接近，而平方调频信号和正弦调频信号的时频支撑区存在交叉。

2、仿真内容及结果

分别对WVD时频分析方法、STFT时频分析方法和本发明方法的时频分析结果进行仿真分析，得到相应的时间-频率域二维分布结果以及某一时刻的时间切片(瞬时频谱)。其中：

图2为不同方法得到的时间-频率域二维分布结果。

图2(a)是多分量信号的WVD方法时频分析结果，图2(b)是多分量信号的STFT方法时频分析结果，图2(c)是采用本发明方法的多分量信号时频分析结果。

由图2(a)可见由于交叉项的存在，根据WVD的时频分析结果无法对有用信号的时频支撑区进行有效区分，由图2(b)可见STFT的能量聚集性较差，时频分析的二维分辨率较低，由图2(c)可见相应的时频分析结果中不存在交叉项，且具有较高的时频分辨率，各信号分量的时频支撑区清晰可分。

图3为不同时频分析结果在某一时刻的切片(瞬时频谱)。

图3(a)是WVD方法时频分析结果对应的时间切片，图3(b)是STFT方法时频分析结果对应的时间切片，图3(c)是采用本发明方法所得时频分析结果对应的时间切片。

由图3(a)可见WVD的瞬时频率分辨率较高，但是由于交叉项的存在，导致有用信号的频率无法识别，图3(c)为采用所提方法时，瞬时频谱曲线随迭代次数的收敛过程，其中第3次迭代时相应功率谱估计值的平均变化量已降至其均值的1.28％，可以发现：经过2～3迭代操作，即可获得较为稳定的瞬时频谱估计结果，与图3(b)可中STFT所得时间切片对比可知，所提方法的瞬时频率分辨率得到了大幅提高。

实施例2：

采用分块截取的一种高分辨率线性时频分析方法实现步骤如下：

步骤1.对待分析信号进行分块截取，获取子段观测数据。

1a)设Δω为期望达到的时频分析瞬时频率分辨率，则由测不准原理可知，短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt满足下列关系式：

ΔtΔω≥0.5

出自《现代信号处理》，清华大学出版社，2002，第362页，

为了在时频分析结果中获取最高的时间分辨率，短时傅里叶变换线性时频分析方法中的分析窗函数宽度应该取最小值，即：

按照上式计算时，时频分析过程的计算速度最低；为了提升时频分析过程中计算速度，可以在范围内增大分析窗函数宽度。

1b)取分析窗函数宽度为短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt的10％，设待分析信号的采样周期为T_s，则本发明中分析窗函数宽度对应的采样点数为：

为保证时频分析中获取瞬时频谱时傅立叶变换的计算效率，将以采样点数表示的分析窗函数宽度取为：

其中，表示向上取整操作；

1c)设待分析信号长度(采样点数)为N，截取方式采用分块截取，整个分块截取过程中共截取得到L段观测数据为：

其中，表示向下取整操作；

令数据向量x_l(l＝1,…,L)为本发明分析窗函数在l时刻截取的子段观测数据，则x_l可表示为：

x_l＝[x₁,x₂,…,x_M]^T

其中[·]^T表示转置操作且x_k(k＝1,…,M)为子段观测数据x_l中的第k个采样点；

步骤2.构造频域导向矢量和导向矩阵。

2a)构造如下频域导向矢量：

其中，i＝1,…,N；

2b)构造如下导向矩阵：

A(f)＝[a(f₁),a(f₂),…,a(f_N)]

其中，f＝[f₁,f₂,…,f_N]；

步骤3.采用迭代自适应方法对子段观测数据x_l进行高分辨率瞬时频率估计。

3a)初始化子段观测数据x_l的相关矩阵：

其中，[·]^H表示共轭转置操作；

3b)初始化功率矩阵：设功率矩阵近似为N阶对角阵则其对角元素为：

3c)据功率矩阵对相关矩阵进行更新：

3d)根据更新后的相关矩阵重新计算功率矩阵：

3e)重复进行3c)、3d)步操作直至功率矩阵对角线元素估计值的平均变化量小于其均值的5％，可认为此时的功率矩阵估计结果已趋于稳定，其对角线元素构成的数据矢量即为相应的高分辨率瞬时频谱估计结果，其中diag(·)表示取矩阵对角线元素构成数据矢量的操作。

步骤4.对滑动截取得到的L段观测数据向量x_l(l＝1,…,L)依次采用步骤3所述的迭代自适应方法进行高分辨率瞬时频率估计，得到相应的高分辨率瞬时频谱估计结果按顺序排列成矩阵，即为最终的高分辨率线性时频分布结果：

实施例1与实施例2区别在于步骤1c)中所述的获取子段观测数据的截取方法不同：采用实施例1中所述的滑动截取方法，其优势是相比实施例2中所述的分块截取方法，可以获得更高的时频分析精度；采用实施例2中所述的分块截取方法，其优势是相比实施例1中所述的滑动截取方法，可以获得更快的时频分析速度。

Claims (4)

1.一种高分辨率线性时频分析方法，其特征在于：包括以下步骤：

步骤1：对待分析信号进行截取，获取子段观测数据，具体过程如下：

1a)设Δω为期望达到的时频分析瞬时频率分辨率，由测不准原理可知，短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt满足下列关系式：

ΔtΔω≥0.5

1b)取分析窗函数宽度为短时傅里叶变换线性时频分析方法中的分析窗函数宽度Δt的10％，设待分析信号的采样周期为T_s，则本发明中分析窗函数宽度对应的采样点数为：

为保证时频分析中获取瞬时频谱时傅立叶变换的计算效率，将以采样点数表示的分析窗函数宽度取为：

其中，表示向上取整操作；

1c)设以采样点数表示的待分析信号长度为N，整个截取过程中共截取得到L段观测数据，令数据向量x_l(l＝1，…，L)为分析窗函数在l时刻截取的子段观测数据，则x_l可表示为

x_l＝[x₁，x₂，…，x_M]^T

其中[·]^T表示转置操作且x_k(k＝1，…，M)为子段观测数据x_l中的第k个采样点；

步骤2：构造频域导向矢量和导向矩阵；

2a)构造如下频域导向矢量：

其中，e^j.表示复指数，f_i表示分析频点，i＝1，…，N；

2b)构造如下导向矩阵：

A(f)＝[a(f₁)，a(f₂)，…，a(f_N)]

其中，f＝[f₁，f₂，…，f_N]；

步骤3：采用迭代自适应方法对子段观测数据x_l进行高分辨率瞬时频率估计；

3a)初始化子段观测数据x_l的相关矩阵：

其中，[·]^H表示共轭转置操作；

3b)初始化功率矩阵：设功率矩阵为N阶对角阵则对角阵的对角元素为：

3c)据功率矩阵对相关矩阵进行更新，得到更新后的相关矩阵

3d)根据更新后的相关矩阵重新计算功率矩阵重新计算后，对角阵的对角元素为：

其中，

3e)重复进行3c)、3d)步操作直至功率矩阵对角线元素估计值的平均变化量小于其均值的5％，可认为此时的功率矩阵估计结果已趋于稳定，其对角线元素构成的数据矢量即为相应的高分辨率瞬时频谱估计结果，其中diag(·)表示取矩阵对角线元素构成数据矢量的操作；

步骤4：对滑动截取得到的L段观测数据向量x_l(l＝1，…，L)依次采用步骤3所述的迭代自适应方法进行高分辨率瞬时频率估计，得到相应的高分辨率瞬时频谱估计结果按顺序排列成矩阵，即为最终的高分辨率线性时频分布结果：

2.根据权利要求1所述的高分辨率线性时频分析方法，其特征在于：步骤1a)中Δt的选取范围如下：为保证时频分析结果的时间分辨率，短时傅里叶变换线性时频分析方法中的分析窗函数宽度应取最短，即取为：

3.根据权利要求1所述的高分辨率线性时频分析方法，其特征在于：步骤1c)中的截取方式采用滑动截取：整个滑动过程中共截取得到L段观测数据为：

L＝N-M。

4.根据权利要求1所述的高分辨率线性时频分析方法，其特征在于：步骤1c)中的截取方式采用分块截取：整个分块截取过程中共截取得到L段观测数据为：

其中，表示向下取整操作。