三角波调频信号参数快速精确估计方法

，

(机电动态控制重点实验室，陕西 西安 710065)

0 引言

近年来，国内外学者针对三角波调频信号的参数估计问题进行了大量的研究。文献[1]和文献[2]提出一种基于 Radon-Ambiguity (RAT)变换与分数阶傅里叶变换相结合的参数估计方法，该方法相比于采用二维搜索的Radon-Wigner变换，以较小的运算量实现了对信号初始频率的估计，但利用RAT变换估计信号调频斜率的运算量非常大。文献[3]和文献[4]采用高阶累积量与多相滤波器组相结合的方法实现对信号的参数估计，该方法无需先验知识且能较好地抑制高斯噪声，然而该方法的估计精度会随着多相滤波器组数的增加而下降。文献[5]提出一种分数阶傅里叶变换与聚类分析相结合的参数估计方法，该方法克服了信号尖峰必须高于噪声幅度的限制，可以实现低信噪比条件下信号的检测与参数估计，但该方法要求截获信号的起点位于调制周期的起点，因此不适合无先验知识的信号参数估计。文献[6]提出一种基于短时傅里叶变换，解线性调频和最大似然估计相结合的参数估计方法，估计精度高，但引入最大似然估计校正精度，运算量非常大。文献[7]和文献[8]提出了基于最佳频率分辨率窗宽的高斯窗短时傅里叶变换与相关解调相结合的参数估计方法，该方法实时性强，但并未使用最佳高斯窗进行短时傅里叶变换，且在进行相关解调时粗估了信号线性调频段的起点和终点，导致该方法参数估计精度不足。本文针对Radon-Ambiguity变换、解线性调频、最大似然估计及短时傅里叶变换等三角波调频信号参数估计方法均无法兼顾运算量小和估计精度高的问题，提出了基于信号特征的最佳高斯窗短时傅里叶变换与分段相关解调相结合的参数估计方法。

1 三角波调频信号参数估计方法

1.1 基于短时傅里叶变换与最大似然估计的三角波调频信号参数估计方法

对于三角波调频信号的参数估计问题，文献[6]提出了基于短时傅里叶变换和最大似然估计的参数估计方法，该方法首先对信号解线性调频并做短时傅里叶变换，求出中心频率，调制带宽和调制周期的粗略估计值，最后在时频平面进行最大似然搜索，求出上述参数的精确估计值。方法流程如图1所示。

1.2 基于短时傅里叶变换与相关解调的三角波调频信号参数估计方法

对于三角波调频信号的参数估计问题，文献[7]提出了基于最佳频率分辨率窗宽的高斯窗短时傅里叶变换与相关解调相结合的参数估计方法，该方法首先根据时频分析的目的，选择高斯窗作为短时傅里叶变换的窗函数，其次使用基于瞬时频率的窗宽倍增搜索方法求出最佳频率分辨率窗宽，然后对三角波调频信号进行短时傅里叶变换，求得信号的中心频率和调制带宽，再使用相关解调求得信号的调制周期。方法流程如图2所示。

高斯窗gs(t)的定义如下：

(1)

式(1)中，σ为高斯窗的标准差，σ的取值将对高斯窗的时频分辨率产生影响。由Heisenberg不确定性原理可知，提高高斯窗的时频分辨率，可以减小信号的参数估计误差。1.2节方法虽然求得了最佳高斯窗宽，但在对算法进行Matlab仿真时选择了Matlab自带的高斯窗函数，该函数使用固定的标准差σ1=2.5进行短时傅里叶变换，并未考虑信号波形特性与高斯窗标准差的关系，因此选择固定标准差σ1=2.5的高斯窗短时傅里叶变换并不能使窗内信号的时频聚集性达到最佳，这将导致参数估计精度不足。所以选择固定标准差σ1=2.5的高斯窗并非最佳高斯窗，要想使用最佳高斯窗，不仅需要求得最佳高斯窗宽，还需要根据信号波形特性求出高斯窗最佳时频聚集性标准差。

2 基于信号特征的三角波调频信号参数快速精确估计方法

本文提出了基于信号特征的最佳高斯窗短时傅里叶变换与分段相关解调相结合的三角波调频信号参数估计方法。该方法首先选择固定的标准差σ1=2.5，使用1.2节提出的方法求得三角波调频信号的调频斜率k，根据k求出高斯窗最佳时频聚集性标准差σbest，进而选择σbest并使用1.2节提出的基于最佳频率分辨率窗宽的高斯窗短时傅里叶变换求出信号的中心频率fc和调制带宽ΔF，其次利用分段相关解调确定三角波调频信号严格线性调频段的起点和终点，最后在严格线性调频段内使用相关解调估计出信号的调制周期T。方法流程如图3所示。

2.1 基于信号特征的最佳高斯窗短时傅里叶变换

基于信号特征的高斯窗最佳时频聚集性标准差推导步骤如下：

给定线性调频信号xLFM(t)：

xLFM(t)=ej2π(f0t+0.5kt2)

(2)

引入变量β，对式(2)配平方，整理得式(3)，

xLFM(t)=ej2π[f0t+0.5k(t-β)2-0.5kβ2+ktβ]

(3)

任选窗函数L(t)，对式(3)进行短时傅里叶变换，得式(4)，

(4)

对式(4)进行整理，得式(5)

(5)

定义函数γ(t)，

(6)

定义函数ρ(f)，

(7)

定义加窗的线性调频信号λ(t)：

λ(t)=ejπkt2L(t)

(8)

由ρ(f)的定义可知，函数ρ(f)是λ(t)的傅里叶变换，将式(8)代入式(6)可得，

(9)

根据傅里叶变换的时移和频移性质，整理式(9)可得，

γ(t)=e-j2πf′t·FT{ej2π[(f0+kt)t′]λ(t′)}=
{e-j2πf′t·FT[λ(t′)]}|f′=f-f0-kt

(10)

将式(7)代入式(10)可得，

γ(t)=e-j2π(f-f0-kt)t·ρ(f-f0-kt)

(11)

将式(11)代入式(5)可得，

STFTx(t,f)=e-jπkt2γ(t)=
e-jπkt2e-j2π(f-f0-kt)t·ρ(f-f0-kt)

(12)

则线性调频信号xLFM(t)的谱图SPECx(t,f)为：

SPECx(t,f)=|e-jπkt2e-j2π(f-f0-kt)t·
ρ(f-f0-kt)|2=|ρ(f-f0-kt)|2

(13)

由式(13)可知，谱图SPECx(t,f)相当于λ(t)频谱的模平方ρ(f)2沿直线f=f0+kt的排列。因此ρ(f)2的宽度决定了线性调频信号xLFM(t)谱图的时频聚集性。

将式(1)定义的高斯窗函数gs(t)代入式(7)的ρ(f)中，可得，

(4)

由于信号x(t)的α转角分数阶傅里叶变换公式为：

(15)

其中，

(16)

因此，高斯窗函数gs(t)的α转角分数阶傅里叶变换为：

(17)

比较式(17)与式(14)可知，当cotα=2πk且ucscα=2πf时，

(18)

对式(18)进行整理得，

(19)

取cotα=2πk且ucscα=2πf，进行三角函数代换，整理式(19)并代入式(18)可得，

(20)

整理式(20)可得，

(21)

整理式(21)可得，

(22)

由式(22)可知，ρ(f)2为高斯型函数，因此其宽度由z(σ)决定：

(23)

z(σ)min=k/2π

(24)

2.2 分段相关解调

分段相关解调具体步骤如下：

1)将区间[tmin,tmax]内的信号均分为r段，为计算方便，取r为4的整数倍，各段信号编号依次为1,2,…，r，为对每段信号进行相关解调，求得每段信号的调频斜率依次为ki,1≤i≤r。

2)计算各段信号调频斜率的均值：因区间[tmin,tmax]端点处的信号为非严格线性调频信号，为减小误差，使用区间中间编号为第0.25r+1小段到编号为第0.75r小段的信号计算每段信号调频斜率的均值ka：

(25)

3)计算各段信号调频斜率ki与调频斜率均值ka的绝对值差Δki：

Δki=ki-ka,1≤i≤20

(26)

综上，本文提出了基于信号特征的最佳高斯窗短时傅里叶变换与分段相关解调相结合的参数估计方法，该方法根据三角波调频信号的波形特性，通过理论推导给出了高斯窗短时傅里叶变换最佳时频聚集性标准差，选择该标准差进行短时傅里叶变换对信号的中心频率和调制带宽进行了参数估计；使用分段相关解调确定了信号线性调频段的起点和终点，求出了信号的严格线性调频区间，在严格线性调频区间内进行相关解调，估计出了信号的调制周期。

3 运算量分析与仿真验证

3.1 运算量分析

假定信号线性调频段内信号点数为M点，采样信号共包含c个线性调频段，以最佳频率分辨率窗宽qbest进行短时傅里叶变换，设步长为a，则初次估计时，短时傅里叶变换的运算量V1为：

(27)

初次估计时相关解调的运算量P1为：

(28)

二次估计时，依然使用最佳频率分辨率窗宽qbest进行短时傅里叶变换，则二次估计时短时傅里叶变换的运算量V2为：

V2=V1

(29)

假定将近似线性调频区间分为r段，令U=M/(2r)，则分段相关解调的运算量P2为：

(30)

因为窗宽倍增搜索方法和相关解调的运算量远小于短时傅里叶变换本身的运算量，而分段相关解调的运算量又小于相关解调的运算量，因此若采样信号包含c个线性调频段，则本文方法的总运算量V为：

(31)

对比1.1节方法，其运算量主要为解线性调频，短时傅里叶变换和最大似然估计的运算量，依然假定线性调频段内信号点数为M点，采样信号包含c个线性调频段，则解线性调频的运算量X为：

X=M(4+2lbM)

(32)

以窗宽q进行短时傅里叶变换，设步长为b，则短时傅里叶变换的运算量Y为：

(33)

最大似然估计的运算量Z为：

Z=M(M+4+2lbM)

(34)

1.1节方法总的运算量R为：

(35)

由以上分析可知，1.2节方法和本文方法的运算量均由M决定，而1.1节方法的运算量由M2决定，因此本文方法运算量与1.2节方法的运算量在同一数量级，远小于1.1节方法的运算量，适用于信号的实时处理。

3.2 实验仿真

为验证本文方法的参数估计精度，选择三组三角波调频信号，在相同条件下，分别使用1.1节方法，1.2节方法和本文方法对信号的中心频率，调制带宽，调制周期进行参数估计。

给定三组三角波调频信号：

信号1：中心频率fc=270 MHz，调制带宽ΔF=40 MHz，调制周期T=20 μs。

信号2：中心频率fc=290 MHz，调制带宽ΔF=30 MHz，调制周期T=30 μs。

信号3：中心频率fc=310 MHz，调制带宽ΔF=20 MHz，调制周期T=40 μs。

仿真条件设定：三组信号长度均为5个周期，采样率均为4倍采样，仿真时均加入高斯白噪声，令信噪比从0～20 dB，每个信噪比进行1 000次蒙特卡洛实验，均采用归一化均方误差(NRMSE)作为参数估计精度的误差衡量标准。

图4分别给出了不同信噪比下，采用1.1节方法，1.2节方法以及本文方法时，信号1、信号2、信号3中心频率的估计误差；图5分别给出了不同信噪比下，采用1.1节方法，1.2节方法以及本文方法时，信号1、信号2、信号3调制带宽的估计误差。图6分别给出了不同信噪比下，采用1.1节方法，1.2节方法以及本文方法时，信号1、信号2、信号3调制周期的估计误差。

以上仿真结果表明，在相同信噪比下，本文方法相比于1.1节方法，中心频率的估计误差下降了一个数量级，调制带宽的估计误差下降了50%左右，调制周期的估计误差下降了40%左右；本文方法相比于1.2节方法，中心频率的估计误差下降了一个数量级，调制带宽的估计误差下降了70%左右，调制周期的估计误差下降了50%左右。

4 结论

本文针对三角波调频信号的参数估计问题，提出了基于信号特征的最佳高斯窗短时傅里叶变换与分段相关解调相结合的参数估计方法。该方法根据三角波调频信号的波形特征确定了高斯窗最佳窗宽和最佳时频聚集性标准差，进而使用短时傅里叶变换估计出了信号的中心频率和调制带宽，其次使用分段相关解调确定了信号的严格线性调频段，并估计出了信号的调制周期。理论分析和仿真表明本文方法以较小的运算量大幅降低了三角波调频信号中心频率，调制带宽和调制周期估计的均方误差，实现了对三角波调频信号参数的快速精确估计。

参考文献:

[1]袁伟明，王敏，吴顺君.三角波调频信号的检测与参数估计[J].电波科学学报，2005，20(5)：594-597.

[2]WANG MSH, CHAIN AK. Linear frequency-modulated signal detection using time-frequency analysis [J]. IEEE Transactions on Signal Processing, 2011, 46(3)： 471-486.

[3]冯志红，赖涛，赵拥军.低信噪比下对称三角线性调频信号参数估计[J].电波科学学报，2012，27(3)：520-525.

[4]Stank L. Quadratic and higher order time-frequency analysis based on the short time Fourier transform [J]. Signal Processing and its Applications, 2014, 10(14), 581-588.

[5]刘锋，徐会法，陶然.基于FRFT的对称三角线性LFMCW信号检测与参数估计[J].电子与信息学报，2011，33(8)： 364-370.

[6]华旭东，臧小刚，唐斌.一种三角波调频信号参数估计实时系统的设计[J].信息技术，2007，25(7)： 1009-1018.

[7]燕天，邹金龙，马捷.基于短时傅里叶变换与相关解调的参数估计方法[J].探测与控制学报，2016，38(2)： 18-22.

[8]Jones D L. Signal dependent time-frequency analysis using Gaussian window [J]. Applied Signal Processing, 2014,12(5)：379-388.