基于期望最大化算法的双音频率估计方法研究*

袁 恩，刘 鹏，陈隆亮，张文宇，常 青

(陆军工程大学，江苏 南京 210007)

0 引言

双音信号在多频连续波测距雷达[1]、无线传感器网络定位[2-3]、无线电信号导航[4]等方面具有广泛应用。文献[4]针对当前卫星导航信号易受干扰的问题，提出一种应用于航空导航的抗干扰测距导航方法。该方法在高速跳频信号的每一跳发送不同的双音信号，构建时频矩阵用于解算地面站与飞机之间的伪距。在上述应用中，都需要估计双音信号的频率，但是，文献[4]对双音信号频率估计方法提出了新的挑战：(1)高速跳频信号驻留时间短，使得信号采样时间较短；(2)为了提高伪距解算的性能，需要提高解算带宽，那么就需要时频矩阵中第一行的双音频率的间隔尽可能小。此外，在基于双音信号的无线传感器网络定位中，为了定位需要多个节点同时发送双音信号，不同节点发送的双音信号要保持一定间隔以防止相互影响。当频谱资源紧张时，双音信号频率间隔可能会比较小。

双音信号的频率估计问题是一个经典的问题，具有广泛的研究[5-9]。文献[5]给出了多音信号参数的近似的最大似然估计量，并利用基于离散傅里叶变换(Discrete Fourier Transform,DFT)的方法解算。但是，当音数大于1时，基于DFT的方法会使估计值存在偏差。虽然加窗可以降低这种偏差，但是却会降低频率的分辨率。基于旋转因子不变法(Estimating Signal Parameters Viarotational Invariance Techniques,ESPRIT)、多重信号分类(Multiple Signal Classification,MUSIC)的频率估计方法[7-8]具有高分辨率的特性，但是当采样样本较少时，协方差的估计误差会比较大；另外，在多音情况下ESPRIT等方法对频率估计的均方误差(Mean-Square Error,MSE)也无法达到克拉美罗下界(Cramer-Rao Lower Bound,CRLB)[7]。文献[7]提出NS-ESPRIT虽然能够使频率估计的MSE达到CRLB，但是，频率子带范围的选取对结果影响很大，而文献[7]并没有给出频率子带范围的有效选取方法。

文献[9]提出了一种高效的迭代算法，它利用期望最大化(Expectation Maximization,EM)算法结合IEEE 1057标准实现对多音实信号的参数估计。该方法首先基于DFT的方法得到多音信号参数的初估计；然后将多音信号的采样数据分解为多个单音采样数据；最后利用IEEE 1057标准采用最小二乘拟合的方法迭代计算各个单音的参数。上述过程中，多音信号采样数据的分解和各个单音参数的估计迭代进行，得到最终的估计值。当双音信号频率间隔足够大、采样次数足够多时，该方法频率估计的MSE可以达到CRLB。但是，文献[9]并没有考虑双音间隔较小、采样数据较短的情况，同时文献[9]中的方法是针对实信号的。

针对双音信号频率间隔小、采样数据少情况下估计双音信号频率面临的问题，本文在文献[9]的基础上提出基于EM算法的双音频率估计算法，主要改进在于以下几点：(1)EM算法中初值对结果有重要影响，针对双音间隔小的情况，选择具有高分辨率特性的MUSIC算法作为双音频率初值的获取方法；(2)针对样本数据少的问题，在MUSIC算法中采用修正的前向-后向线性预测(Modified Forward-Backward Linear Prediction,MFBLP)方法估计协方差[10]；(3)对于单音频率的估计，采用高效的二分搜索方法[11]。

本文提出的算法在EM算法的初估计时采用MUSIC算法，因此，在文中简称MEM算法。

1 信号模型

复双音信号可表示为：

(1)

其中，ai为双音信号的幅度，fi为双音信号频率，θi为双音信号相位。那么，双音信号的未知参数为：

(2)

双音信号经过采样频率为fs的模数转换器采样后，得到信号样本为：

0≤n≤N-1

(3)

s=[s(0)s(1)…s(N-1)]T

(4)

其中，w(n)为方差为σ2的高斯白噪声，Ts=1/fs为信号的采样间隔，N为样本数，NTs为总采样时间。文献[4]中，高速跳频信号单跳驻留时间短，令驻留时间为Tr，那么总采样时间NTs≤Tr，这使得在采样频率一定的条件下，样本数N较小。s为得到的样本向量。

双音信号可以表示成两个单音信号的和，则公式(3)可以表示为：

(5)

xi(n)=aiexp(j(2πfiTsn+θi))+wi(n)

(6)

xi=[xi(0)xi(1)…xi(N-1)]T

(7)

其中，xi表示单音信号的样本序列，wi(n)是噪声w(n)的分量，可表示为：

wi(n)=βiw(n)

(8)

文献[9]指出βi可以是任意值，但是必须满足β1+β2=1。此外，各个单音信号的未知参数可表示为Ai=[aifiθi]T。

2 MEM双音频率估计算法

2.1 算法描述

MEM双音频率估计算法大致可以分为3个步骤：

(1)估计未知参数A的初值；

第(2)步和第(3)步是一个迭代过程，即第(3)步的估计结果用于下一次迭代时第(2)步的估计。这种迭代算法称为EM算法[9]。

初值的估计对最终的结果具有重要的影响，文献[9]提出可用DFT计算初值，但是要求双音信号满足一定的条件，即双音频率的间隔Δf≥3/T，其中T表示信号采样时间。文献[9]并没有考虑双音频率间隔较小的情况，而文献[4]面临双音频率间隔较小的情况，因此需要选择合适的初值估计方法。本文在2.2节中对初值估计方法进行详细描述。

由于双音信号采样时引入了噪声，在式(3)中，高斯白噪声w(n)是未知的。在双音未知参数A的初估计后，令A的初估计为A(0)= [a1(0)f1(0)θ1(0)a2(0)f2(0)θ2(0)]T，那么，在第(2)步对于第j次迭代，可以得到噪声的估计为：

(9)

那么，此时单音信号的样本的估计值可表示为：

(10)

(11)

此时，得到两个单音信号的样本的估计值为：

(12)

通过上述方法将双音信号的样本值分解成两个单音信号的样本值，那么，在第(3)步可以采用单音信号的估计方法估计单音信号的参数。采用二分搜索的方法[11]估计单音信号的频率，单音信号相位和幅度可利用文献[12]给出的估计方法进行估计。

将第(3)步得到的两个单音信号的参数的估计值作为下一次迭代时的输入，在迭代时设置门限Δf，当满足下列条件时，迭代结束：

(13)

2.2 初值估计方法

算法第(1)步中的初估计值对算法的最终结果有重要的影响，特别是当双音信号的频率间隔比较小时体现得更为明显。文献[9]提出采用DFT的方法进行初值估计。文献[11]提出的二分搜索方法是在DFT的基础上进行迭代计算，其针对单音信号能够获得较高的估计精度，可以作为一种EM算法中估计初值的方法。但是，基于DFT的方法会使估计值存在偏移，特别是当双音频率间隔较小时，这种偏移会增大。

MUSIC算法具有高分辨率的特性。令双音信号的采样数据s的协方差为R，对R进行特征分解，特征值较小的N-2个特征向量构成的矩阵为UN，那么双音信号的频率估计值为下列函数的两个最大的峰值：

(14)

(15)

其中，[]H表示共轭转置。通常协方差R可通过计算样本的自相关矩阵进行估计，但是，当采样数据较少时协方差矩阵的估计会存在较大的误差。借助MFBLP方法对协方差的估计进行修正，可提高频率估计的精度。此时，协方差R的估计值为：

(16)

(17)

其中，m为MFBLP方法的阶数，[]*表示共轭。

2.3 双音频率估计的CRLB

将CRLB作为衡量估计方法性能的标准。单音信号频率估计的CRLB由文献[12]给出。如果不考虑双音间的影响，可以将双音信号看成两个独立的单音信号的叠加，那么，双音信号频率估计的CRLB为：

(18)

3 仿真分析

根据文献[4]中仿真的高速跳频信号的单跳驻留时间，在仿真中信号采样时间T为50 μs，基于此，文献[9]要求双音间的间隔大于3/T=60 kHz，而本文提出的MEM算法则不受这个限制。在仿真中，采样频率fs为250 kHz，那么样本数N为12，双音信号的幅度a1=a2=1，双音信号的初相θ1和θ2随机生成，双音信号频率f1=10 kHz，f2则选择不同的频率。迭代时的门限值Δf设置为1 Hz。

在仿真中，为了分析MEM算法的性能，对比三种不同方法的双音频率估计性能。第一种方法是本文提出的MEM算法。文献[10]指出当MFBLP方法中阶数m为3N/4时性能最佳，因此取m=9。第二种方法初值采用基于二分搜索的DFT方法进行估计[11]，其他步骤与第一种方法相同，称这种方法为初值采用二分搜索的EM算法。第三种方法采用MUSIC算法直接估计双音频率。

门特卡罗仿真次数为104次，结果如图1所示。在图1中，在不同信噪比(Signal-to-Noise Ratio,SNR)下，仿真三种频率估计算法的均方根误差(Root-Mean-Square Error,RMSE)，而将CRLB的平方根CRLB1/2作为衡量估计性能的标准。图1(a)中f2为30 kHz，双音的频率间隔为20 kHz。频率估计的RMSE存在门限效应，初值采用二分搜索的EM算法在这种情况下已经失效，其频率估计的RMSE都处于无信息区。二分搜索方法是基于DFT的，当双音频率间隔较小时，DFT方法可能会无法分辨双音信号。MUSIC算法具有高分辨特性，在上述情况下，它虽然有效并且门限效应的门限与MEM算法相当，但是，此时频率估计的RMSE与CRLB1/2相差较大。MEM算法频率估计的RMSE在高信噪比下明显优于MUSIC算法，并且与CRLB1/2非常接近。MEM算法频率估计的RMSE没有达到CRLB1/2的一个重要原因是当双音非常接近时，公式(18)不能准确表达CRLB，由于双音间的相互影响，实际的CRLB值比公式(18)的计算结果大。

图1 双音频率估计的均方根误差

图1(b) 中f2为40 kHz，双音的频率间隔为30 kHz。在这种情况下，初值采用二分搜索的EM算法频率估计的RMSE的门限明显高于其他两种方法。MUSIC算法与MEM算法相比，频率估计的RMSE略高于后者。MEM算法频率估计的RMSE与CRLB1/2更加接近。图1(c)进一步提高了f2的值将其调整为50 kHz，此时双音的频率间隔为40 kHz。这时，初值采用二分搜索的EM算法和MEM算法频率估计的RMSE几乎相当，而MUSIC算法频率估计的RMSE比前两者略高。MEM算法频率估计的RMSE已经达到CRLB1/2。可见，随着双音频率间隔的提高，三种算法的估计性能逐渐接近，但是在双音间隔较小的情况下，MEM算法在三者中的估计性能最好。

4 结论

本文针对双音信号间隔小、采样时间短时估计双音频率时面临的问题，提出基于EM算法的双音频率估计方法，称该方法为MEM算法。MEM算法在频率初估计阶段采用MUSIC算法，并且针对采样时间短的情况采用MFBLP方法估计协方差以提高MUSIC算法估计性能。通过仿真分析了MEM算法的性能，并与初值采用二分搜索的EM算法、MUSIC算法进行了比较，主要结果有以下几点：(1)当频率间隔为1/T时，即频率间隔较小时，MEM算法具有较好的频率估计性能，MUSIC算法频率估计的RMSE与MEM算法有较大差距，而初值采用二分搜索的EM算法完全失效；(2)随着双音频率间隔的增大，三者频率估计的RMSE逐渐接近。可见，在双音频率间隔较小的情况下，采用MEM可有效估计双音信号频率。