基于时频分析频率重心的跳周期估计方法与流程

本发明属于跳频信号侦察技术领域，具体涉及一种基于时频分析频率重心的跳周期估计方法。

背景技术：

跳频信号属于典型的非平稳信号，国内外学者大多采用时频分析方法得到清晰的时频图后，再估计其跳频参数。但目前的时频分析方法存在聚集性和交叉项的矛盾，当时频聚集性较高时，会产生较多的交叉项；当交叉项较少时，时频聚集性较低；无法达到二者的同时最优。同时，为适应信道传输，跳频信号多经过升余弦等低通成型滤波，基带信号频率波形发生变化。目前针对跳频信号的跳周期估计方法，多采用时频分析窗内的能量最大值所在频率作为该窗载频，对该载频序列进行快速傅里叶(fastfouriertransform,fft)变换，得到跳周期估计值；这种方法对于经升余弦低通滤波后的跳频信号的跳周期估计性能较差。

跳频通信由于其优良的低截获概率和抗干扰性，在军事通信对抗领域得到广泛应用；其中，跳频信号的参数估计，尤其是跳周期估计已成为研究的重点。跳频信号属于典型的非平稳信号，国内外学者大多采用时频分析方法得到清晰的时频图后，再估计其跳频参数。时频分析方法主要包括以短时傅里叶变换(shorttimefouriertransform，stft)、gabor变换、小波变换和s变换等为代表的线性变换和以魏格纳威尔分布(wigner-villedistribution,wvd)、伪wvd(pseudowigner-villedistribution,pwvd)、平滑伪wvd(smoothedpseudowigner-villedistribution,spwvd)及以cohen类为核函数的cohen类非线性变换；但两类分析方法存在聚集性和交叉项的矛盾。s.barbarossa针对单一跳频信号，应用wvd时频分析方法得到时频图，以该时频窗内时频能量最大值处的频率作为该窗载频，并给出了跳频信号的跳周期、跳频率和跳时的估计方法，在频点较少时有很好的估计性能。在此研究基础上，国内外众多学者为减弱wvd分布交叉项对信号参数估计的影响，采用pwvd、spwvd和cohen类等通过在时、频域平滑、加窗抑制非线性交叉项。基于线性变换的时频分析方法能够完全消除交叉项的影响，国内学者陈利虎等在《跳频信号的时频分析》中采用stft方法和张曦等在《基于短时哈特莱变换的跳频信号时频分析》中采用短时哈特莱变换(shorttimehartleytransform,stht)方法，但存在聚集性不高的问题。对于聚集性和交叉项的矛盾，部分学者采用组合时频分析和重排时频分析方法达到二者综合效果，但同时引入了计算量大的问题。在得到跳频信号时频图的基础上，大多采用该时频分析窗内的时频能量最大值所在频率作为该窗载频，对该载频序列进行快速傅里叶(fastfouriertransform,fft)变换，得到跳周期估计值；同时陈利虎在其学位博士论文《跳频信号的侦察技术研究》中指出针对慢跳频系统，以同一跳内不同窗的频率重心表征该窗的载波频率，平均其重心差值，得到跳周期的估计值。冯涛等在《一种组合时频分布在跳频信号参数估计中的应用》中通过一组带通滤波器将跳频信号的各个频率分量取出，分别计算每个信号分量的wvd，并将各wvd累加求和，得到新的时频分布和参数估值，其关键在于如何确定带通滤波器数并设计匹配各分量的带通滤波器。张朝阳等在《多跳频信号的盲分离与参数盲估计》中针对接收到的多个未知任何先验参数的跳频信号，提出一种先分离各个信号再分别进行时频分析来估计跳频参数的方法，首先采用特征矩阵联合近似最优(jade)算法分离跳频信号，再利用多窗口重叠的spwvd来估计出跳频信号的跳周期等参数。沙志超等在《基于稀梳重构的跳频信号时频分析方法》和《基于时频稀疏性的跳频信号时频图修正方法》中提出采用稀疏重构的方法得到跳频信号的时频图，继而进行参数估计，但计算量较大。同时，以上的时频分析针对的跳频信号模型都是理想的矩形跳频信号模型，但实际跳频信号的射频频谱不可避免的会受到基带信号影响。但沙志超在《基于时频稀疏性的跳频信号时频图修正方法》中只给出了包括基带信号在内的跳频信号模型，并未深入分析基带信号频谱对跳频信号参数估计的影响。

技术实现要素：

本发明所要解决的技术问题是提供一种能够大大减弱基带信号频谱影响的基于时频分析频率重心的跳周期估计方法。

为解决上述技术问题所采用的技术方案是：一种基于时频分析频率重心的跳周期估计方法，其包括如下步骤：

(1)计算经升余弦低通滤波后的超短波跳频信号sn(τ)的第n个时频分析时刻第k个频率间隔点的时频分析值stftn(n,k)；

在(式1)中，

s为跳频信号功率；th为跳周期；tr为跳周期的频率转换时间；t0为起跳时刻；fn为跳频瞬时频率；θ为跳频信号相位；δt为时间采样间隔；δf为频率采样间隔；n＝1,2,3,…；k＝1,2,3,…；g(τ-k(th-tr)-t0)为经基带成型滤波之后的等效低通信号，g^*(τ-i·δt)为g(τ-i·δt)的共轭；

在(式1)中，若基带成型滤波为升余弦低通成型滤波，则经升余弦低通成型滤波后的等效低通信号表达式为

那么(式1)为

在(式3)中，b为基带信号带宽，α为升余弦滚降因子；

(2)第n个时频分析时刻即为第n个时频分析窗，计算第n个时频分析窗的频率重心

在(式4)中，fs为时频采样率；n为时频分析窗长，即stft的时频采样点数；stftn²(n,k)为第n个时频分析时刻第k个频率间隔点的时频分析值，即能量值；

(3)计算第n个时频分析窗的频率重心和第n-1个时频分析窗的频率重心的差值，即相邻时频分析窗的频率重心差值δfn-1：

(式5)中，n＝2,3,…；

(4)计算一跳内包含的时频分析窗数目nh：

(式6)中，δf为跳频率间隔，nh初始值为0；

若相邻时频分析窗的频率重心差值δfn-1小于跳频间隔δf，则此相邻时频分析窗处在同一跳周期内；反之若相邻时频分析窗的频率重心差值δfn-1大于跳频间隔δf，则出现了频率跳变；

(5)计算跳周期th的估计值

(式7)中，为nh的统计平均值，t为窗时间长度，

进一步地，所述时频分析窗长n的取值范围是16，32，64，128，512，1024，…，32768。

进一步地，所述时频采样率fs的取值范围是200mhz到500mhz。

进一步地，所述时频采样率fs为250mhz。

本发明的有益效果是：本发明针对基于时频分析能量最大值的跳周期估计方法对经升余弦滤波的跳频信号估计性能下降问题，提出了基于时频分析频率重心的跳周期估计方法，本发明以时频分析窗的频率重心的能量值表征该窗的时频分析结果，克服了传统方法中时频分析聚集性与交叉项矛盾的缺陷，达到了二者的最优，并大大减弱了基带成型滤波对跳频信号的频谱扩展和平坦的影响，避免了能量最大值法对于基带成型滤波的载频估计发散问题，得到了较好的估计性能，本发明对于跳频信号的跳周期估计更加实用和准确，同时为跳频信号的侦察估计提供了参考和借鉴。

附图说明

图1为理想矩形基带信号的时域波形。

图2为经升余弦滤波后基带信号的时域波形。

图3为理想矩形基带信号与升余弦滤波后基带信号的频域波形。

图4为利用频率重心法对升余弦滤波跳频信号的时频分析结果。

图5为利用能量最大值法对升余弦滤波跳频信号的时频分析结果。

图6为不同时频分析估计方法的th估计性能对比。

具体实施方式

本发明针对基于时频分析能量最大值的跳周期估计方法对升余弦滤波后的跳频信号估计性能下降问题，提出基于时频分析频率重心的跳周期估计方法。首先给出跳频信号的数学模型，详细分析基带成型滤波对跳频信号的频谱影响和能量最大值估计方法性能下降的原因，并根据常用时频分析方法，提出基于时频分析频率重心的跳周期估计算法及步骤，最后进行仿真验证，得出对于经升余弦滤波的跳速较快的慢跳频信号，基于stft频率重心的跳周期估计方法更适用的结论。

对于经升余弦等基带成型滤波的跳频信号，其中心载频处的时频能量值不一定最大，以时频窗内能量最大值所在频率无法表征该窗的真实载频，因此本发明以时频分析的频率重心表征该窗的载频值，通过对比前后窗的频率重心差值与跳频率间隔，确定一跳内包含的时频窗数目；统计平均窗数目与每个窗时间长度的乘积，得到跳周期估计值。基于时频分析频率重心的跳周期估计算法步骤如下：

(1)计算信号sn(τ)的第n个时刻的时频分析值stftn(n,k)；

式中，s为跳频信号功率；th为跳周期，th的取值范围是0.5毫秒到100毫秒之间，对于超短波跳频信号th为5毫秒；tr为跳周期的频率转换时间，tr与电台的频率合成器有关，对于跳速较慢的信号tr为0.5毫秒；t0为起跳时刻，它由开机时间与跳周期关系决定，0＜t0＜th；fn为跳频瞬时频率，它由电台的使用频段决定，对于超短波电台，fn的取值范围是30mhz～88mhz；θ为跳频信号相位，0＜θ＜2π；δt为时间采样间隔，它由stft的采样率决定；δf为频率采样间隔，它由stft的采样率和取样点数决定；n＝1,2,3,…为第n个时频分析时刻；k＝1,2,3,…为第k个频率间隔点；g(τ-k(th-tr)-t0)为经基带成型滤波之后的等效低通信号，g^*(τ-i·δt)为g(τ-i·δt)的共轭；经升余弦低通成型滤波信号后的等效低通信号表达式为

第n个时刻的时频分析值stftn(n,k)表示为

式中，b为基带信号带宽，α为升余弦滚降因子。

(2)平面薄片的重心公式

设在xoy平面上有n个质点，它们分别位于点(x1,y1),(x2,y2),…,(xn,yn)处，质量分别为m1,m2,…,mn，由力学知识可知，该质点系的重心坐标为

其中，为该质点系的总质量；

在平面重心公式中，xi为第i块微元的x轴坐标i＝1,2,…,n；mi为第i块微元的重量；

借鉴平面重心公式，计算第n个时频分析窗的频率重心

式中，fs为时频采样率，n为时频分析窗长即时频采样点数，stftn²(n,k)为第n个时频分析窗第k个频率间隔点的能量值。

(3)计算第n-1个时频分析窗与第n个时频分析窗的频率重心差值

式中，n＝2,3,…

(4)根据频率重心差值δfn-1，计算一跳内包含的时频分析窗数目nh

式中，δf为跳频率间隔，nh初始值为0。

也就是说，若相邻时频窗的频率重心差值δfn-1小于跳频间隔δf，则此相邻窗处在同一跳周期内；反之若相邻时频窗的频率重心差值δfn-1大于跳频间隔δf，出现了频率跳变。

(5)对nh做统计平均，并求其与窗时间长度t的乘积，得到跳周期估计值

式中，为nh为统计平均值，窗时间长度t为

式中，n为stft的时频采样点数，一般取16，32，64，128，256，…；fs为时频采样率，对于超短波跳频的时频分析系统一般为250mhz，具体由系统性能决定。

需要说明的是，若基于其他时频分析频率重心法估计跳周期的话，只需将算法的stftn(n,k)值，换为相应的时频分析结果，算法后续处理不变。

跳频信号的时频分析：

跳频信号是载波随机跳变的信号，其数学模型为

其中，s为跳频信号功率，th为跳周期，tr为跳周期的频率转换时间，t0为起跳时刻，fk为跳频瞬时频率，θ为跳频信号相位，g(τ-k(th-tr)-t0)为经基带成型滤波之后的等效低通信号。

跳频信号的参数集为跳周期、频率转换时间、起跳时刻和瞬时频率，对跳频信号的参数估计，即是估计参数集{th,tr,fk,t0}，但th的正确估计是其他参数估计的前提。同时，跳频信号的参数估计依赖于跳频信号的时频图，但跳频信号的时频图会不可避免的受到基带信号的频谱影响；因为跳频信号的射频频谱，只是将基带信号的频谱进行了射频搬移，其形状完全由基带信号的频谱形状决定，若不设法消除基带信号的频谱影响，只是以时频分析窗能量最大值处的频率值作为该窗的载频值，对未经基带成型滤波的理想信号尚可，但无法正确估计经基带成型滤波之后的跳频信号。

跳频信号的基带频谱：

目前跳频通信常用的基带成型滤波器为升余弦滤波器，其时域模型为

其中，b为基带信号带宽，α为升余弦滚降因子。

升余弦成型滤波可理解为对基带码元内插和升余弦加权之后，进行矩形滤波，其时间响应由两部分组成：理想矩形成型滤波器的时域因子sinc(4bτ)和升余弦滚降的时域因子理想矩形基带信号的时域波形如图1所示；升余弦滤波后基带信号的时域波形如图2所示；理想基带信号与升余弦滤波后基带信号的频域波形如图3所示。

由图1-3可以看出，经过升余弦滤波后，时域波形出现畸变，频谱旁瓣大幅度下降。这是因为升余弦滤波结构通过对时域波形的加权叠加，将前后码元关联起来，破坏了码元的独立性，使时域波形出现畸变，并且引入了码间串扰(intersymbolinterference,isi)；以时域波形的串扰，换取减小理想矩形滤波旁瓣能量过大的问题。对于理想矩形信号，其频谱能量集中于零频(对射频信号而言，集中于载频)，所以现在针对跳频信号的时频分析结果多为山峰状(三维)或等高线状(二维)，能够以时频窗内能量最大值所在频率表征该窗载频；但实际跳频信号必然会经过基带成型滤波，经升余弦滤波后，基带频谱变得平坦，使得通带内能量最大值处不一定能对准零频(对射频信号而言，不一定能对准载频)，所以针对升余弦滤波后的跳频信号，通过对时频能量最大值所在频率做fft得到跳周期估计值的方法性能会大幅度下降。

常用时频分析方法：跳频信号是典型的非平稳信号，对于跳频信号的频谱，应同时从时域和频域观察；时频分析能够从时间和频率两个方面同时观测信号的能量密度大小，以及它们随着时间和频率的变化；所以时频分析作为跳频信号的分析工具越来越重要，下面简单介绍常用的时频分析方法。

短时傅里叶变换：stft的基本思想是在传统傅立叶变换的框架中，把非平稳信号看成是一系列短时平稳信号的叠加，而短时性则通过时域上的加窗来实现，并通过平移参数来覆盖整个时域。其定义为

其中，h^*(τ-t)为窗函数的共轭。

stft的含义可解释为：在时域用窗函数h^*(τ-t)截取信号s(τ)，对截取结果进行fft，即可得到在τ时刻该段信号的频谱情况；不断地移动窗函数h^*(τ-t)的中心位置，即可得到不同时该的频谱值；这些频谱值集合，即是stfts(τ,f)。

cohen类时频分析：cohen类时频分析把众多二次时频分布统一成一种形式，可改善stft等线性时频分析方法的时频聚集性，其定义为

其中，φ(ν,τ)为分布的核函数。

核函数φ(ν,τ)决定了cohen分布可能的形式及时频分析的性能；当φ(ν,τ)＝1时，cohen分布变成wvd分布；当φ(ν,τ)＝h(ν)时，h(ν)为平滑窗函数，此时cohen分布变成pwvd分布，即未加窗的wvd和窗函数在频域的卷积；当φ(ν,τ)＝h(ν)l(τ)时，相当于wvd在时域和频域同时做平滑滤波，变为spwvd分布。

因为cohen类时频分析采用二次非线性变换，在改善聚集性的同时引入了交叉项干扰，如何选择性能更优的核函数是减弱cohen类时频分析交叉项影响的重点。

性能对比：陈利虎在其学位论文《跳频信号的侦察技术研究》和熊良才等在《choi-williams分布参数优化及其应用》中以信息熵为指标定量分析了几种时频分析的聚集性和交叉项性能，如表1所示。

表1几种时频分析的信息熵

在表1中，n为分析时段内的采样点数，h，g为相应的窗长度。

由表1可知stft的熵值最大，其聚集性相对较差；对于wvd、pwvd和spwvd等时频分析方法，由于增加了时、频域滤波等相关处理，提高了其时频聚集性。

同时较好的时频分析通常会消耗更多的运算量，陈利虎在《跳频信号的侦察技术研究》和《跳频信号的时频分析》同时给出了各种时频分析的运算量比较，如表2所示。

表2时频分析的运算量比较

在表2中，n，h和p定义与表1一致。

由表2可知，对于相同长度的采样点数和时频窗长，wvd的运算量最小，但其交叉项干扰较大，当跳频频点和信号数量较多时，不能得到清晰的时频图。stft的运算量相对较小；经过增大时频分析窗长度，可得到聚集性较好的时频图。同时对于跳频信号，需要在一个跳周期时间内完成信号的时频分析和参数估计等操作，运算量是主要考虑因素；所以运算量较小的stft方法应用较多。

跳周期估计算法：

对于经升余弦等基带成型滤波的跳频信号，其中心载频处的时频能量值不一定最大，故而以时频窗内能量最大值所在频率可能无法表征该窗的真实载频，提出以时频分析的频率重心表征该窗的载频值；通过对比前后窗的频率重心差值与跳频频率间隔，确定一跳内包含的时频窗数目；统计平均窗数目与每个窗时间长度的乘积，得到跳周期估计值。基于时频分析频率重心的跳周期估计算法步骤如下：

(1)计算信号s(τ)的时频分析值stftn²(n,k)；

(2)借鉴平面重心公式，计算第n个时频分析窗的频率重心

其中，fs为时频采样率，n为时频分析窗长，stftn²(n,k)为第n个时间窗第k个频率分量处的能量值(即谱图值)。

(3)计算前后两个窗的频率重心差值

(4)计算一跳内包含的时频分析窗数目nh

式中，δf为跳频率间隔。该式意义为若相邻时频窗的频率重心差值δfn-1小于跳频率间隔δf，则此相邻窗处在同一跳周期内；反之，出现了频率跳变。

(5)统计平均nh与窗时间长度t的乘积，得到跳周期估计值

需要说明的是，若基于其他时频分析频率重心法估计跳周期的话，只需将stftn²(n,k)值换为相应的时频分析结果，算法后续处理不变。

仿真实验：

为验证跳周期估计算法的正确性，搭建simulink仿真实验，其仿真参数如表3所示。

表3仿真参数

实验1频率重心法与能量最大值法对升余弦滤波跳频信号的时频分析结果影响：

频率重心法对升余弦滤波跳频信号的时频分析结果如图4所示，能量最大值法对升余弦滤波跳频信号的时频分析结果如图5所示。可以看出，对升余弦滤波后的跳频信号而言，频率重心法比能量最大值法有更好的时频聚焦性。频率重心法以当前时频分析窗的频率重心所在频率作为该时频分析窗的时频分析结果，正是因为重心的相对稳定性，保证了相同跳周期内的每个时频分析窗的分析结果的稳定性。能量最大法以每个时频分析窗的能量最大处所在频率作为该时频分析窗的时频分析结果，且因为升余弦滤波后，跳频信号频谱变得平坦，频谱发生扩散(如图2和图3所示)，时频能量最大值处不一定能对准该窗的载频，所以在相同跳周期内的每个时频分析窗的分析结果出现不稳定状态，即频率扩散情况。频率重心法以该窗内的频率重心所在载频作为该窗的载频，可以减弱频谱扩展和平坦带来的影响。对频率重心法而言，由于统计因素也会出现个别窗的重心计算偏差，但该偏差远远小于跳频率间隔，在载频判断时可以归为一跳信号，不影响后续的th估计，同时从频率重心所在处的时频能量值也可看出，重心所在载频并不是能量最大值所在载频处，也证明了应用频率重心法比能量最大值法更适合于分析升余弦滤波后的跳频信号。

实验2不同时频分析估计方法的th估计性能：

不同时频分析估计方法的th估计性能对比如图6所示。

如图6所示为应用stft时频能量(即谱图)、wvd和spwvd的时频分析结果，然后应用频率重心法估计跳周期的归一化方差。由图6可以看出，当信噪比(signaltonoiseratio,snr)较小时，基于stft的频率重心法估计性能较好，在snr较大时，基于wvd和spwvd的频率重心法估计性能较好。这是因为基于时频分析频率重心法对跳周期的估计性能主要受到三个因素制约：一是噪声干扰频率重心的正确位置，二是时频分析的交叉项，三是时频分析的聚集性。在snr较小时，前两个是主要影响因素，而对于stft不受交叉项干扰影响，此时其估计性能较好；而spwvd通过时、频域加窗，消除了部分交叉项干扰，此时其估计性能次之；而wvd方法交叉项干扰严重，此时其估计性能较差。当snr较大时，后两个是主要影响因素，并且在snr充分大后，与噪声的交叉项干扰可忽略，此时估计性能主要是时频聚集性影响，并达到时频聚集性的极限界；stft的时频聚集性无法达到测不准原理的性能界，二次时频分析如wvd和spwvd可以达到该性能界，所以在snr较大时，基于stft频率重心的估计性能差于基于wvd和spwvd频率重心的估计性能。

实验3不同时频分析估计方法的运算开销：

由前面的基于时频分析频率重心的跳周期估计算法步骤可知，对于不同的基于时频分析频率重心法估计跳周期的算法，在得到时频分析结果之后的运算开销是一样，都经过了h次乘法，h+3n/h次加法。不同时频分析估计方法的运算量不同，主要取决于各种时频分析的运算量不同。由表2可知，wvd的运算量要远远小于spwvd的运算量，对于慢跳频系统(时频窗内采样点数小于跳周期内采样点数)，跳周期的采样点数n大于时频窗的采样点数h，所以对于stft的运算量远远小于wvd，这是因为为保证慢跳频系统的时间分辨率，h要远远小于n。所以对于跳速快的慢跳频系统，基于stft的频率重心估计跳周期更加适用。

本发明针对基于时频分析能量最大值的跳周期估计方法对经升余弦滤波的跳频信号估计性能下降问题，提出了基于时频分析频率重心的跳周期估计方法。该方法大大减弱了基带成型滤波对跳频信号的频谱扩展和平坦的影响，以某时频分析窗内的频率重心表征该窗的载频，避免了能量最大值法对于基带成型滤波的载频估计发散问题，得到了较好的估计性能；通过对比不同时频分析的估计性能和运算量，得到对于跳速较快的慢跳频系统，基于stft频率重心估计跳周期的算法更加适用的结论。所提算法和结论为经基带成型滤波的跳频信号及实际跳频信号的参数估计提供了参考和借鉴。

以上所述实施方式仅为本发明的优选实施例，而并非本发明可行实施例的穷举。对于本领域一般技术人员而言，在不背离本发明原理和精神的前提下对其所作出的任何显而易见的改动，都应当被认为包含在本发明的权利要求保护范围之内。