一种基于参数预估与能量聚焦准则的自适应时频分析方法

本发明属于信号处理，具体涉及一种基于参数预估与能量聚焦准则的自适应时频分析方法。

背景技术：

1、由于非平稳信号的频率随时间变化，传统的一维时间或频率分析都无法有效地分析其特征。因此，为了能同时在时间维与频率维观察信号特征，时频分析被提出并应用于非平稳信号的分析。短时傅里叶变换是一种传统的时频分析方法，该算法由于算法复杂度低，在实际工程中得到了广泛的应用。但窗长的使用限制了其时间分辨率与频率分辨率，使得二者无法同时提高。为同时获取较高的时间分辨率与频率分辨率，学者提出了wigner-ville分布(wigner-ville distribution,wvd)。尽管该算法实现了高分辨率的需求，但在多分量和非线性调频信号的分析中，会将交叉项引入时频分布结果中，从而影响信号特征提取的准确性。为了抑制交叉项带来的影响，有许多算法被提出，主要可以分为如下两类：

2、(1)将信号变换到模糊域中，利用交叉项远离模糊域原点而自项靠近模糊域原点的性质，设计相应的模糊域滤波核进行滤波，再将滤波结果反变换到时频域中，代表算法有拓展修正b分布和高斯径向核分布等。这种思路的优势在于：计算量较小；线性调频信号的模糊域特征明显，易于设计核函数；多分量信号的交叉项远离模糊域原点，也可以通过设计低通的模糊域滤波器，滤除信号间的交叉项并保留信号自项。但其不足在于：非线性调频信号的模糊域特征聚焦性差，信号自项不易于分辨，难以设计理想的核函数。

3、(2)直接在时频域进行二维滤波，代表的算法有自适应方向时频分布(adaptivedirectional time-frequency distribution,adtfd)。这类算法的优势在于：通过设计合理的时频域滤波器，能够直接在信号的时频域上进行二维滤波，最大程度保留信号信息。但不足在于：二维时频滤波往往具有较大的计算量。

技术实现思路

1、发明目的：本发明目的在于通过参数预估与能量聚焦准则设计时频域滤波核的角度，在减少计算量的同时，对单分量信号与多分量信号、线性调频信号与双曲调频信号都能获得高分辨率的时频分布结果。

2、技术方案：为解决上述技术问题，本发明提出一种基于参数预估与能量聚焦准则的自适应时频分析方法，该方法包括以下步骤：

3、步骤1、获取待处理信号的采样数据序列s(a)，a为时间量离散索引；

4、步骤2、计算待处理信号采样数据序列s(a)的模糊函数az(m,n)与wigner-ville分布(简称wvd)ρwvd(a,b)，其中，m为频移量离散索引，n为时延量离散索引，a为时间量离散索引，b为频率量离散索引；

5、步骤3、利用radon变换对信号模糊函数az(m,n)进行检测，预估信号模糊域自项初始径向角集合ψ；

6、步骤4、利用自项初始径向角集合ψ设计时频域滤波核，计算时频域滤波结果及其能量聚焦程度；

7、步骤5、以步骤4计算得到的能量聚焦程度作为依据，通过能量聚焦准则对自项初始径向角集合ψ进行拓展，得到自项拓展后径向角集合φafter，并通过φafter计算信号的最优时频分布ρ(i,j)，其中，i为时频域时间量离散索引，j为时频域频率量离散索引。

8、进一步的，步骤1中，采用如下方法获取待处理信号的采样数据序列s(a)，具体包括以下步骤：从传感器接收l个采样点的采集数据作为待处理信号的采样数据序列s(a)，a＝0,1,…,l-1，其中，a为时间量离散索引，所述的l为包含整个信号持续时间所对应的采样点个数，l＝2r，r为大于等于10的整数。

9、进一步的，步骤2中，采用如下方法计算待处理信号采样数据序列s(a)的模糊函数az(m,n)和wigner-ville分布ρwvd(a,b)，具体包括以下步骤：

10、步骤2.1、计算待处理信号采样数据序列s(a)的离散傅里叶变换s(α)：

11、

12、其中，α为s(α)的频率量离散索引，j为虚数单位，即

13、步骤2.2、利用待处理信号的离散傅里叶变换s(α)生成s(a)的离散解析信号z(a)的傅里叶变换z(β)：

14、

15、其中，β为z(β)的频率量离散索引；

16、步骤2.3、对z(β)进行离散逆傅里叶变换得到待处理信号的离散解析信号z(a)：

17、

18、其中，a为时间量离散索引；

19、步骤2.4、计算离散解析信号z(a)的自相关函数kz(a,n)：

20、

21、其中，表示的共轭形式，n为时延量离散索引，n＝0,1,…,l-1，a＝0,1,…,l-1；

22、步骤2.5、依据离散解析信号z(a)的自相关函数kz(a,n)，计算待处理信号的模糊函数az(m,n)：

23、

24、其中，m为频移量离散索引，n为时延量离散索引；

25、步骤2.6、依据离散解析信号z(a)的自相关函数kz(a,n)，计算待处理信号的wigner-ville分布ρwvd(a,b)：

26、

27、其中，a为时间量离散索引，b为频率量离散索引。

28、进一步的，步骤3中，采用如下方法利用radon变换对信号模糊函数az(m,n)进行检测，预估信号模糊域自项初始径向角集合ψ，具体包括以下步骤：

29、步骤3.1、对预估信号模糊域自项初始径向角集合ψ的各参数进行初始化，具体包括以下参数的初始化：

30、①radon变换的离散径向角数量k初始化为：21≤k≤181的正整数；

31、②radon变换的离散径向角集合初始化为：

32、θ＝[θ1,θ2,…,θk,…,θk]

33、其中，k为径向角离散索引，k＝1,2,…,k，第k个离散径向角θk初始化为：

34、θk＝-π/2+kπ/k；

35、③信号模糊域自项初始径向角提取阈值rth初始化为：0.1≤rth≤0.5的正实数；

36、步骤3.2、计算待处理信号的模糊函数az(m,n)在径向距离为0处的radon变换，得到radon变换集合r(θ)＝[r(θ1),r(θ2),…,r(θk),…,r(θk)]，其中，第k个离散径向角θk所对应的径向距离为0处的radon变换结果r(θk)为：

37、

38、其中，k为radon变换的离散径向角数量，其值已在步骤3.1中初始化，δ()表示为：

39、

40、其中，x为实数，当x＝0时，δ(x)＝1，当x≠0时，δ(x)＝0；

41、步骤3.3、利用径向距离为0处的radon变换集合r(θ)估计信号模糊域自项初始径向角集合ψ：

42、ψ＝{θ:r(θ)≥rth}

43、其中，rth为信号模糊域自项初始径向角提取阈值，其值已在步骤3.1中初始化，r(θ)＝[r(θ1),r(θ2),…,r(θk),…,r(θk)]为步骤3.2计算得到的radon变换集合，为预估的模糊域自项初始径向角集合，其中，ψ表示radon变换集合r(θ)中大于等于提取阈值rth的radon变换值对应的角度集合，是预估的第p个模糊域自项初始径向角，p＝1,2,…,p，p为模糊域自项初始径向角集合ψ所包含的角度个数。

44、进一步的，步骤4中，采用如下方法利用信号自项初始径向角集合ψ设计时频域滤波核，计算时频域滤波结果及其能量聚焦程度，具体包括以下步骤：

45、步骤4.1、对设计时频域滤波核、计算时频域滤波结果及其能量聚焦程度的各参数进行初始化，具体包括以下参数的初始化：

46、①时频域滤波核的时间维长度参数vt初始化为：1≤vt≤5的正整数；

47、②时频域滤波核的频率维长度参数vf初始化为：1≤vf≤15的正整数；

48、③时频域滤波核长度lγ初始化为：ceil(l/16)≤lγ≤ceil(l/8)的正奇数，其中，ceil()为向上取整函数；

49、④时频域滤波核时间量t(aγ)初始化为：t(aγ)＝-1+aγ(2/lγ)，其中，aγ为滤波核时间量离散索引，aγ＝0,1,…,lγ-1；

50、⑤时频域滤波核频率量f(bγ)初始化为：f(bγ)＝-1+bγ(2/lγ)，其中，bγ为滤波核频率量离散索引，bγ＝0,1,…,lγ-1；

51、步骤4.2、令自项初始径向角的离散索引p＝1；

52、步骤4.3、利用预估的自项初始径向角按如下方式设计时频域滤波核γp(aγ,bγ)：

53、

54、其中，vt为时频域滤波核的时间维长度参数，其值已在步骤4.1中初始化，vf为时频域滤波核的频率维长度参数，其值已在步骤4.1中初始化，tp(aγ)和fp(bγ)分别为第p个自项初始径向角对应的时频域滤波核时间量和频率量，表示为：

55、

56、

57、其中，t(aγ)是时频域滤波核时间量，其值已在步骤4.1中初始化，aγ为滤波核时间量离散索引，aγ＝0,1,…,lγ-1，f(bγ)是时频域滤波核频率量，其值已在步骤4.1中初始化，bγ为滤波核频率量离散索引，bγ＝0,1,…,lγ-1，lγ为时频域滤波核长度，其值已在步骤4.1中初始化；

58、步骤4.4、利用时频域滤波核γp(aγ,bγ)和信号的wigner-ville分布ρwvd(a,b)，按如下方式计算时频域滤波结果ρp(i,j)：

59、

60、其中，i为时频域时间量离散索引，j为时频域频率量离散索引，表示对括号内所有可计算的数值矩阵中的每个元素求和，其中，a和b为满足以下条件的非负整数值：

61、

62、其中，l为包含整个信号持续时间所对应的采样点个数，其值已在步骤1中初始化，lγ为时频域滤波核长度，其值已在步骤4.1中初始化；

63、步骤4.5、计算时频域滤波结果ρp(i,j)的能量聚焦程度ep：

64、

65、其中，平均化的时频域滤波结果ρavg(i,j)利用时频域滤波结果ρp(i,j)计算得到：

66、

67、步骤4.6、令p＝p+1，并判断以下条件是否成立：

68、

69、其中，p为自项初始径向角的离散索引，p为步骤3.3预估得到的自项初始径向角个数，若上述条件成立，则进入步骤5；若条件不成立，则返回步骤4.3。

70、进一步的，步骤5中，采用如下方法以能量聚焦程度为依据，通过能量聚焦准则对自项初始径向角集合ψ进行拓展，得到自项拓展后径向角集合φafter，并通过φafter计算信号的最优时频分布ρ(i,j)，具体包括以下步骤：

71、步骤5.1、对以能量聚焦程度为依据，通过能量聚焦准则对自项初始径向角集合ψ进行拓展，得到自项拓展后径向角集合φafter，并通过φafter计算信号的最优时频分布ρ(i,j)的各参数进行初始化，具体包括以下参数的初始化：

72、①拓展自项径向角的角度步进初始化为：的正实数，其中，c为满足1≤c≤20的正整数；

73、②控制拓展自项径向角方向的标志flag1初始化为：flag1＝1；

74、③终止角度拓展的比例阈值rth初始化为：0.2≤rth≤0.8的正实数；

75、④自项拓展前径向角集合φbefore初始化为：φbefore＝ψ；

76、⑤最大值时频分布ρmax(i,j)初始化为：

77、ρmax(i,j)＝0,i＝0,1,…,l-1,j＝0,1,…,l-1

78、其中，i为时频域时间量离散索引，j为时频域频率量离散索引；

79、⑥最优时频分布ρ(i,j)初始化为：

80、ρ(i,j)＝0,i＝0,1,…,l-1,j＝0,1,…,l-1

81、其中，i为时频域时间量离散索引，j为时频域频率量离散索引；

82、步骤5.2、令自项初始径向角的离散索引p＝1；

83、步骤5.3、利用时频域滤波结果ρp(i,j)，计算最大值时频分布ρmax(i,j)：

84、ρmax(i,j)＝max(ρ(i,j),ρp(i,j)),i＝0,1,…,l-1,j＝0,1,…,l-1

85、其中，max(ρ(i,j),ρp(i,j))表示对ρ(i,j)和ρp(i,j)进行比较，取模值最大的元素值，ρ(i,j)为最优时频分布，其值已在步骤5.1中初始化，ρp(i,j)为步骤4.4中计算得到的时频域滤波结果，i为时频域时间量离散索引，j为时频域频率量离散索引，完成上述步骤后，令ρ(i,j)＝ρmax(i,j)；

86、步骤5.4、计算自项拓展前径向角

87、

88、其中，为步骤3.3中预估得到的第p个自项初始径向角；

89、步骤5.5、基于自项拓展前径向角得到自项拓展后径向角

90、

91、其中，flag1为控制拓展自项径向角方向的标志，其值已在步骤5.1中初始化，为拓展自项径向角的角度步进，其值已在步骤5.1中初始化；

92、步骤5.6、计算自项拓展后径向角的拓展后时频域滤波核γafter(aγ,bγ)：

93、

94、其中，vt为时频域滤波核的时间维长度参数，其值已在步骤4.1中初始化，vf为时频域滤波核的频率维长度参数，其值已在步骤4.1中初始化，tafter(aγ)和fafter(bγ)分别为自项拓展后径向角对应的拓展后时频域滤波核时间量和频率量，表示为：

95、

96、

97、其中，t(aγ)是时频域滤波核时间量，其值已在步骤4.1中初始化，aγ为滤波核时间量离散索引，aγ＝0,1,…,lγ-1，f(bγ)是时频域滤波核频率量，其值已在步骤4.1中初始化，bγ为滤波核频率量离散索引，bγ＝0,1,…,lγ-1，lγ为时频域滤波核长度，其值已在步骤4.1中初始化；

98、步骤5.7、利用拓展后时频域滤波核γafter(aγ,bγ)和信号的wigner-ville分布ρwvd(a,b)，按如下方式计算拓展后时频域滤波结果ρafter(i,j)：

99、

100、其中，i为时频域时间量离散索引，j为时频域频率量离散索引，表示对括号内所有可计算的数值矩阵中的每个元素求和，其中，a和b为满足以下条件的非负整数值：

101、

102、其中，l为包含整个信号持续时间所对应的采样点个数，其值已在步骤1中初始化，lγ为时频域滤波核长度，其值已在步骤4.1中初始化；

103、步骤5.8、计算拓展后时频域滤波结果ρafter(i,j)的能量聚焦程度eafter：

104、

105、其中，平均化的拓展后时频域滤波结果ρ′avg(i,j)利用拓展后时频域滤波结果ρafter(i,j)计算得到：

106、

107、步骤5.9、判断以下条件是否成立：

108、eafter＞ep·rth且且

109、其中，ep为步骤4.5中利用第p个自项初始径向角计算得到的能量聚焦程度，rth为终止角度拓展的比例阈值，其值已在步骤5.1中初始化，φbefore为自项拓展前径向角集合，其值已在步骤5.1中初始化，若上述条件不成立，进入步骤5.10；若条件成立，按如下方式计算自项拓展后径向角集合φafter与最大值时频分布ρmax(i,j)：

110、

111、ρmax(i,j)＝max(ρ(i,j),ρafter(i,j)),i＝0,1,…,l-1,j＝0,1,…,l-1

112、其中，φbefore是自项拓展前径向角集合，其值已在步骤5.1中初始化，为步骤5.5中计算得到的自项拓展后径向角，符号∪表示对符号两边的集合求并集。max(ρ(i,j),ρafter(i,j))表示对ρ(i,j)和ρafter(i,j)进行比较，取模值最大的元素值，ρ(i,j)为最优时频分布，其值已在步骤5.1中初始化，ρafter(i,j)为步骤5.7中计算得到的拓展后时频域滤波结果。完成上述步骤后，令φbefore＝φafter，ρ(i,j)＝ρmax(i,j)，并返回步骤5.5；

113、步骤5.10、判断以下条件是否成立：

114、flag1＞0

115、其中，flag1为控制拓展自项径向角方向的标志，其值已在步骤5.1中初始化，若条件成立，令flag1＝-1，返回步骤5.5；若条件不成立，令flag1＝1，并进入步骤5.11；

116、步骤5.11、令p＝p+1，并判断以下条件是否成立：

117、p＞p

118、其中，p为自项初始径向角的离散索引，p为步骤3.3预估得到的自项初始径向角个数。若条件成立，结束步骤5，得到待检测信号的最优时频分布ρ(i,j)，若不成立，返回步骤5.3。

119、有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

120、1、本发明的时频分析方法采用参数预估，能够有效降低计算量。根据步骤3，通过radon变换对信号模糊域分布进行检测，预估得到信号的自项初始径向角，并以此为基础进行之后的角度拓展，以拓展后的角度设计时频域滤波核。这种参数预估的方法克服了使用所有角度进行滤波核设计的传统时频滤波方法的局限性，仅使用代表信号自项特征的角度来进行滤波器的设计，降低了计算量。

121、2、本发明的时频分析方法采用了能量聚焦准则，能够对信号进行自适应的分析，并获得良好的分析效果。根据步骤5，在拓展角度的过程中，根据能量聚焦准则计算时频域滤波结果的能量聚焦程度，通过不需要人为控制的条件判决，能够自适应地获得所有能够反映信号自项特征的角度，并用于设计时频域滤波核。而通过能量聚焦准则得到的角度，都能较好地反映信号自项特征，因此滤波后的时频分布能够更好地降低噪声与交叉项对信号自项分布辨识的影响。本发明方法基于参数预估与能量聚焦准则对时频域滤波核的方向设计做出优化，减少了运算量，并保证了时频分析的质量，对线性调频与双曲调频信号都具备高分辨率的时频分析效果。