一种基于群时延迭代变权拟合的双曲调频信号参数估计方法

本发明属于信号处理领域，具体涉及一种基于群时延迭代变权拟合的双曲调频信号参数估计方法。

背景技术：

1、双曲调频信号具有的多普勒不变性，使其广泛应用于主动声纳、水声通信等多个领域中。起始频率和周期斜率是能够表征双曲调频信号特性的两个基本参数，如果能够对它们进行精确的估计，就可以恢复所接收到的双曲调频信号，这在电子侦察和对抗等领域有着十分重要的意义。因此，起始频率和周期斜率的快速和精确估计是主动声纳和水声通信等领域的重要研究内容。

2、当前双曲调频信号参数的估计方法主要有：(1)最大似然法，该方法在理论上可以达到克拉美劳下界，但由于不存在解析解，计算复杂度较高；(2)基于时频分析和图像处理的方法，包括：①结合维格纳-威利分布和霍夫变换的维格纳-霍夫变换，该方法为实现较高精度的参数估计，需要较高的计算量；②利用短时傅里叶变换或拟最大似然法的分析，该方法的性能受到短时傅里叶变换窗长选择的影响；③利用短时分数阶傅里叶变换的分析，该方法的时频分辨率相互制约，且具有较高的计算复杂度；(3)频谱分析法，该方法计算量较小，但由于频谱容易受到噪声影响，算法精度受到限制。

技术实现思路

1、技术问题：本发明目的在于提供一种无需参数搜索，可在保证参数快速估计前提下，提高参数估计精度的双曲调频信号参数估计方法。

2、技术方案：为解决上述技术问题，本发明提出一种基于群时延迭代变权拟合的双曲调频信号参数估计方法，该方法包括以下步骤：

3、(1)获取待处理的双曲调频信号数据序列x(n),n＝0,1,…,n-1，n为检测到的双曲调频信号脉宽长度所对应的采样点个数；

4、(2)提取双曲调频信号归一化平滑幅度谱yn(l),l＝0,1...,n/2-1；

5、(3)利用直方图统计，估计所述双曲调频信号的信号带宽离散频点个数nh和归一化平滑幅度谱yn(l)的离散重心频率索引kg；

6、(4)依据信号带宽离散频点个数和离散重心频率索引，提取双曲调频信号归一化平滑幅度谱yn(l)的上升沿特征u(l)与下降沿特征d(l)；

7、(5)依据信号带宽离散频点个数和离散重心频率索引，以及上升沿特征和下降沿特征，估计双曲调频信号离散下限频率索引klb和离散上限频率索引khb；

8、(6)提取双曲调频信号在离散频率索引lsp范围内的群时延序列τ(lsp),lsp=klb,klb+1,…,khb；

9、(7)利用群时延序列的连续性，得到满足正常值条件的双曲调频信号群时延矢量τ0；

10、(8)对满足正常值条件的双曲调频信号群时延矢量τ0进行迭代变权拟合；

11、(9)依据最后一次迭代变权拟合的结果，得到双曲调频信号起始频率的估计值和周期斜率的估计值

12、进一步的，步骤(1)中，采用如下方法获取待处理的双曲调频信号数据序列x(n),n=0,1,…,n-1，具体包括以下步骤：

13、从传感器接收n个采样点的实时采集数据作为待处理的双曲调频信号数据序列x(n),n＝0,1,…,n-1，或从存储器中提取从检测到信号时刻起始的n个采样点的数据作为待处理的双曲调频信号数据序列x(n),n＝0,1,…,n-1，所述的n为检测到的双曲调频信号脉宽长度所对应的采样点个数，取值为n＝2a，a为大于等于5的整数。

14、进一步的，步骤(2)中，采用如下方法提取双曲调频信号归一化平滑幅度谱yn(l),l＝0,1...,n/2-1，具体包括以下步骤：

15、(2-1)对数据序列x(n),n＝0,1,…,n-1做离散傅里叶变换，得到双曲调频信号的离散傅里叶变换结果x(l)；

16、

17、其中，l为离散频率索引，j为虚数单位，即

18、(2-2)计算双曲调频信号的离散傅里叶变换结果x(l)的模值，得到双曲调频信号的幅度谱xa(l)；

19、xa(l)＝|x(l)|，l＝0,1,…,n/2-1；

20、(2-3)对双曲调频信号幅度谱的各平滑参数进行初始化，具体包括以下参数的初始化：

21、①双曲调频信号幅度谱最大平滑迭代次数e初始化为：e＞1的整数；

22、②双曲调频信号幅度谱平滑处理窗长ξ初始化为：3≤ξ≤15的奇数；

23、③双曲调频信号幅度谱平滑误差判决门限εe初始化为：0<εe<1的实数；

24、④双曲调频信号幅度谱平滑迭代次数e初始化为：e=1；

25、⑤第0次平滑处理后的双曲调频信号平滑幅度谱s0(l)初始化为：s0(l)＝xa(l),l＝0,1,…,n/2-1；

26、⑥第0次平滑处理后的双曲调频信号幅度谱平滑误差p(0)初始化为：p(0)＝0；

27、(2-4)对第e-1次平滑处理后的双曲调频信号平滑幅度谱se-1(l)进行平滑处理，得到第e次平滑处理后的双曲调频信号平滑幅度谱se(l)为：

28、

29、其中，ks为第l个离散频点所对应的平滑窗口内的离散频率索引；

30、(2-5)更新第e次平滑处理后的双曲调频信号幅度谱平滑误差p(e)；

31、

32、(2-6)判断是否满足以下迭代平滑条件

33、e＜e且

34、若条件成立，则令e＝e+1，并返回步骤(2-4)；否则进入步骤(2-7)；

35、(2-7)对第e次平滑处理后的双曲调频信号平滑幅度谱se(l)进行归一化处理，得到双曲调频信号归一化平滑幅度谱yn(l)：

36、

37、其中，max[se(l)]代表se(l)的离散频率索引l满足0≤l≤n/2-1范围内的se(l)的最大值。

38、进一步的，步骤(3)中，估计所述双曲调频信号的信号带宽离散频点个数nh和归一化平滑幅度谱yn(l)的离散重心频率索引kg，具体包括以下步骤：

39、(3-1)对估计信号带宽离散频点个数nh和离散重心频率索引kg的各参数进行初始化，具体包括以下参数的初始化：

40、①双曲调频信号归一化平滑幅度谱直方图统计起始值υstart初始化为：0≤υstart＜1的实数；

41、②双曲调频信号归一化平滑幅度谱直方图统计结束值υend初始化为：υstart＜υend≤1的实数；

42、③双曲调频信号归一化平滑幅度谱直方图统计子区间数量c初始化为：c＞5的整数；

43、④双曲调频信号归一化平滑幅度谱直方图统计的第c个子区间hc初始化为：

44、υstart+(c-1)(υend-υstart)/c≤hc＜υstart+c(υend-υstart)/c

45、其中，c为直方图统计子区间离散索引，c＝1,2,…,c；

46、⑤双曲调频信号归一化平滑幅度谱纯噪声判别门限ρnoi初始化为：0＜ρnoi＜0.45的实数；

47、⑥双曲调频信号归一化平滑幅度谱直方图统计子区间离散索引c初始化为：c＝1；

48、(3-2)对所述双曲调频信号归一化平滑幅度谱yn(l)的分布情况进行直方图统计，具体方法为：分别统计双曲调频信号归一化平滑幅度谱yn(l)的值落在各子区间h1,h2,…,hc内的离散频点个数h1,h2,…,hc；

49、(3-3)对双曲调频信号归一化平滑幅度谱yn(l)的值落在第c个子区间hc内的离散频点个数hc，判断hc＞ρnoi*n是否成立，若条件成立，则进入步骤(3-4)；否则进入步骤(3-5)；

50、(3-4)更新双曲调频信号归一化平滑幅度谱yn(l)的值落在第c个子区间hc内的离散频点个数hc＝0；

51、(3-5)判断c＜c是否成立，若条件成立，则令c＝c+1，并返回步骤(3-3)；否则进入步骤(3-6)；

52、(3-6)依据双曲调频信号归一化平滑幅度谱yn(l)的值落在第c个子区间hc内的离散频点个数hc，估计双曲调频信号归一化平滑幅度谱含噪信号判别门限ρsig：

53、

54、ρsig＝max[ρsig-0.1，0.69]

55、其中，代表hc的离散频率索引c满足1≤c≤c范围内的hc的最大值对应的c的值，max[ ]表示取两者中最大的一个；

56、(3-7)统计双曲调频信号归一化平滑幅度谱yn(l),l＝0,1...,n/2-1中满足ρsig≤yn(l)≤1的离散频点个数，即得到信号带宽离散频点个数nh，并搜索yn(l)满足ρsig≤yn(l)≤1的nh个离散频点的离散频率索引计算该nh个离散频点的离散重心频率索引kg：

57、

58、其中，int[]代表四舍五入取整运算，i为信号带宽范围内的nh个离散频点的离散索引。

59、进一步的，步骤(4)中，采用如下方法提取双曲调频信号归一化平滑幅度谱yn(l)的上升沿特征u(l)与下降沿特征d(l)，具体包括以下步骤：

60、(4-1)计算双曲调频信号归一化平滑幅度谱yn(l)的差值序列σa(l)：

61、

62、(4-2)对双曲调频信号归一化平滑幅度谱yn(l)的上升沿特征u(l)与下降沿特征d(l)进行初始化：

63、u(l)＝0,l＝0,1,...,n/2-1

64、d(l)＝0,l＝0,1,...,n/2-1

65、(4-3)对重计算双曲调频信号归一化平滑幅度谱的上升沿与下降沿特征的各参数进行初始化，具体包括以下参数的初始化：

66、①重计算双曲调频信号归一化平滑幅度谱的上升沿与下降沿特征的参考窗长krw初始化为：krw＝int[max[nh*α,kg*β]]，α为0.1≤α≤0.3的实数，β为0.01≤β≤0.1的实数，如果krw小于3，则令krw＝3；其中，int[]代表四舍五入取整运算，max[]表示取两者中最大的一个；

67、②重计算双曲调频信号归一化平滑幅度谱的上升沿与下降沿特征的离散频率索引lec初始化为：lec＝int[krw/2]+1；

68、(4-4)分别统计双曲调频信号归一化平滑幅度谱的差值序列σa(l)在第lec个离散频率索引预设范围内正值和负值的离散频点个数，即在lec-int[krw/2]≤l≤lec+int[krw/2]范围内，所有满足σa(l)＞0的离散频点个数cr(lec)和满足σa(l)＜0的离散频点个数cf(lec)；

69、(4-5)重计算双曲调频信号第lec个离散频点对应的归一化平滑幅度谱的上升沿特征u(lec)与下降沿特征d(lec)：

70、u(lec)＝σa(lec)cr(lec)

71、d(lec)＝σa(lec)cf(lec)

72、(4-6)判断lec＜n/2-int[krw/2]是否成立，若条件成立，则令lec=lec+1，并返回步骤(4-4)；否则，进入步骤(5)。

73、进一步的，步骤(5)中，采用如下方法估计双曲调频信号离散下限频率索引klb和离散上限频率索引khb，具体包括以下步骤：

74、(5-1)双曲调频信号上升沿与下降沿特征的搜索窗长ksw初始化为：ksw＝int[min[nh*γ,kg*μ]]，γ为0.9≤γ≤1.1的实数，μ为0.15≤μ≤0.25的实数，如果ksw小于3，则令ksw＝3；其中，int[ ]代表四舍五入取整运算，min[ ]表示取两者中最小的一个，nh和kg分别是步骤(3)中得到的双曲调频信号归一化平滑幅度谱yn(l)的信号带宽离散频点个数和离散重心频率索引；

75、(5-2)分别搜索上升沿特征u(l)的最大值与下降沿特征d(l)的最小值所对应的离散频率索引，作为离散下限频率索引klb和离散上限频率索引khb：

76、

77、

78、其中，和分别表示离散频率索引l在ω1＝{l|max[kg-ksw,1]≤l≤max[kg-1,1]}范围内搜索上升沿特征u(l)的最大值对应的离散频率索引，以及在ω2＝{l|min[kg+1,n/2-1]≤l≤min[kg+ksw,n/2-1]}范围内搜索下降沿特征d(l)的最小值对应的离散频率索引。

79、进一步的，步骤(6)中，采用如下方法提取双曲调频信号在离散频率索引lsp范围内的群时延序列τ(lsp),lsp＝klb,klb+1,…,khb，具体包括以下步骤：

80、(6-1)依据双曲调频信号的离散傅里叶变换结果x(l)，得到双曲调频信号的频谱相位角φ(l)：

81、φ(l)＝angle[x(l)]，l＝0,1,…,n/2-1

82、其中，angle[x(l)]代表取x(l)的相位角；

83、(6-2)计算双曲调频信号频谱相位角φ(l)的差值序列σp(l)：

84、

85、(6-3)对双曲调频信号频谱相位角解缠的各参数进行初始化，具体包括以下参数的初始化：

86、①双曲调频信号解缠后的频谱相位角差值序列初始化为：其中，lsp＝klb,klb+1,…,khb为频谱相位角解缠离散频率索引；

87、②双曲调频信号频谱相位角解缠最大迭代次数b初始化为：b＞1的整数；

88、③双曲调频信号频谱相位角解缠离散频率索引lsp初始化为：lsp＝klb；

89、④双曲调频信号频谱相位角解缠迭代次数b初始化为：b＝1；

90、(6-4)判断第lsp个离散频点解缠后的频谱相位角差值序列是否成立，若成立，则进入步骤(6-5)；否则，进入步骤(6-6)；

91、(6-5)更新第lsp个离散频点解缠后的频谱相位角差值序列

92、

93、(6-6)判断是否满足以下频谱相位角解缠迭代条件：

94、b＜b且

95、若条件成立，则令b＝b+1，并返回步骤(6-5)；否则令b＝1，并进入步骤(6-7)；

96、(6-7)判断lsp＜khb是否成立，若条件成立，则令lsp＝lsp+1，并返回步骤(6-4)；否则进入步骤(6-8)；

97、(6-8)依据双曲调频信号解缠后的频谱相位角差值序列提取群时延序列τ(lsp)：

98、

99、其中，δf＝fs/n为离散傅里叶变换的频率分辨率，fs为采样频率。

100、进一步的，步骤(7)中，采用如下方法得到满足正常值条件的双曲调频信号群时延矢量τ0，具体包括以下步骤：

101、(7-1)对判断群时延序列连续性的各参数进行初始化，具体包括以下参数的初始化：

102、①双曲调频信号群时延正常值判决门限η初始化为：0.1≤η≤0.3的实数；

103、②双曲调频信号群时延离散频率索引lτ初始化为：lτ＝klb+1；

104、③双曲调频信号群时延正常频点个数a初始化为：a＝0；

105、④双曲调频信号群时延正常频点离散频率索引集合p初始化为空集；

106、(7-2)判断第lτ个离散频点的群时延τ(lτ)是否满足以下正常值条件：

107、|τ(lτ)-τ(lτ-1)|≤ητ(lτ)且|τ(lτ)-τ(lτ+1)|≤ητ(lτ)

108、若条件成立，则令a＝a+1，第a个正常频点离散频率索引pa＝lτ，p＝p∪[pa]；其中，p∪[pa]表示将pa并入集合p中；

109、(7-3)判断lτ＜khb-1是否成立，若条件成立，则令lτ＝lτ+1，并返回步骤(7-2)；否则进入步骤(7-4)；

110、(7-4)依据双曲调频信号群时延正常频点离散频率索引集合p＝[p1 p2 … pa]，得到满足正常值条件的双曲调频信号群时延矢量τ0：

111、

112、其中，a＝1,2,…,a，a为a个双曲调频信号群时延正常频点的离散索引。

113、进一步的，步骤(8)中，采用如下方法对满足正常值条件的双曲调频信号群时延矢量τ0进行迭代变权拟合，具体包括以下步骤：

114、(8-1)对双曲调频信号群时延迭代变权拟合的各参数进行初始化，具体包括以下参数的初始化：

115、①双曲调频信号群时延最大迭代变权拟合次数q初始化为：q＞1的整数；

116、②双曲调频信号群时延迭代变权拟合误差判决门限εq初始化为：0＜εq＜1的实数；

117、③双曲调频信号群时延迭代变权拟合次数q初始化为：q＝1；

118、④第0次迭代变权拟合后的双曲调频信号群时延拟合误差初始化为：

119、⑤第0次迭代变权拟合中的双曲调频信号群时延权重矩阵w0初始化为：

120、其中，表示以为特征值的对角矩阵，

121、(8-2)计算第q-1次迭代变权拟合中的中间参数矩阵bq-1：

122、

123、其中，和分别是第q-1次迭代变权拟合中的中间斜率参数和中间截距参数，是频率辅助矩阵，是第pa个离散频点的频率；

124、(8-3)计算第q-1次迭代变权拟合中得到的群时延矢量的估计值τq-1；

125、

126、其中，是第q-1次迭代变权拟合中第a个正常频点的群时延的估计值；

127、(8-4)计算第q次迭代变权拟合中的a个正常频点的群时延估计误差

128、

129、(8-5)计算第q次迭代变权拟合中的a个正常频点的群时延的权重

130、

131、其中，代表的离散频率索引a满足1≤a≤a范围内的的最大值，代表的离散频率索引a满足1≤a≤a范围内的的最小值；

132、(8-6)得到第q次迭代变权拟合中的双曲调频信号群时延权重矩阵wq：

133、

134、(8-7)更新第q次迭代变权拟合后的双曲调频信号群时延拟合误差

135、

136、(8-8)判断是否满足以下群时延迭代变权拟合条件：

137、q＜q且

138、若条件成立，则令q＝q+1，并返回步骤(8-2)；否则进入步骤(9)。

139、进一步的，步骤(9)中，依据所述第q次迭代变权拟合中的中间斜率参数和中间截距参数采用如下方法得到双曲调频信号起始频率的估计值和周期斜率的估计值

140、

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

142、1、本发明充分利用了双曲调频信号幅度谱和噪声幅度谱的能量分布特征，基于直方图统计对双曲调频信号的带宽进行了初步估计，如步骤3所示，并利用了观测信号幅度谱的重心频率代替常用的峰值频率，从而提高了后续带宽精细估计的稳健性。

143、2、本发明充分利用了双曲调频信号频谱的瞬变特征，对双曲调频信号的带宽进行了精细估计，如步骤4、5中所示，所提取的瞬变特征参量在低信噪比时仍能较好地反映双曲调频信号的频域分布特性，提高了带宽估计的稳健性和精度。

144、3、本发明充分利用了群时延的连续性特征，保留了双曲调频信号已提取的群时延中的正常值，如步骤7所示，提高了双曲调频信号参数估计的稳健性和精度。

145、4、本发明依据迭代变权最小二乘拟合的方法，自适应地调整各频点处群时延所对应的权重，如步骤8所示，改善了群时延的拟合效果，从而进一步提高了双曲调频信号参数估计的稳健性和精度。