CN103795660A - 基于噪声近似处理的双阶段频率估计方法 - Google Patents

Abstract

本发明提供的是一种基于噪声近似处理的双阶段频率估计方法。该方法将频率估计分为频率粗估计和频率精估计两个阶段进行。在频率粗估计阶段，利用傅里叶变换对含有加性干扰噪声的信号进行频率估计，得到一个粗略的估计值。之后利用信号处理相关技术，构造一个以真实频率和粗估计频率的差值为频率的频差信号。在频率精估计阶段，利用噪声近似理论，对信号进行噪声近似处理，之后利用最小二乘估计算法，得到更为精确的频率估计值。本发明能够降低信噪比门限，提高频率估计的精度。

Description

基于噪声近似处理的双阶段频率估计方法

技术领域

本发明涉及的是一种信号处理的参数估计方法，属于无线通信系统的信号检测领域。

背景技术

在信号处理领域和无线通信系统中，接收信号的频率估计问题已经引起了人们的广泛注意。在无线信号的传输过程中，信号会受到反射、散射、干扰、多径和多普勒效应的影响，使得在信号接收机终端接收的信号已经和发射机发射的原始信号有很大的不同，因此只有从未知的信号中正确的估计出信号频率，才能从混合接收信号中解调出需要的有用信号。

最为常用而且简单有效的频率估计方法就是利用傅里叶变换（DFT）将时域信号变成频域信号，再根据其频谱特性估计其频率值。但是由于傅里叶变换的精度与信号的长度成正比，且受到栅栏效应和频率能量泄露的影响，使得傅里叶频率估计方法的精度不是很高，很难满足人们对频率估计高精度的要求。

图1给出了应用DFT进行频率估计得到的信号频谱图。如图1所示，利用DFT得到的信号频谱是离散的。假设信号的采样点数为N，则DFT方法的最小分辨率为

也就是说相邻的两条频谱线之间的间隔为

在图1所示的实验中，信号的实际频率为f=0.2361，该值位于DFT的两条谱线f₁=0.233和f₂=0.239之间，因此导致了使用DFT算法进行频率估计时，会产生一定的估计偏差。

为了提高傅里叶频率估计方法的精度，一些改进的傅里叶频率估计方法被人们提出。这些方法主要包括DFT插值法和DFT比值法，其主要思想是在离真实信号频率最近的两个频谱线之间插入新的谱线，或根据这两个谱线的值的比值对估计值进行修正。但是由于干扰噪声的影响，在信噪比较低的情况下这些方法可能会导致错误的插值区间判断和错误的谱线比值计算，因此不但不能降低原始的估计误差，反而会引入更大的修正误差。

其他一些提高信号频率估计的方法主要包括有子空间算法（MUSIC算法为代表）、最大概率密度（ML）算法和一些迭代DFT算法（IFEIFC算法为代表）。这些算法虽然能在一定程度上提高频率估计的精度，但是他们所需要的计算复杂度很高，因此并不适用于实际的信号接收机的信号处理过程。

在实际的接收信号中，信号处理所需要估计的频率信息存在于在有用信号的相位之中，而由于加性噪声的存在，使得人们无法直接从接收信号中提取有用信号的相位信息。噪声近似理论为人们处理信号噪声提供了一种有效的工具。设信号模型为：

r(n)=Aexp[j(φ₀+wnt_s)]+z(n)   n=1,2,...,N

其中，A是信号的幅度，φ₀是信号的初始相位，w是信号的角频率，t_s是采样周期，N是信号长度，z(n)是服从复高斯分布的加性噪声，其均值为0，方差为σ_z。

于是可得信号的信噪比为SNR_r=A²/σ_z。将信号r(n)用另一种形式描述为：

r(n)=[1+v(n)]Aexp[j(φ₀+wnt_s)]

其中， $v\left(n\right)=\frac{1}{A}z\left(n\right)\exp [-j({\mathrm{\φ}}_0+{\mathrm{wnt}}_s)],$ 其方差为var(v(n))=1/SNR_r。

令v(n)=v_I(n)+jv_Q(n)，则有

$1+v\left(n\right)={\{{[1+v_I\left(n\right)]}^2+v_Q^2\left(n\right)\}}^{1/2}\exp \left[j{\tan }^{-1}\frac{v_Q\left(n\right)}{1+v_I\left(n\right)}\right]$

当信噪比SNR_r>>1时，

$1+v\left(n\right)\≅\exp \left[j{\tan }^{-1}{v}_Q\left(n\right)\right]\≅\exp \left[{\mathrm{jv}}_Q\left(n\right)\right]$

最终信号r(n)可以近似的写成

$r\left(n\right)\≅\mathrm{Aexp}\left[j({\mathrm{wnt}}_s+{\mathrm{\φ}}_0+v_Q\left(n\right))\right]$

从上式中可以看到，经过噪声近似处理以后，原来的加性噪声等效为相应的相位噪声。

将上面的两种技术相结合，可以为频率估计开辟新的途径。首先用DFT进行频率粗估计，然后利用粗估计值进行信号重构和处理，再利用噪声近似原理提取含有粗估计频率偏差的相位信息，利用最小二乘算法，得到频率偏差的估计值。这样做一方面可以解决单纯的DFT算法精度不高的问题，提高频率估计的精度，另一方面又不会增加太多的复杂度使其可以在实际情况中得到应用。

发明内容

本发明的目的在于提供一种能提高精度且降低复杂度的基于噪声近似处理的双阶段频率估计方法。

本发明的目的是这样实现的：

步骤一：频率粗估计

（1）信号的接收

通信系统中的接收信号一般用含有噪声的复正弦信号进行表示，经过接收机终端采样后，实际处理的信号为：

x(n)=Aexp[j(φ₀+2πfnt_s)]+v(n)   n=1,2,...,N

其中，A是信号的幅度，φ₀是信号的初始相位，f是信号的频率，t_s是采样周期，N是信号长度，v(n)是服从高斯分布的加性噪声、其均值为0方差为σ²，频率粗估计阶段的任务是完成对接收信号x(n)的频率粗估计；

（2）DFT频率估计

在接收到信号x(n)后，对时域的接收信号x(n)进行傅里叶变换，得相应的频域信号为

$X\left(k\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}x\left(n\right)\exp (-j\frac{2\mathrm{\π}}{N}\mathrm{nk}),k=\mathrm{1,2},...,N$

信号的频谱信息为P(k)=|X(k)|²，根据信号频谱的特性，在信号真实频率处所对应的频谱值最大，通过寻找频谱值最大的点来估计信号的频率，设第m个点的频谱值最大，根据对应关系信号的频率估计值为

${\hat{f}}_0=\frac{2\mathrm{\π}}{N}m;$

（3）本地参考信号的构建

利用DFT对接收信号x(n)进行频率估计，得到了频率粗估计值

构造一个本地参考信号，所述本地参考信号以估计值

作为参考信号的频率，构建的本地参考信号为：

$y\left(n\right)=\exp \left(j2\mathrm{\π}{\hat{f}}_0{\mathrm{nt}}_s\right),n=\mathrm{1,2},...,N$

本地参考信号包含的是频率粗估计的信息；

（4）频差信号的构建

利用接收信号x(n)与本地参考信号y(n)，构造一个反映粗估计频差的频差信号，通过将接受信号x(n)与构建的本地参考信号y(n)共轭相乘，得到频差信号为：

z(n)=x(n)×y(n)^*

=Aexp[j(φ₀+2πΔfnt_s)]+v′(n)

n=1,2,...,N；

其中，

为粗估计的频率偏差、该值是由于DFT估计精度不高而导致的，v′(n)为加性噪声项，经过上述变换，得到一个以Δf为频率的频差信号z(n)，频差信号z(n)的频率信息即为在粗估计阶段中产生的频率估计偏差，频差Δf的值将在频率精估计阶段进行估计；

步骤二：频率精估计

（1）频差信号的叠加处理

将z(n)按照每M个采样点叠加在一起，得到

$\begin{array}{c}s\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}z\left(n\right)\\ ={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}\left\{\mathrm{Aexp}\right[j({\mathrm{\φ}}_0+2\mathrm{\π\Δ}{\mathrm{fnt}}_s)]+v^{\′}\left(n\right)\}\\ =A\frac{\sin \left(\mathrm{\π\Δ}{\mathrm{fMt}}_s\right)}{\sin \left(\mathrm{\π\Δ}{\mathrm{ft}}_s\right)}\exp \left\{j\right[{\mathrm{\φ}}_0+\mathrm{\π\Δf}(M-1)t_s\\ +2\mathrm{\π\Δ}{\mathrm{fpMt}}_s\left]\right\}+{\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right)\end{array}$

其中，为s(p)的幅度，P为s(p)的长度，其值为向下取整最接近

的整数， $v^{\′\′}\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right)$ 为加性噪声项，

经过叠加处理后得到的信号s(p)是一个信噪比较高的信号；

（2）信号噪声近似处理

引入噪声近似处理操作，将加性噪声转换为相位噪声，然后直接提取信号的相位信息，进而实现对频差Δf的快速高效估计；

当信号的信噪比SNR远大于1时，根据噪声近似理论，此时信号的加性噪声近似等价于相应的相位噪声，于是有

s(p)≈A′exp{j[φ₀′+2πΔfpT_s+e_p]}   p=1,2,...,P

其中，φ₀′=φ₀+πΔf(M-1)t_s，T_s=Mt_s，e_p为利用噪声近似理论得到的相位噪声项；

（3）最小二乘估计

经过噪声近似处理后，s(p)的近似信号中不再含有加性噪声项，直接提取相位，得到的相位信息为θ_p=φ₀′+2πΔfpT_s+e_p，定义向量Θ=[θ₀,θ₁,...,θ_P-1]^T，H₀=[1,1,...,1]^T，H₁=[0,2πT_s,...,2π(P-1)T_s]^T，H=[H₀,H₁]，X=[φ₀′,Δf]^T，E=[e₀,e₁,...,e_P-1]^T，于是得到相位的向量表达式为

Θ=HX+E

根据最小二乘估计算法，令估计误差最小，得到的X估计：

$\hat{X}=\underset{X}{\min }{\left|E\right|}^2=\underset{X}{\min }{|\mathrm{\Θ}-\mathrm{HX}|}^2$

对上式求导，并令导数为零，得到含有频差Δf信息的项X估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^T\mathrm{\Θ}$

从

中提取频差估计值

于是得到最终的信号频率估计为

最优的叠加长度的选择方法为：

定义s(p)经过噪声近似处理后的信号为

s′(p)=A′exp{j[φ₀′+2πΔfpT_s+e_p]}

定义一个比例参数G来衡量噪声近似处理对原始信号的改变情况：

$G=\frac{s\left(p\right)}{s^{\′}\left(p\right)}=\frac{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s\left]\right\}+v^{\′\′}\left(p\right)}{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s+e_p\left]\right\}}$

对上式求均方误差

$\mathrm{varG}=E\left\{{|\frac{s\left(p\right)}{s^{\′}\left(p\right)}-E[\frac{s\left(p\right)}{s^{\′}\left(p\right)}\left]\right|}^2\right\}=\frac{1}{\mathrm{MSNR}}$

该值varG是原始信号s(p)相对于近似信号s′(p)在均方意义上的改变量，因此修正的均方误差表达式为

$\begin{array}{c}E\left[{(f-\hat{f})}^2\right]=(1+\mathrm{varG})\frac{6}{{\left({2\mathrm{\πt}}_s\right)}^{2}N(N^2-M^2)\mathrm{SNR}}\\ =\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\end{array}$

通过求取上式的最小值，得到最优的叠加长度M_opt

$M_{\mathrm{opt}}=\arg \underset{M}{\min }\left\{\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\right\}$

对上式求导并令导数值为零，得最优叠加长度M_opt满足下式

${2\mathrm{SNRM}}_{\mathrm{opt}}^3+{3M}_{\mathrm{opt}}^2-N^2=0$

从上式中求出最优的叠加长度M_opt。

本发明提供了一种高精度且复杂度较低的频率估计方法。本发明弥补了DFT频率估计方法估计精度不高的缺陷，引入了一个频率精估计阶段，通过频率精估计过程对DFT估计算法的估计值进行修正，以提高算法的整体精度。

本发明寻求出一种新的精估计方法，不仅能够获得较高的估计精度，而且能有效的降低现有精估计算法的复杂度。

本发明研究了在信号叠加处理过程中叠加长度对系统的估计性能的影响，提出了针对本发明的精确的均方误差公式，并且基于最小均方误差的思想提出了最优的信号叠加长度。

本发明的目的主要技术要点体现在：

（1）频率粗估计

设接收信号为：

x(n)=Aexp[j(φ₀+2πfnt_s)]+v(n)   n=1,2,...,N

其中，A是信号的幅度，φ₀是信号的初始相位，f是信号的频率，t_s是采样周期，N是信号长度，v(n)是服从高斯分布的加性噪声，其均值为0方差为σ²。

利用傅里叶变换（DFT），将时域的接收信号x(n)变成频域信号X(k)：

$X\left(k\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}x\left(n\right)\exp (-j\frac{2\mathrm{\π}}{N}\mathrm{nk}),k=\mathrm{1,2},...,N$

利用频域信号X(k)可得信号的频谱为

由DFT的原理可知，在信号真实频率处所对应的频谱值最大，因此通过寻找频谱值最大的点来估计信号的频率。设第m个点的频谱值最大，根据对应关系可得信号的频率估计值为

${\hat{f}}_0=\frac{2\mathrm{\π}}{N}m$

由于DFT频率估计的精度不高，导致估计的频率

与信号的实际频率存在一定的频差。为了研究频差的特性，本发明构造了一个本地参考信号，该信号以估计值

作为参考信号的频率，构建的本地参考信号为：

$y\left(n\right)=\exp \left(j2\mathrm{\π}{\hat{f}}_0{\mathrm{nt}}_s\right),n=\mathrm{1,2},...,N$

为了研究频率粗估计的估计偏差，需要构建一个频差信号，该信号以f和

的偏差为频率。通过将接受信号x(n)与构建的本地参考信号y(n)共轭相乘，可得到频差信号为：

z(n)=x(n)×y(n)^*

=Aexp[j(φ₀+2πΔfnt_s)]+v′(n)

n=1,2,...,N

其中，为粗估计的频率偏差，v′(n)为加性噪声项。经过上述变换，得到了一个以Δf为频率的频差信号z(n)。

（2）频率精估计

与接收信号x(n)相比，频差信号z(n)的频率值Δf是一个极小量，利用这个极小的频率值，将频差信号z(n)进行叠加信号处理，可以大大提高信噪比。

将z(n)按照每M个采样点叠加在一起，可得到

$\begin{array}{c}s\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}z\left(n\right)\\ ={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}\left\{\mathrm{Aexp}\right[j({\mathrm{\φ}}_0+2\mathrm{\π\Δ}{\mathrm{fnt}}_s)]+v^{\′}\left(n\right)\}\\ =A\frac{\sin \left(\mathrm{\π\Δ}{\mathrm{fMt}}_s\right)}{\sin \left(\mathrm{\π\Δ}{\mathrm{ft}}_s\right)}\exp \left\{j\right[{\mathrm{\φ}}_0+\mathrm{\π\Δf}(M-1)t_s\\ +2\mathrm{\π\Δ}{\mathrm{fpMt}}_s\left]\right\}+{\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right)\end{array}$

其中，

为s(p)的幅度，P为s(p)的长度，其值为向下取整最接近

的整数， $v^{\′\′}\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right)$ 为加性噪声项。

通过叠加处理之后s(p)的信噪比为而初始的接收信号x(n)信噪比为

因此经过叠加处理之后的信号s(p)的信噪比将近提高了M倍。

对高信噪比的信号，可以进行噪声近似处理，将加性噪声变成相位噪声。在信号s(p)=A′exp{j[φ₀′+2πΔfpT_s]}+v″(p)中，由于加性噪声v″(p)的存在，想要直接获取A′exp{j[φ₀′+2πΔfpT_s]}的相位信息是行不通的。因此对s(p)进行进行噪声近似处理，可得

s(p)≈A′exp{j[φ₀′+2πΔfpT_s+e_p]}   p=1,2,...,P

其中，φ₀′=φ₀+πΔf(M-1)t_s，T_s=Mt_s，e_p为利用噪声近似理论得到的相位噪声项。通过噪声近似处理后，可以看到噪声项跑到了有用信号的相位信息中，信号s(p)中不再含有加性噪声项，此时就可以从近似信号中直接提取信号的相位信息了。

提取s(p)近似信号的相位为θ_p=φ₀′+2πΔfpT_s+e_p，定义向量Θ=[θ₀,θ₁,...,θ_P-1]^T，H₀=[1,1,...,1]^T，H₁=[0,2πT_s,...,2π(P-1)T_s]^T，H=[H₀,H₁]，X=[φ₀′,Δf]^T，E=[e₀,e₁,...,e_P-1]^T，于是得到相位的向量表达式为

Θ=HX+E

根据最小二乘估计算法，令估计误差最小

$\hat{X}=\underset{X}{\min }{\left|E\right|}^2=\underset{X}{\min }{|\mathrm{\Θ}-\mathrm{HX}|}^2$

对上式求导，并令导数为零，可得到含有频差Δf信息的项X估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^T\mathrm{\Θ}$

从

中提取频差估计值

于是得到最终的信号频率估计为

（3）最优的叠加长度选择

定义s(p)经过噪声近似处理后的信号为

s′(p)=A′exp{j[φ₀′+2πΔfpT_s+e_p]}

定义一个比例参数G来衡量噪声近似处理对原始信号的改变情况：

$G=\frac{s\left(p\right)}{s^{\′}\left(p\right)}=\frac{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s\left]\right\}+v^{\′\′}\left(p\right)}{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s+e_p\left]\right\}}$

对上式求均方误差

$\mathrm{varG}=E\left\{{|\frac{s\left(p\right)}{s^{\′}\left(p\right)}-E[\frac{s\left(p\right)}{s^{\′}\left(p\right)}\left]\right|}^2\right\}=\frac{1}{\mathrm{MSNR}}$

该值varG是原始信号s(p)相对于近似信号s′(p)在均方意义上的改变量。因此修正的均方误差表达式为

$\begin{array}{c}E\left[{(f-\hat{f})}^2\right]=(1+\mathrm{varG})\frac{6}{{\left({2\mathrm{\πt}}_s\right)}^{2}N(N^2-M^2)\mathrm{SNR}}\\ =\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\end{array}$

通过求取上式的最小值，可以得到最优的叠加长度M_opt

$M_{\mathrm{opt}}=\arg \underset{M}{\min }\left\{\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\right\}$

对上式求导并令导数值为零，可得最优叠加长度M_opt满足下式

${2\mathrm{SNRM}}_{\mathrm{opt}}^3+{3M}_{\mathrm{opt}}^2-N^2=0$

从上式中可以求出最优的叠加长度M_opt。

附图说明

图1是DFT频率估计的频谱示意图。

图2是频率粗估计阶段的流程图。

图3是频率精估计阶段的流程图。

图4是本发明的估计误差与叠加长度M的对应关系。

图5是本发明与传统的一些估计方法的性能比较图。

图6是在相同复杂度的情况下本发明的性能验证图。

图7是本发明的流程框图。

具体实施方式

下面结合附图对本发明作更详细的描述：

本发明主要包括三个部分：频率粗估计阶段（如图2所示）和频率精估计阶段（如图3所示）和最优叠加长度的设计。频率粗估计阶段主要包括有：信号的接收、DFT频率估计、本地参考信号的构建和频差信号的构建。频率精估计阶段主要包括有：频差信号的叠加处理、信号噪声近似处理和最小二乘估计。最优的叠加长度是根据最小均方误差准则得到的，它可以使本发明的均方误差最小。接下来详细介绍每个流程都是如何操作的。

步骤一：频率粗估计

图2给出了频率粗估计阶段的操作流程图，频率粗估计是本发明的基础，它的可靠性直接决定整体频率估计性能的好坏。下面对其进行描述。

101信号的接收

在实际的通信系统中，接收信号一般用含有噪声的复正弦信号进行表示，经过接收机终端采样后，实际处理的信号为：

x(n)=Aexp[j(φ₀+2πfnt_s)]+v(n)   n=1,2,...,N

其中，A是信号的幅度，φ₀是信号的初始相位，f是信号的频率，t_s是采样周期，N是信号长度，v(n)是服从高斯分布的加性噪声，其均值为0方差为σ²。频率粗估计阶段的任务是完成对接收信号x(n)的频率粗估计。

102DFT频率估计

在接收到信号x(n)后，利用简便快捷的傅里叶变换（DFT）对x(n)进行频率估计。DFT是一种把时域信号变换成频域信号的方法，通过观察频域信号的频谱特性，可以估计出信号的频率值。对时域的接收信号x(n)进行傅里叶变换，可得相应的频域信号为

$X\left(k\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}x\left(n\right)\exp (-j\frac{2\mathrm{\π}}{N}\mathrm{nk}),k=\mathrm{1,2},...,N$

于是信号的频谱信息为P(k)=|X(k)|²，根据信号频谱的特性，在信号真实频率处所对应的频谱值最大，因此通过寻找频谱值最大的点来估计信号的频率。设第m个点的频谱值最大，根据对应关系可得信号的频率估计值为

${\hat{f}}_0=\frac{2\mathrm{\π}}{N}m$

103本地参考信号的构建

利用DFT对接收信号x(n)进行频率估计，得到了频率粗估计值

由于DFT频率估计的精度不高，导致估计的频率与信号的实际频率存在一定的频差，因此需要进一步研究频率粗估计值

与信号频率f的偏差有多大。本发明构造了一个本地参考信号，该信号以估计值

作为参考信号的频率，构建的本地参考信号为：

$y\left(n\right)=\exp \left(j2\mathrm{\π}{\hat{f}}_0{\mathrm{nt}}_s\right),n=\mathrm{1,2},...,N$

本地参考信号包含的是频率粗估计的信息。

104频差信号的构建

利用接收信号x(n)与本地参考信号y(n)，可以构造一个反映粗估计频差的频差信号。通过将接受信号x(n)与构建的本地参考信号y(n)共轭相乘，可得到频差信号为：

z(n)=x(n)×y(n)^*

=Aexp[j(φ₀+2πΔfnt_s)]+v′(n)

n=1,2,...,N

其中，

为粗估计的频率偏差，该值是由于DFT估计精度不高而导致的，v′(n)为加性噪声项。经过上述变换，得到了一个以Δf为频率的频差信号z(n)。频差信号z(n)的频率信息即为在粗估计阶段中产生的频率估计偏差，频差Δf的值将在频率精估计阶段进行估计。

步骤二：频率精估计

图3给出了频率精估计阶段的操作流程图。精估计阶段是为了估计在频率粗估计阶段产生的频差Δf。具体的操作步骤将在下面进行描述。

105频差信号的叠加处理

在频率粗估计阶段，得到了一个新的含有加性噪声的频差信号z(n)。与接收信号x(n)相比，频差信号z(n)的频率值是一个极小量，利用这个极小的频率值，将频差信号z(n)进行叠加信号处理，可以大大提高信噪比。

将z(n)按照每M个采样点叠加在一起，可得到

$\begin{array}{c}s\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}z\left(n\right)\\ ={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}\left\{\mathrm{Aexp}\right[j({\mathrm{\φ}}_0+2\mathrm{\π\Δ}{\mathrm{fnt}}_s)]+v^{\′}\left(n\right)\}\\ =A\frac{\sin \left(\mathrm{\π\Δ}{\mathrm{fMt}}_s\right)}{\sin \left(\mathrm{\π\Δ}{\mathrm{ft}}_s\right)}\exp \left\{j\right[{\mathrm{\φ}}_0+\mathrm{\π\Δf}(M-1)t_s\\ +2\mathrm{\π\Δ}{\mathrm{fpMt}}_s\left]\right\}+{\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right)\end{array}$

其中，

为s(p)的幅度，P为s(p)的长度，其值为向下取整最接近

的整数， $v^{\′\′}\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right)$ 为加性噪声项。

由于Δf为极小量，于是有 $A^{\′}=A\frac{\sin \left(\mathrm{\π\Δ}{\mathrm{fMt}}_s\right)}{\sin \left(\mathrm{\π\Δ}{\mathrm{ft}}_s\right)}\≈\mathrm{MA},$ 另外由 $v^{\′\′}\left(p\right)={\mathrm{\Σ}}_{n=\mathrm{pM}}^{(p+1)M-1}{v}^{\′}\left(n\right),$ 可得v″(p)的方差为Mσ²。可以看到通过叠加处理之后s(p)的信噪比为而初始的接收信号x(n)信噪比为

因此可以看到与接收信号x(n)相比，经过叠加处理之后的信号s(p)的信噪比将近提高了M倍。由于M的取值一般是一个比较大的整数，所以经过叠加处理后得到的信号s(p)是一个信噪比较高的信号。

106信号噪声近似处理

在得到了一个高信噪比的信号s(p)后，可以采用噪声近似处理操作，将加性噪声转化为相位噪声，大大简化后续的对于频差Δf的精估计

在信号s(p)=A′exp{j[φ₀′+2πΔfpT_s]}+v″(p)中，由于加性噪声v″(p)的存在，想要直接获得有用信号的相位信息是行不通的。而一旦获取了信号的相位信息就可以快速简便的完成对频差Δf的估计。为此本发明引入了噪声近似处理操作，将加性噪声转换为相位噪声，然后可以直接提取信号的相位信息，进而实现对频差Δf的快速高效估计。

当信号的信噪比SNR远大于1时，根据噪声近似理论，此时信号的加性噪声可以近似等价于相应的相位噪声，于是有

s(p)≈A′exp{j[φ₀′+2πΔfpT_s+e_p]}   p=1,2,...,P

其中，φ₀′=φ₀+πΔf(M-1)t_s，T_s=Mt_s，e_p为利用噪声近似理论得到的相位噪声项。通过噪声近似处理后，可以看到噪声项跑到了有用信号的相位信息中，信号s(p)中不再含有加性噪声项，此时就可以从近似信号中直接提取信号的相位信息了。

107-108最小二乘估计

经过噪声近似处理后，s(p)的近似信号中不再含有加性噪声项，可以直接提取相位，得到的相位信息为θ_p=φ₀′+2πΔfpT_s+e_p。定义向量Θ=[θ₀,θ₁,...,θ_P-1]^T，H₀=[1,1,...,1]^T，H₁=[0,2πT_s,...,2π(P-1)T_s]^T，H=[H₀,H₁]，X=[φ₀′,Δf]^T，E=[e₀,e₁,...,e_P-1]^T，于是得到相位的向量表达式为

Θ=HX+E

根据最小二乘估计算法，令估计误差最小，可以得到的X估计：

$\hat{X}=\underset{X}{\min }{\left|E\right|}^2=\underset{X}{\min }{|\mathrm{\Θ}-\mathrm{HX}|}^2$

对上式求导，并令导数为零，可得到含有频差Δf信息的项X估计值为：

$\hat{X}={\left(H^{T}H\right)}^{-1}{H}^T\mathrm{\Θ}$

从中提取频差估计值

于是得到最终的信号频率估计为

最优的叠加长度

定义s(p)经过噪声近似处理后的信号为

s′(p)=A′exp{j[φ₀′+2πΔfpT_s+e_p]}

定义一个比例参数G来衡量噪声近似处理对原始信号的改变情况：

$G=\frac{s\left(p\right)}{s^{\′}\left(p\right)}=\frac{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s\left]\right\}+v^{\′\′}\left(p\right)}{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s+e_p\left]\right\}}$

对上式求均方误差

$\mathrm{varG}=E\left\{{|\frac{s\left(p\right)}{s^{\′}\left(p\right)}-E[\frac{s\left(p\right)}{s^{\′}\left(p\right)}\left]\right|}^2\right\}=\frac{1}{\mathrm{MSNR}}$

该值varG是原始信号s(p)相对于近似信号s′(p)在均方意义上的改变量。因此修正的均方误差表达式为

$\begin{array}{c}E\left[{(f-\hat{f})}^2\right]=(1+\mathrm{varG})\frac{6}{{\left({2\mathrm{\πt}}_s\right)}^{2}N(N^2-M^2)\mathrm{SNR}}\\ =\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\end{array}$

通过求取上式的最小值，可以得到最优的叠加长度M_opt

$M_{\mathrm{opt}}=\arg \underset{M}{\min }\left\{\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\right\}$

对上式求导并令导数值为零，可得最优叠加长度M_opt满足下式

${2\mathrm{SNRM}}_{\mathrm{opt}}^3+{3M}_{\mathrm{opt}}^2-N^2=0$

从上式中可以求出最优的叠加长度M_opt。

性能分析

结合附图对本发明的性能进行分析。

（1）分析叠加长度对本发明性能的影响

图4给出了本发明的估计误差与叠加长度M的对应关系。从图4中可以开看到，本发明的估计误差不是随着叠加长度M的增加或减小呈现单调性变化，不是M的值越大估计误差就越小，也不是M的值越小估计误差就越小。当信噪比较低的时候，如图4中的SNR=-13dB时，估计误差与M值的关系就更为明显。可以看到估计误差与M值形成了近似的开口向上的抛物线形态，因此存在一个最优的叠加长度M_opt使得估计误差最小。

为此本发明经过研究提出了一种寻找最优叠加长度M_opt的方法。当信号的信噪比远大于1时，采用噪声近似处理之后的频率估计方法得到的均方误差为：

$E\left[{(f-\hat{f})}^2\right]=\frac{6}{{\left({2\mathrm{\πt}}_s\right)}^{2}N(N^2-M^2)\mathrm{SNR}}$

其中，E[]表示求平均值操作。从上式中，可以看到估计误差应该是随着M的值减小而减小，但是通过实验仿真可以看到实际情况并非如此。本发明经过深入研究，发现上面的均方误差表达式存在欠缺，它忽略了噪声近似处理阶段的噪声近似误差。虽然当信噪比很高时，噪声近似引入的误差很小，但是观察上面的均方误差表达式，同样可以看到由于N>>M，因此相对于N²来说，M²也是个很小的量。

为此，本发明给出了一个更加客观合理的均方误差表达式。定义s(p)经过噪声近似处理后的信号为

s′(p)=A′exp{j[φ₀′+2πΔfpT_s+e_p]}

定义一个比例参数G来衡量噪声近似处理对原始信号的改变情况：

$G=\frac{s\left(p\right)}{s^{\′}\left(p\right)}=\frac{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s\left]\right\}+v^{\′\′}\left(p\right)}{A^{\′}\exp \left\{j\right[{\mathrm{\φ}}_0^{\′}+2\mathrm{\π\Δ}{\mathrm{fpT}}_s+e_p\left]\right\}}$

对上式求均方误差

$\mathrm{varG}=E\left\{{|\frac{s\left(p\right)}{s^{\′}\left(p\right)}-E[\frac{s\left(p\right)}{s^{\′}\left(p\right)}\left]\right|}^2\right\}=\frac{1}{\mathrm{MSNR}}$

该值varG是原始信号s(p)相对于近似信号s′(p)在均方意义上的改变量。因此修正的均方误差表达式为

$\begin{array}{c}E\left[{(f-\hat{f})}^2\right]=(1+\mathrm{varG})\frac{6}{{\left({2\mathrm{\πt}}_s\right)}^{2}N(N^2-M^2)\mathrm{SNR}}\\ =\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\end{array}$

通过求取上式的最小值，可以得到最优的叠加长度M_opt

$M_{\mathrm{opt}}=\arg \underset{M}{\min }\left\{\frac{6(1+\mathrm{MSNR})}{{\left({2\mathrm{\πt}}_s\right)}^2\mathrm{MN}(N^2-M^2)\mathrm{SNR}}\right\}$

对上式求导并令导数值为零，可得最优叠加长度M_opt满足下式

${2\mathrm{SNRM}}_{\mathrm{opt}}^3+{3M}_{\mathrm{opt}}^2-N^2=0$

从上式中可以求出最优的叠加长度M_opt。

（2）与传统的一些算法的性能比较

图5给出了本发明与传统的一些估计方法的性能比较图。如图5所示，将本发明与傅里叶法、傅里叶插值法、迭代傅里叶法、傅里叶与子空间结合法以及克拉美罗下界进行对比。克拉美罗下界是用来描述频率估计算法理论上可以达到的最小均方误差。

从图5可以看到，本发明提出的算法性能可以达到理论上的最优估计性能，并且明显好于傅里叶法、傅里叶插值法和傅里叶与子空间结合法法。与迭代傅里叶算法相比，本发明的性能在信噪比-14dB到-8dB的区间上稍微差一些。但是必须指出的是，迭代傅里叶算法是一种迭代算法，通过多次重复操作才达到这样的估计性能，而本发明没有经过迭代运算，因此本发明的复杂度比迭代傅里叶算法小很多，更有利于实际的应用。

（3）综合复杂度和估计精度的分析

从图4中可以看到，本发明的性能远远好于除迭代傅里叶算法之外的其他算法，因此在这一节中从复杂度和估计精度两方面合理的比较这两种算法。

迭代傅里叶法是一种迭代的DFT插值算法，设其迭代次数为Q，由于其最少需要迭代两次，所以Q≥2。用快速傅里叶变换（FFT）来实现DFT是一种快速且复杂度低的操作，用FFT来估计频率的复杂度为

在迭代傅里叶算法在迭代过程中还需要2QN次乘法和2QN次加法，因此迭代傅里叶算法的复杂度为

本发明在DFT估计时复杂度为

在其后的信号处理和精估计阶段的复杂度为

因此本发明的复杂度为

除去和迭代傅里叶算法相同的DFT阶段的复杂度，比较两种额外的复杂度

$\frac{4\mathrm{QN}}{2N+20N/M}=\frac{2Q}{1+10/M}\≥\frac{4}{1+10/M}$

由于M值一般远大于10，因此上式近似大于4。于是我们可以看到迭代傅里叶算法的额外复杂度近似为本发明额外复杂度的4倍。

为了比较在复杂度相似情况下两种算法的性能，重新设置了一下本发明的参数值。把本发明在DFT阶段的长度N变为1.5N，其他条件不变。令N=1024，M=128，Q=2，比较一下两种算法总体的复杂度为

$\frac{\frac{9N}{4}{\log }_2\left(\frac{3N}{2}\right)+2N+20N/M}{\frac{3N}{2}{\log }_2\left(N\right)+4\mathrm{QN}}\≈1.1$

可以看到此时两种算法的复杂度十分接近，图6给出了复杂度相同的时候两种算法的性能比较。从图6中可以看到本发明完全胜过迭代傅里叶算法，在信噪比-14dB到-12dB时，本分明还能获得非常理想的性能，而迭代傅里叶算法已经完全失效了。本发明能够不仅能保证很高的估计精度，另外还能明显的降低频率估计算法的信噪比门限值。

Claims (2)

1.一种基于噪声近似处理的双阶段频率估计方法，其特征是： 

步骤一：频率粗估计 

（1）信号的接收 

通信系统中的接收信号一般用含有噪声的复正弦信号进行表示，经过接收机终端采样后，实际处理的信号为： 

x(n)=Aexp[j(φ₀+2πfnt_s)]+v(n)   n=1,2,...,N 

其中，A是信号的幅度，φ₀是信号的初始相位，f是信号的频率，t_s是采样周期，N是信号长度，v(n)是服从高斯分布的加性噪声、其均值为0方差为σ²，频率粗估计阶段的任务是完成对接收信号x(n)的频率粗估计； 

（2）DFT频率估计 

在接收到信号x(n)后，对时域的接收信号x(n)进行傅里叶变换，得相应的频域信号为 

信号的频谱信息为P(k)=|X(k)|²，根据信号频谱的特性，在信号真实频率处所对应的频谱值最大，通过寻找频谱值最大的点来估计信号的频率，设第m个点的频谱值最大，根据对应关系信号的频率估计值为 

（3）本地参考信号的构建 

利用DFT对接收信号x(n)进行频率估计，得到了频率粗估计值

构造一个本地参考信号，所述本地参考信号以估计值

作为参考信号的频率，构建的本地参考信号为： 

本地参考信号包含的是频率粗估计的信息； 

（4）频差信号的构建 

利用接收信号x(n)与本地参考信号y(n)，构造一个反映粗估计频差的频差信号，通过将接受信号x(n)与构建的本地参考信号y(n)共轭相乘，得到频差信号为： 

z(n)=x(n)×y(n)^*

=Aexp[j(φ₀+2πΔfnt_s)]+v′(n) 

n=1,2,...,N； 

其中，

为粗估计的频率偏差、该值是由于DFT估计精度不高而导致的，v′(n)为加 性噪声项，经过上述变换，得到一个以Δf为频率的频差信号z(n)，频差信号z(n)的频率信息即为在粗估计阶段中产生的频率估计偏差，频差Δf的值将在频率精估计阶段进行估计； 

步骤二：频率精估计 

（1）频差信号的叠加处理 

将z(n)按照每M个采样点叠加在一起，得到 

其中，

为s(p)的幅度，P为s(p)的长度，其值为向下取整最接近

的整数， 

为加性噪声项， 

经过叠加处理后得到的信号s(p)是一个信噪比较高的信号； 

（2）信号噪声近似处理 

引入噪声近似处理操作，将加性噪声转换为相位噪声，然后直接提取信号的相位信息，进而实现对频差Δf的快速高效估计； 

当信号的信噪比SNR远大于1时，根据噪声近似理论，此时信号的加性噪声近似等价于相应的相位噪声，于是有 

s(p)≈A′exp{j[φ₀′+2πΔfpT_s+e_p]}   p=1,2,...,P 

其中，φ₀′=φ₀+πΔf(M-1)t_s，T_s=Mt_s，e_p为利用噪声近似理论得到的相位噪声项； 

（3）最小二乘估计 

经过噪声近似处理后，s(p)的近似信号中不再含有加性噪声项，直接提取相位，得到的相位信息为θ_p=φ₀′+2πΔfpT_s+e_p，定义向量Θ=[θ₀,θ₁,...,θ_P-1]^T，H₀=[1,1,...,1]^T，H₁=[0,2πT_s,...,2π(P-1)T_s]^T，H=[H₀,H₁]，X=[φ₀′,Δf]^T，E=[e₀,e₁,...,e_P-1]^T，于是得到相位的向量表达式为 

Θ=HX+E 

根据最小二乘估计算法，令估计误差最小，得到的X估计： 

对上式求导，并令导数为零，得到含有频差Δf信息的项X估计值为： 

从

中提取频差估计值

于是得到最终的信号频率估计为

。

2.根据权利要求1所述的基于噪声近似处理的双阶段频率估计方法，其特征是所述的频差信号的叠加处理的最优的叠加长度的选择方法为： 

定义s(p)经过噪声近似处理后的信号为 

s′(p)=A′exp{j[φ₀′+2πΔfpT_s+e_p]} 

定义一个比例参数G来衡量噪声近似处理对原始信号的改变情况： 

对上式求均方误差 

该值varG是原始信号s(p)相对于近似信号s′(p)在均方意义上的改变量，因此修正的均方误差表达式为 

通过求取上式的最小值，得到最优的叠加长度M_opt

对上式求导并令导数值为零，得最优叠加长度M_opt满足下式 

从上式中求出最优的叠加长度M_opt。 

Images