CN104502699A - 基于数据延拓和Hilbert变换的频率估计方法 - Google Patents

Abstract

本发明涉及信号处理领域，特别是信号的频率估计。本发明的适用对象为任意同频信号，包括以下步骤：首先，对信号进行数据延拓，获得数据延拓信号；其次，对数据延拓信号进行Hilbert变换，获得解析信号；然后，计算解析信号的自相关，并进行时域平均处理，获得时域平均信号；最后，根据时域平均信号生成参考复正弦信号，利用时域平均信号和参考复正弦信号构造误差函数，求解误差函数最小值，其对应的频率即为信号的频率估计值。本发明涉及的正弦信号的频率估计方法实现简单、抗噪性能好、精度高，消除了信号信息融合对相位差的依赖性，提高了信号的频率估计精度。

Description

基于数据延拓和Hilbert变换的频率估计方法

技术领域

本发明涉及信号处理领域，特别是信号的频率估计。

背景技术

信号的频率估计是从含有噪声的采样信号中检测出信号的频率值，是现代信号处理中最具有实用价值的技术之一。信号的频率估计主要是通过时域变换或者频谱分析获得信号的频率估计值。随着科学技术的发展，信号的频率估计广泛应用于雷达、通信、声纳、设备故障诊断和仪器仪表等领域，具有重要的研究意义和应用价值。

近年来涌现了多种信号的频率估计方法，主要有比值法、相位差法、两步自相关法、时域平均法和频谱融合法等。

(1)比值法(参考文献[1]：D.Rife，R.Boorstyn.Single tone parameter estimation fromdiscrete-time observations[J].IEEE Transactions on Information Theory，1974，20(5)：591-598.)，该方法利用主瓣内最大谱线和次大谱线幅值谱的比值求得频率偏移量，进而获得信号频率估计值。方法简单，计算量小，但是当被估计频率接近量化频率时，频率估计性能较差。

(2)相位差法(参考文献[2]：丁康，钟舜聪，朱小勇.通用的离散频谱相位差校正方法[J].电子学报，2003，31(1)：142-145.)，该方法利用信号两次FFT频谱最大谱线相位差包含频率偏移量信息的特点进行频谱校正的，相位差法的普适性强，适用于任何对称窗函数，算法简单，实时性强，抗噪性能有待提高。

(3)两步自相关法(参考文献[3]：K.Lui，H.C.So.Two-stage autocorrelation approach foraccurate single sinusoidal frequency estimation[J].Signal Processing，2008，88(7)：1852-1857.)，通过信号的两步自相关提高信噪比，进而改善信号频率估计精度，但是在信噪比比较低时边界效应明显，并且信号频率接近0.5π频率估计性能较差。

(4)时域平均法(参考文献[4]：纪铁军，任丽军，赵爱明.不用同步信号的相干平均弱信号检测法[J].哈尔滨理工大学学报，2002，7(3)：35-37.)，该方法利用信号自相关获得信号的自相关序列，消除多段信号初相位不同对多段信号信息融合的影响，再通过时域变换或频谱分析获得信号的频率估计值。该方法原理简单，易于实现，但是由于其为一种有偏估计，所以信噪比较高时频率估计精度较差。

(5)频谱融合法(参考文献[5]：肖玮.频率估计的多段正弦信号频谱融合法及LFMCW雷达测距系统验证[D].博士学位论文，后勤工程学院，2012.)，该方法利用多段信号间信号频率、相位存在关联的特点通过相位积累近似达到同等长度相位连续信号的频谱分析效果，通过分段信号功率谱和相位累积频谱的频谱相关提高信噪比来改善信号频谱分析精度。该方法适用范围广，但方法复杂，计算量大，并对多段信号的初相位依赖性较强。

综上所述，信号的频率估计方法具有重要研究意义和应用价值，但现有方法存在诸多问题，需要提出一种频率估计精度高、适用范围广、计算量小和抗噪性能好的频率估计方法。

发明内容

本发明的目的是提出一种信号的频率估计方法，提高信号的频率估计精度，改善实时性和普适性。

本发明信号频率估计方法说明如下：

以正弦信号为例进行说明。

设定正弦信号为M段单一频率实正弦信号，

x(m，n)＝s(m，n)+z(m，n)，n＝1，2…N(m)   (1)

式中，ω为正弦信号的圆周频率，A(m)、N(m)分别为第m段信号的幅值、初相位和信号长度，z(m，n)是均值为零、方差为σ²的高斯白噪声。

为有效抑制信号Hilbert变换的端点效应，对第m段信号数据延拓D(m)个采样点，获得数据延拓信号{x_e(m，n)}。第m段数据延拓信号为

{x_e(m，n)}＝{x(m，1)，x(m，2)，…x(m，N)…x(m，N_e(m))}   (2)

式中，n＝1，2…N_e(m)，N_e(m)＝N(m)+D(m)，第m段信号第N(m)+d点的数值为

$x(m,N\left(m\right)+d)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}[2\cos \left(b{\mathrm{\ω}}_1\right)x(m,N\left(m\right)+d-b)-x(m,N\left(m\right)+d-2b)]---\left(3\right)$

ω₁表示信号的粗略估计值，d∈[0，D(m)]，数据延拓长度D(m)满足

$|\frac{{\mathrm{\ω}}_1{N}_e\left(m\right)}{2}-\mathrm{round}(\frac{{\mathrm{\ω}}_1{N}_e\left(m\right)}{2}/\mathrm{\π})\mathrm{\π}|/\mathrm{\π}\≤{\mathrm{\δ}}_A---\left(4\right)$

对数据延拓信号{x_e(m，n)}进行Hilbert变换，得到解析信号{x_h(m，n)}。

计算解析信号{x_h(m，n)}的自相关，得到自相关信号{r(m，k)}为

$r(m,k)=\underset{n=1}{\overset{N_e\left(m\right)-k}{\mathrm{\Σ}}}{x}_h(m,n+k)x_h^{*}(m,n),k=[1,K]---\left(5\right)$

对自相关信号{r(m，k)}进行时域平均，获得时域平均信号λ为

$\mathrm{\λ}\left(k\right)=\frac{1}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}(N_e\left(m\right)-k)}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}r(m,k)={\mathrm{\λ}}_R\left(k\right)+j{\mathrm{\λ}}_I\left(k\right)---\left(6\right)$

根据时域平均信号生成参考复正弦信号为

$h\left(k\right)=e^{\mathrm{jk}\widehat{\mathrm{\ω}}}=h_R\left(k\right)+j{h}_I\left(k\right)---\left(7\right)$

式中，h_R(k)为参考复正弦信号h第k个元素h(k)的实部，h_I(k)为参考复正弦信号第k个元素h(k)的虚部。

利用时域平均信号λ和参考复正弦信号h构造误差函数为

$\begin{array}{c}J(p,\widehat{\mathrm{\ω}})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[\mathrm{p\λ}\left(k\right)-h\left(k\right)]{[\mathrm{p\λ}\left(k\right)-h\left(k\right)]}^{*}\\ =\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\{{[p{\mathrm{\λ}}_R\left(k\right)-h_R\left(k\right)]}^2+{[p{\mathrm{\λ}}_I\left(k\right)-h_I\left(k\right)]}^2\}\end{array}---\left(8\right)$

由式(8)可知误差函数当且仅当h(k)＝pλ(k)时等于零，此时所以误差函数的最小值对应的频率即为频率估计值。为获得信号频率估计值的解析解，在信号频率的粗略估计值ω₁对h_R(k)和h_I(k)进行泰勒级数展开，并保留一次项为：

$\left\{\begin{array}{c}{h}_R\left(k\right)=\cos \left(\mathrm{k\ω}\right)=\cos \left(k{\mathrm{\ω}}_1\right)-k\sin \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)\\ h_I\left(k\right)=\sin \left(\mathrm{k\ω}\right)=\sin \left(k{\mathrm{\ω}}_1\right)+k\cos \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)\end{array}\right.---\left(9\right)$

将式(9)代入式(8)得

$J(p,\widehat{\mathrm{\ω}})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\{{[p{\mathrm{\λ}}_R\left(k\right)-\cos \left(k{\mathrm{\ω}}_1\right)+k\sin \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)]}^2+{[p{\mathrm{\λ}}_I\left(k\right)-\sin \left(k{\mathrm{\ω}}_1\right)-k\cos \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)]}^2\}---\left(10\right)$

当取得最小值时，有

$\left\{\begin{array}{c}\frac{\∂J(p,\widehat{\mathrm{\ω}})}{\∂p}=2\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\left\{\begin{array}{c}{\mathrm{\λ}}_R\left(k\right)[p{\mathrm{\λ}}_R\left(k\right)-\cos \left(k{\mathrm{\ω}}_1\right)+k\sin \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)]\\ +{\mathrm{\λ}}_I\left(k\right)[p{\mathrm{\λ}}_I\left(k\right)-\sin \left(k{\mathrm{\ω}}_1\right)-k\cos \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)]\end{array}\right\}=0\\ \frac{\∂J(p,\widehat{\mathrm{\ω}})}{\∂\widehat{\mathrm{\ω}}}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\left\{\begin{array}{c}k\sin \left(k{\mathrm{\ω}}_1\right)[p{\mathrm{\λ}}_R\left(k\right)-\cos \left(k{\mathrm{\ω}}_1\right)+k\sin \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)]\\ -k\cos \left(k{\mathrm{\ω}}_1\right)[p{\mathrm{\λ}}_I\left(k\right)-\sin \left(k{\mathrm{\ω}}_1\right)-k\cos \left(k{\mathrm{\ω}}_1\right)(\widehat{\mathrm{\ω}}-{\mathrm{\ω}}_1)]\end{array}\right\}=0\end{array}\right.---\left(11\right)$

根据式(11)得到信号的频率估计值为

$\widehat{\mathrm{\ω}}={\mathrm{\ω}}_1+\frac{c_1{a}_{21}-c_2{a}_{11}}{a_{12}{a}_{21}-a_{22}{a}_{11}}---\left(12\right)$

式中， $a_{11}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[{\mathrm{\λ}}_R^2\left(k\right)+{\mathrm{\λ}}_I^2\left(k\right)];a_{12}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[k\sin \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_R\left(k\right)-k\cos \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_I\left(k\right)];$

$a_{21}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[k\sin \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_R\left(k\right)-k\cos \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_I\left(k\right)];a_{22}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[{\left[k\sin \left(k{\mathrm{\ω}}_1\right)\right]}^2+{\left[k\cos \left(k{\mathrm{\ω}}_1\right)\right]}^2];$

$c_1=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[{\mathrm{\λ}}_R\left(k\right)\cos \left(k{\mathrm{\ω}}_1\right)+{\mathrm{\λ}}_I\left(k\right)\sin \left(k{\mathrm{\ω}}_1\right)];$ c₂＝0。

附图说明

下面根据附图和具体实施方式对本发明进一步阐述。以M＝4，信号类型为正弦信号为例进行说明。

图1为基于数据延拓和Hilbert变换的频率估计方法的基本思想。

图中：1表示M段正弦信号；2表示数据延拓信号；3表示解析信号；4表示自相关信号；5表示时域平均信号；6表示参考复正弦信号；7表示误差函数；8表示频率估计值；9表示数据延拓；10表示Hilbert变换；11表示自相关；12表示初相位匹配，13表示频率匹配；14表示计算误差函数最小值。

图2为正弦信号时域示意图。

图中：15表示第一段信号；16表示第二段信号；17表示第三段信号；18表示第四段信号。

图3为时域平均信号时域示意图。

图中：19表示时域平均信号实部；20表示时域平均信号实部理论值；21表示时域平均信号虚部；22表示时域平均信号虚部理论值。

具体实施方式

本发明的具体实施方式如下：

第一步：利用M段正弦信号获得信号的粗略频率估计值ω₁为

${\mathrm{\ω}}_1=\frac{2\mathrm{\π}}{N}(k_0+\frac{a|X(m_{\max },k_0+a)}{\left|X(m_{\max },k_0)\right|+\left|X(m_{\max },k_0+a)\right|})---\left(13\right)$

其中：k₀为对应FFT最大谱线的标号，a＝±1表示修正方向，当|X(m_max，k₀+1)|＞|X(m_max，k₀-1)|时，a＝1；当|X(m_max，k₀+1)|≤|X(m_max，k₀-1)|时，a＝-1，m_max表示第m_max段信号是M段正弦信号中最长的一段信号。

第二步：对第m段信号数据延拓D(m)个采样点，获得数据延拓信号{x_e(m，n)}。第m段数据延拓信号为

{x_e(m，n)}＝{x(m，1)，x(m，2)，…x(m，N)…x(m，N_e(m))}   (14)

第三步：对数据延拓信号{x_e(m，n)}进行Hilbert变换，获得其解析信号{x_h(m，n)}。

第四步：计算解析信号{x_h(m，n)}的自相关，并进行时域平均处理，获得时域平均信号λ为

$\mathrm{\λ}\left(k\right)=\frac{1}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}(N_e\left(m\right)-k)}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N_e\left(m\right)-k}{\mathrm{\Σ}}}{x}_h(m,n+k)x_h^{*}(m,n)---\left(15\right)$

第五步：根据时域平均信号λ生成参考复正弦信号h，利用时域平均信号λ和参考复正弦信号h构造误差函数求解误差函数最小值，其对应的频率即为信号的频率估计值。

$\widehat{\mathrm{\ω}}={\mathrm{\ω}}_1+\frac{c_1{a}_{21}-c_2{a}_{11}}{a_{12}{a}_{21}-a_{22}{a}_{11}}---\left(16\right)$

Claims (1)

1.基于数据延拓和Hilbert变换的频率估计方法，其特征在于：适用对象为任意同频信号。

该方法包括以下步骤：

第一步：利用M段同频信号获得信号的粗略频率估计值ω₁为

${\mathrm{\ω}}_1=\frac{2\mathrm{\π}}{N}(k_0+\frac{a\left|X(m_{\max },k_0+a)\right|}{\left|X(m_{\max },k_0)\right|+\left|X(m_{\max },k_0+a)\right|})\·\·\·\left(1\right)$

其中：k₀为对应FFT最大谱线的标号，a＝±1表示修正方向，当|X(m_max，k₀+1)|＞|X(m_max，k₀-1)|时，a＝1；当|X(m_max，k₀+1)|≤|X(m_max，k₀-1)|时，a＝-1，m_max表示第m_max段信号是M段信号中最长的一段信号。

第二步：对第m段信号数据延拓D(m)个采样点，获得数据延拓信号{x_e(m，n)}。第m段数据延拓信号为

{x_e(m，n)}＝{x(m，1)，x(m，2)，...x(m，N)...x(m，N_e(m))}  (2)

式中，n＝1，2...N_e(m)，N_e(m)＝N(m)+D(m)，第m段信号第N(m)+d点的数值为

$x(m,N\left(m\right)+d)=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\Σ}}}[2\cos \left({\mathrm{b\ω}}_1\right)x(m,N\left(m\right)+d-b)-x(m,N\left(m\right)+d-2b)]\·\·\·\left(3\right)$

式中，d∈[0，D(m)]，数据延拓长度D(m)满足

$|\frac{{\mathrm{\ω}}_1{N}_e\left(m\right)}{2}-\mathrm{round}(\frac{{\mathrm{\ω}}_1{N}_e\left(m\right)}{2}/\mathrm{\π})\mathrm{\π}|/\≤{\mathrm{\δ}}_A\·\·\·\left(4\right)$

第三步：对数据延拓信号{x_e(m，n)}进行Hilbert变换，获得其解析信号{x_h(m，n)}。

第四步：计算解析信号{x_h(m，n)}的自相关并进行时域平均处理，获得时域平均信号λ为

$\mathrm{\λ}\left(k\right)=\frac{1}{\underset{m=1}{\overset{m}{\mathrm{\Σ}}}(N_e\left(m\right)-k)}\underset{m=1}{\overset{m}{\mathrm{\Σ}}}\underset{n=1}{\overset{N_e\left(m\right)-k}{\mathrm{\Σ}}}{x}_h(m,n+k)x_h^{*}(m,n)\·\·\·\left(5\right)$

第五步：根据时域平均信号λ生成参考复正弦信号h，利用时域平均信号λ和参考复正弦信号h构造误差函数求解误差函数最小值，其对应的频率即为信号的频率估计值。

$\widehat{\mathrm{\ω}}={\mathrm{\ω}}_1+\frac{c_1{a}_{21}-c_2{a}_{11}}{a_{12}{a}_{21}-a_{22}{a}_{11}}\·\·\·\left(6\right)$

式中， $a_{11}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[{\mathrm{\λ}}_R^2\left(k\right)+{\mathrm{\λ}}_1^2\left(k\right)];$ $a_{12}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[k\sin \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_R\left(k\right)-k\cos \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_1\left(k\right)];$

$a_{21}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[k\sin \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_R\left(k\right)-k\cos \left(k{\mathrm{\ω}}_1\right){\mathrm{\λ}}_1\left(k\right)];$ $a_{22}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[{\left[k\sin \left({\mathrm{k\ω}}_1\right)\right]}^2+{\left[k\cos \left({\mathrm{k\ω}}_1\right)\right]}^2];$

$c_1=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}[{\mathrm{\λ}}_R\left(k\right)\cos \left({\mathrm{k\ω}}_1\right)+{\mathrm{\λ}}_I\left(k\right)\sin \left(k{\mathrm{\ω}}_1\right)];$ c₂＝0。