基于插值短时分数阶傅里叶变换- 变权拟合的线性调频信号参数估计

曹伟浩，姚直象，夏文杰，闫肃

(1.海军工程大学 电子工程学院, 湖北 武汉 430033;2.国家海洋技术中心漳州基地筹建办公室，福建 厦门 361000; 3.91388部队，广东 湛江 524000)

0 引言

线性调频(LFM)信号是主动声纳应用较广泛的一种信号形式，由于水声环境的复杂性及水中可利用的信号频带有限，对噪声污染的LFM信号进行瞬时频率估计变得非常重要。LFM信号高精度的瞬时频率估计，不仅在水下目标探测、水声通信多普勒补偿等方面扮演了重要角色，而且在雷达、通信等领域得到了广泛应用。此外，提高LFM信号参数估计的精度，对水下高精度目标识别也具有重要意义。

短时傅里叶变换(STFT)[1-2]和维格纳- 威利分布(WVD)[3-4]是两种常用的对LFM信号进行时频分析的方法，但是STFT对于宽带信号频率估计效果并不理想，WVD是双线性变换，时频分析存在交叉项问题。分数阶傅里叶变换(FRFT)对LFM信号具有较好的能量聚集性，利用此变换对LFM信号进行时频分析的应用前景广阔[5-10]。

本文针对参数估计问题，利用声纳、雷达等领域常用的一种频谱修正算法——Rife插值法，提出一种基于Rife插值短时FRFT(STFRFT)方法来估计LFM信号瞬时频率，通过变权拟合(VWSF)估计信号的初始频率和调频率，进一步推导了初始频率和调频率估计的克拉美- 罗限(CRLB)，分析了CRLB随参数变化的规律。

1 时频分析和参数估计方法原理

1.1 短时分数阶傅里叶变换

LFM信号作为水声信号处理领域中一种常用的信号，其表达式如下：

x(t)=Aexp {-j(2πf0t+πμt2)}，

式中：A为信号幅度；f0和μ分别为信号初始频率和调频率；t为信号持续时间。实际处理中，信号采用离散形式，故给出离散FRFT的表达式[10]为

(1)

式中：Aα=exp {-jπsgn(sinα)/4+jα/2}/|sinα|1/2，α=(p-1)π/2，p为FRFT最优阶数，0<|p|<2；x为量纲归一化后新的横坐标；Δx为量纲归一化宽度；X1(·)为信号的傅里叶变换；N为采样点数，|m0|≤N，|n0|≤N.

若对信号进行加窗处理，则由(1)式可得第i个短时窗内离散FRFT的表达式为

(2)

式中：i=1,2,…,K，K为窗移动的总次数，ΔK为移动步进；m=0,1,2,…,M-1；h(·)为长度为M的短时窗；αi=(pi-1)π/2，pi为第i个窗内FRFT的最优阶数。通常考虑到估计精度与计算量等因素，移动步进取窗长的0.1～0.3倍之间比较合适。

图1 LFM信号归一化功率谱密度

FRFT作为傅里叶变换的一种推广，直观上可以认为是Chirp基分解，非常适用于处理Chirp信号。对一个窗内的LFM信号分别作傅里叶变换和FRFT，得到归一化功率谱密度如图1所示。从图1中可以明显看出，经过傅里叶变换后频谱幅度在一定带宽内起伏变化，而经过FRFT后频谱幅度位于某一频率，表明FRFT对LFM信号具有更好的能量聚集性。

1.2 基于STFRFT的LFM瞬时频率估计

对信号进行瞬时频率估计时，在一定频率范围内，短时窗内的信号可近似认为是单频信号。STFT中通常取信号频谱幅度的均值作为一个窗内的频谱幅度，而实际频谱幅度高于此均值。因此，可以令第i个短时窗内FRFT幅度最大值对应的频率为该窗内LFM信号瞬时频率fi的估计值i，

i=kiΔf，

由于FRFT满足Parseval关系，可得

(3)

式中：Efi、Eti分别为信号频域和时域能量；Ni为第i个短时窗内时域采样点数；Nαi为第i个短时窗内分数阶域采样点数，Nαi=Nfisinαi=μT2sinαi，Nfi为第i个短时窗内频域采样点数，αi=piπ/2，pi为第i个短时窗内FRFT的最优阶数，μ为调频率，T为短时窗时长；Ai为第i个短时窗内信号的频谱幅度。由(3)式可以得到第i个短时窗内信号经过FRFT后的频谱幅度为

式中：fs为采样率。

对于离散频谱，考虑栅栏效应后的频谱幅度为

(4)

当信号频率不是频率分辨率的整数倍时，不可避免地会存在频谱泄露情况，例如1个主瓣内会存在2个峰值。采用较常用的Rife插值法进行频谱修正，该方法的频差表示[11]为

δi=

分别采用STFT、插值STFT(ISTFT)、STFRFT和插值STFRFT(ISTFRFT)4种方法估计LFM信号的瞬时频率，结果如图2所示。由图2可见，整体而言，STFRFT的方法较STFT方法更优，对于STFRFT方法而言，ISTFRFT方法比STFRFT方法估计的频率更加接近实际频率。

图2 LFM信号瞬时频率估计

在得到LFM信号每一时刻的瞬时频率后，利用曲线拟合方法估计其初始频率和调频率。最小二乘拟合(CLSF)是将观测值赋予相同的权值，使观测值与拟合值的误差平方和达到最小。当Rife插值估计的频率存在异常点时，将每次观测值的影响同等看待显然不合理，应该对误差平方和较大的项赋予较小的权值，误差平方和较小的项赋予较大的权值。因此采用VWSF方法估计信号的初始频率和调频率。

1.3 基于VWSF的LFM初始频率和调频率估计

式中：f′i=a1+a2t，a1、a2为2个待估计参量；Wi为第i个时刻估计瞬时频率的加权值。分别求RW(f)关于a1、a2的偏导并令其偏导为0，可解得2个参量的值[12]如下：

以上两个公式中的权值可由预先设定的精度控制实现。

分别利用CLSF和VWSF两种方法估计的瞬时频率结果如图3所示。从图3中可以看出，当估计频率出现异常点时，利用常规CLSF方法得到的曲线偏离实际频率，而VWSF方法能够很好地拟合实际频率。

图3 LFM信号瞬时频率拟合曲线

2 初始频率和调频率估计的CRLB研究

2.1 插值STFRFT频率估计的CRLB

在一定条件下，任何估计量都存在一个均方误差下限，即CRLB.对于未知确定参量f0和μ的均方误差等于估计量的方差，该估计为无偏估计。因此，求得估计量的方差限即为其CRLB[13].

本文采用Rife插值法，参照文献[14]中的方法，给出一个窗内插值STFRFT方法估计频率的均方误差为

(5)

假设估计的LFM信号瞬时频率为

f′i=f′0+μ′i+wi，

式中：f′0和μ′分别为估计的初始频率和调频率；wi为高斯白噪声，其方差为σ2.则瞬时频率的似然函数为

P(x|f′0,μ′)=

关于未知参量的Fisher信息矩阵为

式中：Pf′0f′0为函数P对f′0的2阶偏导；Pf′0μ′、Pμ′f′0为函数P对f′0和μ′的2阶偏导；Pμ′μ′为函数P对μ′的2阶偏导。由此可得

(6)

2.2 初始频率估计的CRLB分析

根据(6)式可以得到初始频率估计的CRLB为

(7)

在输入信噪比一定的情况下，取δ=0.5，初始频率的CRLB与滑窗次数、FRFT的最优阶数有关。对(7)式的K求偏导，有

(8)

考虑实际情况，当K大于1时(8)式恒小于0，表明初始频率的估计精度随滑窗次数增加而提高，即CRLB随滑窗次数的增加而减小。同理，对α求偏导，有

(9)

由于阶数p是以4为周期的参数，只需考察区间p∈[-2,2]即可，根据对称性，取阶数范围为p∈[0,2]，则α∈[0,π]。当α∈[0，π/2]时cosα>0，(9)式为正值；当α∈[π/2,π]时cosα<0，(9)式为负值。因此初始频率估计精度随阶数的增加先降低、后提高，即CRLB随阶数的增加先增大、后减小。当阶数为1时，初始频率估计性能最低，此时CRLB达到最大。当阶数不变时，初始频率估计性能与滑窗次数呈正比，但是滑窗次数的确定还需兼顾频率分辨率、计算量等因素。因此，分析CRLB的变化规律对窗长等参数选取具有一定的指导意义。

初始频率估计的CRLB仿真结果归一化后如图4所示。其中，图4(a)的阶数取值为[1，2]。由图4(a)可以看出CRLB随滑窗次数增加呈递减趋势。图4(b)的阶数取值为[0，2]，由图4(b)可以看出CRLB随阶数的增加先增大、后减小，阶数为1时估计性能最低。这一结论与理论分析结果一致。

2.3 调频率估计的CRLB分析

根据(6)式可以得到调频率估计的CRLB为

(10)

式中：τ=ΔK/fs，即短时窗移动步进对应时长。

在输入信噪比一定的情况下，取δ=0.5，移动步进取窗长的一半时，调频率估计的CRLB与滑窗次数、信号时长和FRFT的最优阶数有关。调频率的CRLB表达式为

(11)

式中：Tx为信号时长，为便于分析，取Tx=1 s，此时调频率估计的CRLB与滑窗次数、FRFT的最优阶数有关，分析方法与初始频率估计的CRLB分析相同。分别求(11)式关于K和α的偏导，有

(12)

(13)

当K>1时，(12)式的值恒小于0，故调频率的估计精度随滑窗次数的增加而提高，即CRLB随滑窗次数的增加而减小。(13)式的正负由cosα决定，变化规律与初始频率估计的CRLB情况一致。

对调频率估计的CRLB仿真结果归一化后如图5所示。同样地，图5(a)的阶数取值为[1，2]。由图5(a)可以明显看出，调频率估计的CRLB随滑窗次数增加而减小。图5(b)的阶数取值为[0,2]。由图5(b)可以看出，调频率估计的CRLB随阶数的增加先增大、后减小，阶数为1时估计性能最低。这一结论与理论分析结果一致。

图5 归一化调频率估计的CRLB

3 实验验证

3.1 仿真实验

仿真信号采用LFM信号，初始频率为600 Hz，调频率为800 Hz/s，采样率为4 096 Hz，脉宽为1 s，噪声为高斯白噪声，信噪比变化区间为-12～0 dB.经过5 000次蒙特卡洛仿真，得到估计的调频率均值及其均方根误差如图6所示；估计的初始频率均值及其均方根误差如图7所示。图6和图7中，VWSF-ISTFT表示利用插值STFT估计信号瞬时频率，然后通过VWSF估计参数；VWSF-ISTFRFT表示利用插值STFRFT估计信号瞬时频率，然后通过VWSF估计参数。

图6 调频率估计

从图6(a)可以看出，输入信噪比为-11 dB时，VWSF-ISTFRFT方法估计的调频率已经接近真值，与VWSF-ISTFT方法相比，所要求的输入信噪比降低了1 dB；相比于CLSF-ISTFRFT方法，VWSF-ISTFRFT方法所要求的信噪比降低了约3 dB.从图6(b)可以看出，输入信噪比为-11 dB时，VWSF-ISTFRFT方法估计的调频率精度已经接近CRLB;与VWSF-ISTFT方法相比，相同估计精度下，VWSF-ISTFRFT方法要求的输入信噪比降低了1 dB；达到相同的估计精度时，VWSF方法同CLSF方法相比，要求的输入信噪比降低了约3 dB.

图7 初始频率估计

从图7(a)可以看出，输入信噪比为-11 dB时，VWSF-ISTFRFT方法估计的初始频率已经接近真值，与VWSF-ISTFT方法相比信噪比降低了1 dB；相比于CLSF-ISTFRFT方法，VWSF-ISTFRFT方法所要求的信噪比降低了约3 dB.从图7(b)可以看出，相同估计精度下，VWSF-ISTFRFT方法估计初始频率所要求的输入信噪比较VWSF-ISTFT方法降低了1 dB，VWSF方法所要求的输入信噪比较CLSF方法降低了约3 dB.

针对被动声纳检测到的脉冲信号参数估计问题，为了证明本文方法的普适性，在被动声纳工作频率范围内随机选取参数进行仿真。根据文献[15]可知，现有的被动声纳工作的大致频率范围为0.1～10 kHz，调频信号通常作为窄带信号应用，脉宽由几十毫秒到几秒，在合理的参数取值范围内取3组参数进行仿真。信号均采用LFM信号，噪声为高斯白噪声，信噪比变化范围为-12～0 dB，其他仿真条件及得到的结果如下。

1)仿真条件1.采样率为2 kHz，初始频率为50 Hz，带宽为80 Hz，脉宽为4 s，经过5 000次蒙特卡洛仿真，得到调频率和初始频率的估计误差如图8和图9所示。

图8 调频率估计误差(仿真条件1)

图9 初始频率估计误差(仿真条件1)

2)仿真条件2.采样率为8 kHz，初始频率为2 kHz，带宽为100 Hz，脉宽为1 s，经过5 000次蒙特卡洛仿真，得到调频率和初始频率的估计误差如图10和图11所示。

图10 调频率估计误差(仿真条件2)

图11 初始频率估计误差(仿真条件2)

3)仿真条件3.采样率为20 kHz，初始频率为8 kHz，带宽为400 Hz，脉宽为0.4 s，经过5 000次蒙特卡洛仿真，得到调频率和初始频率的估计误差如图12和图13所示。

图12 调频率估计误差(仿真条件3)

图13 初始频率估计误差(仿真条件3)

从上述3个不同条件的仿真结果中均可以看出，输入信噪比为-12 dB时，VWSF-ISTFRFT方法仍有较高的估计精度，与基于STFT方法相比，VWSF-ISTFRFT方法的初始频率和调频率估计精度更高，估计性能更优；达到相同的估计精度时，VWSF-ISTFRFT方法所要求的输入信噪比更低。

3.2 海试验证

2017年12月在某地海试，发射信号为LFM信号，采样率为64 kHz，脉宽为128 ms，归一化后初始频率为0.96，信噪比为-8 dB.分别采用VWSF-ISTFT、CLSF-ISTFRFT和VWSF-ISTFRFT方法3种方法估计的瞬时频率如图14所示。从图14中可以看出，VWSF-ISTFRFT方法估计的瞬时频率最接近实际频率。将初始频率和调频率归一化，估计的参数如表1所示。从表1中可以看出，VWSF-ISTFRFT方法估计的值最接近真实值，CLSF-ISTFRFT方法估计参数误差最大。

图14 LFM信号瞬时频率估计

从海试数据的处理结果中可以看出，VWSF-ISTFRFT方法估计的瞬时频率最接近实际频率；从表1中数据可以看出，VWSF-ISTFRFT方法估计的参数精度最高，与CLSF方法相比，VWSF方法的参数估计精度由百分之一量级提高到千分之一量级，提高了一个数量级。

表1 估计参数及相对误差

注：括号内的值为参数估计相对误差。

4 结论

本文针对低信噪比下LFM信号参数估计问题，利用Rife插值STFRFT方法估计LFM信号的瞬时频率，通过VWSF估计初始频率和调频率，并推导分析了初始频率和调频率估计的CRLB影响因素。通过仿真分析和海试数据验证，得到主要结论如下：

1)插值STFRFT方法估计的LFM信号瞬时频率较插值STFT精度更高，有着更优的估计性能，估计精度相同时，对输入信噪比的要求前者较后者更低。

2)利用瞬时频率进行参数估计时，VWSF方法可以消除插值异常点的影响，与CLSF方法相比，得到的参数估计精度更高。

3)初始频率和调频率估计的CRLB受滑窗次数和FRFT阶数影响，随滑窗次数的增加而减小，随FRFT阶数的增加先增大、后减小，阶数为1时，CRLB达到最大值，此时参数的估计性能最低。