密级：无多线性调频信号瞬时频率估计迭代算法二炮工程大学士官学院 作者 于鹏鹏 黄向阳 艾名舜摘要：针对多线性调频信号的瞬时频率估计问题提出一种快速算法，该算法以特征子空间跟踪算法为基础，结合矩阵线性变换和多项式方程求根得到参数估计。

该算法的优点是计算量小，其计算量仅与短时傅里叶变换相当；频率分辨力较高；多信号情况下不存在交叉项问题；当多信号的功率差异达到14dB 时仍能有效估计瞬时频率。

由于采用了矩阵求逆的步骤，该算法在低信噪比环境下性能较差。

仿真实验显示在信噪比不低于6dB 时本文算法具有明显的优越性。

关键词：线性调频 瞬时频率 时频分析 一、引言线性调频 (Linear Frequency Modulation, LFM) 信号在雷达、声纳、通信等领域有着广泛的应用，由于瞬时频率随时间变化，LFM 信号具有非平稳特性，因此通常采用时频分析的方法对其进行分析及参数估计。

短时傅里叶变换是一种简单的时频分析方法，但是时频聚集性较差；Wigner-Ville 分布[1]（WVD ）的时频聚集性较好，但由于采用了二次型变换，在多LFM 信号情况下不可避免地存在交叉项，为信号参数估计造成了一定的困难；在Cohen 类时频分布[2]的框架下各种核函数被设计出来用于抑制交叉项，自适应核函数[3-4]的提出进一步提高了交叉项的抑制能力，然而性能较优的时频分析方法计算量也较大，因此在一定程度上较低了此类算法的实用性。

上述方法都是描述信号功率在时频平面上的分布，即信号的功率谱，其频率分辨率受限于信号时窗长度的倒数，这个限制被称为“瑞利限”。

超分辨算法利用信号特征子空间的正交性得到信号在频域上的“伪谱”，使有限长信号的频率分辨率能够突破“瑞利限”，从而获得更优的参数估计，但由于传统的超分辨频率估计算法的计算量较大，该类算法很少被用于估计非平稳信号参数。

本文提出一种基于子空间跟踪的信号瞬时频率估计算法，该算法利用数据投影实现信号特征子空间的跟踪，对特征子空间矩阵进行线性变换后得到多项式系数，进而利用多项式方程求根的方法获得信号瞬时频率的估计。

本文算法得到的是信号在时频平面上的 “伪谱”，不仅具有较好的时频聚集性，而且在多LFM 信号情况下不存在交叉项的问题，更重要的是，本文算法的计算量仅与短时傅里叶变换相当，因此是一种快速算法。

二、信号模型考虑一维时间序列S (t )由M 个调频信号线性叠加而成1()()(),1,2,...,Mm m m t A t t t T ==+=∑S s n(1)这里21()exp(2())2m m m t j f t k t π=-+s ，m =1,2,…,M ， A m 、f m 和m k 分别表示第m 个信号的幅度、起始频率和调频斜率。

T 表示有限长采样点数，设采样频率为f s ，测向无模糊范围不大于12sf 。

n (t )表示通道噪声，这里假设为零均值高斯白噪声，设等间隔采样，将N 个连续的采样点构成的向量称为一个快拍，N > M ，忽略噪声，t 0时刻的快拍向量0()t y 可以表示为[]0000022(1)1111122()(),(1),...,(1)(),(),...,()[(),(),...,()]m m t tM M Mj f t j f N t m m m m m m m m m t t M M t t t t t t N A t A t e A t e A t A t A t ππ=-∆--∆======--+⎡⎤=⋅⋅⎢⎥⎣⎦=⋅∑∑∑y S S S s s s s s s F(2)其中，F 是包含当前瞬时频率的矩阵，表达式为[]12TM =F F F F L(3)0022(1)[1,,,],1,2,...,m m j f t j f N t T m e e m M ππ-∆--∆==F L(4)t ∆是采样时间间隔，f m0表示第m 个信号在t 0时刻的瞬时频率。

以下将快拍()t y 作为采样单位，即每次采样得到一个数据向量，其原理如图1所示。

单快拍数据的协方差矩阵表示为()()()H t t t =⋅R y y ，这是一个秩为1的病态矩阵，直接对其进行特征值分解对测频没有帮助，但是在迭代模式的子空间跟踪当中采用单快拍协方差矩阵有利于降低运算量。

[ ]()t Y =()t S (1)t -S (2)t -S (1)t N -+S ()t S [ ](1)t -Y =(1)t -S (2)t -S (3)t -S ()t N -S图1 采样快拍示意图三、算法描述定理1：设G 为对称非负定的矩阵，N N C ⨯∈G ，其特征值满足120N λλλ≥≥>，每个特征值对应的特征向量表示为12N ,,...,u u u 。

{}nU 为N L ⨯维的正交矩阵序列，定义迭代式()orthnorm((1)) 1,2,...t t t =⋅-=U G U (5)“orthnorm”表示正交归一化，通常可用Gram- Schmidt 方法实现正交化。

如果初始矩阵0U 满足012[,,...,]T L U u u u 非奇异，则有：12lim ()[,,...,]L t t →∞=U u u u定理证明见文献[5]。

在该定理的基础上，文献[6]首次提出数据投影法(DPM)，其表达式为()orthnorm{(())(1)} 1,2,...t t t t μ=-⋅⋅-=U I R U (6)设信号数目M 已经事先估计得到，构造初始矩阵U (0)为()N N M ⨯-维随机矩阵，利用上述迭代，U (t )能够逼近信号噪声子空间的基向量矩阵U n ，从而为超分辨频率估计奠定基础。

这里μ被称为“步长因子”，当满足0<μ<<max 1/λ时迭代是稳定的，max λ是矩阵的R (t )的最大特征值。

MUSIC 算法[7]是超分辨谱估计的经典算法，根据MUSIC 算法的原理，瞬时频率矢量m F 与信号噪声子空间的基向量满足正交关系，即m n ⋅=F U 0 (7)U 是一个秩为N-M 的列满秩矩阵，可以认为其后N-M 行向量2U 线性无关，前M 行向量1U 可用后N-M 行向量线性表示，将矩阵U 分割为12⎡⎤=⎢⎥⎣⎦U U U ，根据线性变换可得11122--⎡⎤⋅=⋅=⎢⎥⎣⎦U U U U U I (8)线性变换不改变正交性，因此仍然有m ⋅=f U 0(9)采用上述线性变换的好处是可以在不影响参数估计结果的情况下去掉原DPM 算法中的正交化步骤，从而降低计算量。

去掉正交化的DPM 表达式为()(())(1) 1,2,...t t t t μ=-⋅⋅-=U I R U (10)此时，U 不能逼近噪声子空间基向量矩阵U n ，而是得到一个列向量线性无关的矩阵U 0，并且有0n =⋅U U Ω，Ω是可逆的旋转矩阵。

将0U 分割为01002⎡⎤=⎢⎥⎣⎦U U U ，则有111101021212----⋅=⋅=⋅U U U ΩΩU U U(11)可见，当采用线性变换时，非正交的矩阵U 0与噪声子空间基向量矩阵U n 是等价的，因此对于方程(9)而言，式(10)与式(6)也是等价的，但前者的计算量为O （N （N-M ）2）阶，而后者计算量仅为3N （N-M ）次复数乘法，达到了同类子空间跟踪算法计算量的下限，因此大大降低了运算量。

定义02m j ftz e π-∆=，瞬时频率矢量m F 可以写成21[1,,,...,]N m z z z -=F ，令ij U 表示U 的第j 行第i 列的元素，i =1,2,…,N ，j =1,2,…,N-M ，则方程(9)可以展开为一组多项式方程11111211()1000M M i i i MM i i i M N i i N M i z U z zU z z U z -=+-=--=⎧+⋅=⎪⎪⎪+⋅=⎪⎨⎪⎪⎪+⋅=⎪⎩∑∑∑n (12)原理上这N-M 个多项式方程是等价的，每个方程都有M 个有效根z m ，m =1,2,…,M 落在复平面的单位圆上，利用下式可以得到M 个信号瞬时频率的估计0()/2m m f angle z π=-(13)考虑到解多项式方程计算量与方程阶数成正比，因此可以只用第一个多项式方程求解参数。

本文算法主要步骤简述如下：1. 根据信号数目M 和采样快拍维数N 构造N ⨯L 初始矩阵U (0)，其中L=N-M ；2. 利用当前时刻采样快拍数据，根据式(10)计算U (t )，这一步骤的计算量约为3NL 次复数乘法和NL 次减法；3. 利用线性变换式(8)得到U ，并利用矩阵U 的列向量构造多项式方程，这一步骤的计算量约为ML 2+O (L 2)；4.在多项式方程组(12)中任选一个进行求解，选取M 个有效根，根据式(13)得到信号瞬时频率估计，这一步骤的计算量约为O (M 2+M )；5. 若信号采样结束，算法停止；否则转到步骤2。

四、数值仿真目前的时频表示除了短时傅里叶变换和 WVD 等少数方法外，其他方法都谈不上快速算法[9]，考虑到WVD 在多LFM 信号情况下存在严重的交叉项，与本文算法也没有可比性，因此本实验中只将短时傅里叶变换与本文算法作比较。

设两线性调频信号参数为(A 1 , f 1, 1k )=(1, 260MHz, -7.5⨯1013Hz/s 2)，(A 2 , f 2, 2k )=(0.5, 180MHz, 4⨯1013Hz/s 2)，采样频率f s =109Hz ，采样时长T =1u s 。

本文算法中设快拍长度为N =8，20.510μ-=⨯，作为对比，短时傅里叶变换的时窗长度为50n s ，两种方法得到的时频图如图2、图3所示。

123456789108time / nsF r e q u e n c y / H z图2 本文算法得到的信号时频图（SNR=6dB ） 图3短时傅里叶变换得到的时频图（SNR=6dB ）实验结果显示，两种算法得到的时频图都不存在交叉项，而计算效率大致相当，但本文算法的时频聚集性明显优于短时傅里叶变换。

另外，当两信号功率相差较大时，采用传统算法小功率信号往往被大信号所掩盖而难以估计其参数，对小信号的检测和参数估计之前需要滤除大信号，当两信号参数接近时，滤除大信号的同时小信号能量往往也被衰减了，从而增加了小信号检测的难度。

本文算法对信号功率的差异不敏感，在两信号功率相差14dB的情况下本文算法仍能有效估计两信号的瞬时频率（图4），而短时傅里叶变换中小信号几乎完全被淹没（图5）。

8time / nsFrequency/Hz图4 本文算法得到的信号时频图（SNR=12dB）图5短时傅里叶变换得到的时频图（SNR=12dB）本文算法的特点是计算效率高，信号时频聚集性好且不存在交叉项影响，但是在信噪比较低时性能不稳定，此时必须使其与降噪算法配合使用。