CN112162153A - 基于相位直线拟合的正弦波信号频率估计方法 - Google Patents

Abstract

本发明公开了一种基于相位直线拟合的正弦波信号频率估计方法。该方法步骤为：对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值f_e；利用f_e对输入的复解析数字信号进行下变频处理；对下变频后的数据进行低通滤波处理；数据量若大于设定阈值则进行数据抽取处理；对所得数据进行解相位角运算，并对所得相位角数据进行预处理；进行直线拟合运算，根据结果计算频率估计误差Δf_e，正弦波信号的频率估计值调整为f_e＝f_e+Δf_e；若Δf_e大于阈值，用Δf_e对输入复解析数字信号进行再次下变频处理，提取信号的相位信息重复前述步骤；否则输出正弦波频率的最终估计值f_e。本发明在低信噪比条件下具有更优的估计性能。

Description

基于相位直线拟合的正弦波信号频率估计方法

技术领域

本发明涉及数字信号处理技术领域，特别是一种基于相位直线拟合的正弦波信号频率估计方法。

背景技术

在很多工程应用领域，如雷达/声呐信号处理、通信信号处理、无源定位、计量与测试等领域，都需要对正弦波信号进行高精度频率估计。正弦波信号频率估计的最常用方法是FFT算法，但FFT算法的频率估计性能除受检测信噪比制约之外，还受制于其数字分辨力f_s/M，f_s、M分别为采样频率和FFT运算的点数。当f_s一定时，数字分辨力每提高一倍，数据存贮量和计算量就应相应地增加一倍，因此FFT算法通常不适用于需要对正弦波信号进行高精度频率估计的场合。

为解决此问题，人们又提出了Rife算法及其各种改进型算法(下方中统称为M-Rife算法)。M-Rife算法以FFT算法为基础，在无需提高FFT算法数字分辨力的情况下，可明显改善频率估计性能。但M-Rife算法也是有局限性的。其局限性表现在以下两个方面：1)受噪声影响，在确定最大谱线两侧的次大谱线时，有可能出现误判；2)在较低输入信噪比条件下，真实频率并不一定位于与最大谱线相邻的两根谱线之间的某个位置。当出现上述两种情况时，M-Rife算法反而会恶化频率估计性能。

正弦波频率估计的另一种常用方法是间隔多个采样点提取两个采样数据间的相位差，然后除以这两个数据采样点间的时间间隔以获得频率的粗估计值，再通过对多个粗估计值进行平均运算以提高频率估计精度。该种频率估计方法的估计性能通常劣于FFT法和M-Rife算法，只适用于较高信噪比的条件下。

发明内容

本发明的目的在于提供一种基于相位直线拟合的正弦波信号频率估计方法，从而在较低检测信噪比条件下实现更优的频率估计性能。

实现本发明目的的技术解决方案为：一种基于相位直线拟合的正弦波信号频率估计方法，包括以下步骤：

步骤1，对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值，记作f_e；

步骤2，利用正弦波信号频率的初估计值f_e，对输入的复解析数字信号进行下变频处理；

步骤3，对下变频后的数据进行低通滤波处理；

步骤4，判断经低通滤波处理后的数据量是否大于设定阈值：若大于阈值则对经低通滤波处理后的数据进行数据抽取处理，然后进入步骤5；否则，直接进入步骤5；

步骤5，对所得数据进行解相位角运算；

步骤6，对所得相位角数据进行预处理；

步骤7，对经预处理后的相位角数据进行直线拟合运算，再根据直线似合的结果计算频率估计误差Δf_e，正弦波信号的频率估计值调整为f_e＝f_e+Δf_e；

步骤8，判断频率估计误差Δf_e是否大于给定阈值：若Δf_e大于给定阈值，用Δf_e对经步骤4处理后的输入复解析数字信号进行再次下变频处理，然后转到步骤5；否则，停止正弦波频率迭代估计过程，并输出正弦波频率的最终估计值f_e。

进一步地，步骤1所述对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值，记作f_e，具体如下：

记含噪复解析正弦波信号s₀(t)的采样频率、采样点数分别为f_s、N；对s₀(t)的采样信号序列s₀(n)作M点FFT运算得s₀(t)的幅度谱|S(k)|，n＝1,2,3,…,N-1，M≥N，k＝1,2,3,…,M-1；再计算幅度谱中最强谱线所对应的频率值记作f_e。

进一步地，步骤2所述利用正弦波信号频率的初估计值f_e，对输入的复解析数字信号进行下变频处理，具体如下：

利用正弦波信号频率的初估计值f_e对采样信号序列s₀(n)进行下变频处理，即

得下变频后的采样信号序列

n＝1,2,3,…,N-1。

进一步地，步骤3所述对下变频后的数据进行低通滤波处理，具体如下：

对下变频后的采样信号序列

进行低通滤波处理，得下变频及低通滤波后的采样数据序列

n＝0,1,2,…,N-1；低通滤波处理时，低通滤波器的带宽根据步骤1中f_e的最大估计误差及采样频率f_s共同决定，即低通滤波器的带宽大于f_s的0.04倍，并大于f_e的最大估计误差。

进一步地，步骤4所述数据抽取处理，具体如下：

对下变频及低通滤波后的采样数据序列

进行数据抽取运算，得到新的采样数据序列

n'＝0,1,2,…,Q-1，Q为进行数据抽取后所得采样数据序列的长度；

假设数据抽取因子为η，则进行数据抽取后所得采样数据序列的采样频率f_s'＝f_s/η，f_s'大于低通滤波器带宽的2倍；当不进行数据抽取操作时，η＝1。

进一步地，步骤5所述对所得数据进行解相位角运算，具体如下：

计算数据序列

的相位信息序列θ(n')；记数据序列

其中，s_I(n')、s_Q(n')分别为

的实部和虚部，则

-π≤θ(n')≤π，n'＝0,1,2,…,Q-1。

进一步地，步骤6所述对所得相位角数据进行预处理，具体如下：

首先，判断所提取的相位数据序列是否发生无周期性的相位折叠现象，当出现此种情况时，令

n'＝0,1,2,…,Q-1，并重新提取相位数据序列；

其次，判断是否发生了周期性的相位折叠现象：当发生周期性的相位折叠现象时，取最大没发生相位折叠的相位数据段，用于后续的相位直线拟合运算；

最后，从最大没发生相位折叠的相位数据段中剔除不符合要求的相位估计值，具体方法是：先计算Θ(n')＝θ(n'+1)-θ(n')，n'＝0,1,2,…,Q-1，再求Θ(n')的均值和标准差，当Θ(n')与其均值之差大于其标准差的3倍时，用Θ(n'-1)或Θ(n'+1)代替Θ(n')。

进一步地，步骤7所述对经预处理后的相位角数据进行直线拟合运算，具体如下：

假设所得拟合的直线为

则Δf_e即为频率估计误差；

为拟合后的相位信息序列，θ₀为拟合出的常数。

进一步地，步骤8所述判断频率估计误差Δf_e是否大于给定阈值：若Δf_e大于给定阈值，用Δf_e对经步骤4处理后的输入复解析数字信号进行再次下变频处理，然后转到步骤5；否则，停止正弦波频率迭代估计过程，并输出正弦波频率的最终估计值f_e，具体如下：

当Δf_e大于给定的阈值时，令f_e＝f_e+Δf_e，并利用Δf_e对数据序列

进行下变频处理，即

其中，n'＝0,1,2,…,Q-1，f_s'为进行数据抽取后的采样频率；当Δf_e小于给定的门限值时，输出频率估计值f_e。

本发明与现有技术相比，其显著优点为：(1)将FFT算法与相位直线拟合相结合，克服了M-Rife算法的局限性；(2)基于相位直线拟合的正弦波频率高精度迭代估计，在较低检测信噪比条件下具有更优的频率估计性能。

附图说明

图1为本发明基于相位直线拟合的正弦波信号频率估计方法的具体流程图。

图2为本发明实施例中第一次仿真步骤5所提取的相位数据序列示意图。

图3为本发明实施例中第二次仿真步骤5所提取的相位数据序列示意图。

图4为本发明实施例中第三次仿真步骤5所提取的相位数据序列示意图。

图5为本发明实施例中第四次仿真步骤5所提取的相位数据序列示意图。

图6为本发明实施例与FFT算法、M-Rife算法在几种输入信噪比情况下的估计性能对比图。

具体实施方式

本发明一种基于相位直线拟合的正弦波信号频率估计方法，包括以下步骤：

步骤1，对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值，记作f_e；

步骤2，利用正弦波信号频率的初估计值f_e，对输入的复解析数字信号进行下变频处理；

步骤3，对下变频后的数据进行低通滤波处理；

步骤4，判断经低通滤波处理后的数据量是否大于设定阈值：若大于阈值则对经低通滤波处理后的数据进行数据抽取处理，然后进入步骤5；否则，直接进入步骤5；

步骤5，对所得数据进行解相位角运算；

步骤6，对所得相位角数据进行预处理；

步骤7，对经预处理后的相位角数据进行直线拟合运算，再根据直线似合的结果计算频率估计误差Δf_e，正弦波信号的频率估计值调整为f_e＝f_e+Δf_e；

步骤8，判断频率估计误差Δf_e是否大于给定阈值：若Δf_e大于给定阈值，用Δf_e对经步骤4处理后的输入复解析数字信号进行再次下变频处理，然后转到步骤5；否则，停止正弦波频率迭代估计过程，并输出正弦波频率的最终估计值f_e。

进一步地，步骤1所述对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值，记作f_e，具体如下：

记含噪复解析正弦波信号s₀(t)的采样频率、采样点数分别为f_s、N；对s₀(t)的采样信号序列s₀(n)作M点FFT运算得s₀(t)的幅度谱|S(k)|，n＝1,2,3,…,N-1，M≥N，k＝1,2,3,…,M-1；再计算幅度谱中最强谱线所对应的频率值记作f_e。

进一步地，步骤2所述利用正弦波信号频率的初估计值f_e，对输入的复解析数字信号进行下变频处理，具体如下：

利用正弦波信号频率的初估计值f_e对采样信号序列s₀(n)进行下变频处理，即

得下变频后的采样信号序列

n＝1,2,3,…,N-1。

进一步地，步骤3所述对下变频后的数据进行低通滤波处理，具体如下：

对下变频后的采样信号序列

进行低通滤波处理，得下变频及低通滤波后的采样数据序列

n＝0,1,2,…,N-1；低通滤波处理时，低通滤波器的带宽根据步骤1中f_e的最大估计误差及采样频率f_s共同决定，即低通滤波器的带宽大于f_s的0.04倍，并大于f_e的最大估计误差。

进一步地，步骤4所述数据抽取处理，具体如下：

对下变频及低通滤波后的采样数据序列

进行数据抽取运算，得到新的采样数据序列

n'＝0,1,2,…,Q-1，Q为进行数据抽取后所得采样数据序列的长度；

假设数据抽取因子为η，则进行数据抽取后所得采样数据序列的采样频率f_s'＝f_s/η，f_s'大于低通滤波器带宽的2倍；当不进行数据抽取操作时，η＝1。

进一步地，步骤5所述对所得数据进行解相位角运算，具体如下：

计算数据序列

的相位信息序列θ(n')；记数据序列

其中，s_I(n')、s_Q(n')分别为

的实部和虚部，则

-π≤θ(n')≤π，n'＝0,1,2,…,Q-1。

进一步地，步骤6所述对所得相位角数据进行预处理，具体如下：

首先，判断所提取的相位数据序列是否发生无周期性的相位折叠现象，当出现此种情况时，令

n'＝0,1,2,…,Q-1，并重新提取相位数据序列；

其次，判断是否发生了周期性的相位折叠现象：当发生周期性的相位折叠现象时，取最大没发生相位折叠的相位数据段，用于后续的相位直线拟合运算；

最后，从最大没发生相位折叠的相位数据段中剔除不符合要求的相位估计值，具体方法是：先计算Θ(n')＝θ(n'+1)-θ(n')，n'＝0,1,2,…,Q-1，再求Θ(n')的均值和标准差，当Θ(n')与其均值之差大于其标准差的3倍时，用Θ(n'-1)或Θ(n'+1)代替Θ(n')。

进一步地，步骤7所述对经预处理后的相位角数据进行直线拟合运算，具体如下：

假设所得拟合的直线为

则Δf_e即为频率估计误差；

为拟合后的相位信息序列，θ₀为拟合出的常数。

进一步地，步骤8所述判断频率估计误差Δf_e是否大于给定阈值：若Δf_e大于给定阈值，用Δf_e对经步骤4处理后的输入复解析数字信号进行再次下变频处理，然后转到步骤5；否则，停止正弦波频率迭代估计过程，并输出正弦波频率的最终估计值f_e，具体如下：

当Δf_e大于给定的阈值时，令f_e＝f_e+Δf_e，并利用Δf_e对数据序列

进行下变频处理，即

其中，n'＝0,1,2,…,Q-1，f_s'为进行数据抽取后的采样频率；当Δf_e小于给定的门限值时，输出频率估计值f_e。

下面结合附图及具体实施例对本发明作进一步详细描述。

实施例

图1示出了本发明采用的单一正弦波频率估计方法的具体处理步骤。本发明中待处理含噪复解析正弦性信号可表示成式(1)所示形式。

s_r(t)＝Aexp(j2πf₀t+θ)+w_r(t),(0≤t≤T) (1)

式(1)中，f₀、A、θ、T和w_r(t)分别为正弦波信号的频率、信号幅度、初相角、信号持续时间和观测噪声，其中，f₀为待估计量。文中假设观测噪声w_r(t)为零均值复高斯白噪声。待处理信号s_r(t)的相位可表示成式(2)所示形式(本发明虽然是对采样数字信号进行处理，但文中仍用连续时间变量t表征信号相位与时间间的对应关系，这并不影响对本发明的描述)。

φ(t)＝2πf₀t+θ+w(t),(0≤t≤T) (2)

式(2)中，w(t)为由w_r(t)所产生的相位噪声。当w_r(t)为均值为零的随机噪声时，w(t)也是均值为零的随机噪声。当w(t)可以忽略不计时，理论上讲，φ(t)是随时间t线性变化的。但考虑到φ(t)只能在(-π,π)范围内取值的约束条件，则在整个观测时间内，φ(t)实际上应是周期线性变化的，除非在整个观测时间内待处理正弦波信号仅在一个周期内变化。当w(t)不能忽略不计时，除非w(t)足够强以至于淹没掉φ(t)，否则φ(t)的这种变化规律仍应能观测到。

当对式(1)中所示信号进行数字采样，并进行FFT运算时，幅度谱中最强谱线所对应的频率即为f₀的粗估计值f_e。利用f_e对输入信号进行下变频及低通滤波处理后，再提取的信号相位信息可表示成式(3)所示形式。

φ_LP(t)＝2πΔft+θ+w_LP(t),(0≤t≤T) (3)

式(3)中，Δf＝f₀-f_e。对输入信号进行下变频及低通滤波处理带来两方面的好处：一是大大降低了相位估计值φ_LP(t)中的噪声信号分量w_LP(t)；另一方面，因频率估计误差Δf通常远小于f₀，所以，所提取的信号相位信息或许不再是周期变化的，至少变化周期显著增大。

图2、图3、图4及图5分别示出了某输入信噪比下四次仿真所提取的相位估计值序列(信号采样频率、采样点数分别为2000Hz、256点；在对正弦波信号的频率进行粗估计时，采用256点FFT运算；数字下变频后所采用的低通滤波器的截止频率约为80Hz)。图2中所提取的相位估计值序列没有出现相位折叠现象，也没出现明显超差的相位估计值，可以直接对该相位估计值序列进行直线拟合以计算频差的估计值Δf。图3中所提取的相位估计值序列既出现了局部相位折叠现象，也出现了明显超差的相位估计值。此种情况下，可先取发生局部相位折叠前的相位数据序列进行首次直线拟合运算。图4中所提取的相位估计值序列出现了周期性相位折叠现象，但没出现明显超差的相位估计值。此种情况下，既可先取持续时间最长的没发生相位折叠的局部相位数据序列进行首次直线拟合运算，也可采取解相位折叠的方法先对相位数据序列进行线性化，再进行直线拟合运算，或者直接根据采样数据长度及相位折叠的周期数选对频率估计误差进行粗估计。图5中所提取的相位数据序列因基本上都位于-π或π附近，所以发生了较频繁的相位折叠现象。当出现此种情况时，应调整经下变频及低通滤波后的信号的初相位，并重新提取相位数据序列。

得到频率估计差的首次估计值Δf_e后，若Δf_e小于或等于某一给定门限值，则输出频率估计值f_e，频率估计过程结束。当Δf_e大于某一给定门限值时，一方面将频率估计值按式(4)进行调整，另一方面利用Δf_e对已下变频并低通滤波处理后的输入采样数据序列再次进行下变频处理，然后，重新提取相位数据序列，并进行相位数据序列的预处理、相位直线拟合、频率误差估计等过程，直至频率估计误差小于或等于某一给定门限值。通常，经过两次频率误差估计值的迭代运算即可结束频率估计过程。

f_e＝f_e+Δf_e (4)

图6示出了本发明一实施例与FFT算法、M-Rife算法在几种输入信噪比情况下的估计性能对比图(信号采样频率、采样点数分别为2000Hz、256点；在对正弦波信号的频率进行粗估计时，M-Rife算法及本发明均采用256点FFT运算；FFT算法采用1024点FFT运算；正弦波信号的频率在100Hz附近随机设置)。由图6可见：在给定的输入信噪比范围内本发明所提供频率估计方法的估计性能均优于M-Rife算法；当输入信噪比低于-5dB时，本发明所提供频率估计方法的估计性能虽劣于高数字分辨力的FFT算法，但仍具有待数据处理量较小的优势。

Claims (9)

1.一种基于相位直线拟合的正弦波信号频率估计方法，其特征在于，包括以下步骤：

步骤1，对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值，记作f_e；

步骤2，利用正弦波信号频率的初估计值f_e，对输入的复解析数字信号进行下变频处理；

步骤3，对下变频后的数据进行低通滤波处理；

步骤4，判断经低通滤波处理后的数据量是否大于设定阈值：若大于阈值则对经低通滤波处理后的数据进行数据抽取处理，然后进入步骤5；否则，直接进入步骤5；

步骤5，对所得数据进行解相位角运算；

步骤6，对所得相位角数据进行预处理；

步骤7，对经预处理后的相位角数据进行直线拟合运算，再根据直线似合的结果计算频率估计误差Δf_e，正弦波信号的频率估计值调整为f_e＝f_e+Δf_e；

步骤8，判断频率估计误差Δf_e是否大于给定阈值：若Δf_e大于给定阈值，用Δf_e对经步骤4处理后的输入复解析数字信号进行再次下变频处理，然后转到步骤5；否则，停止正弦波频率迭代估计过程，并输出正弦波频率的最终估计值f_e。

2.根据权利要求1所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤1所述对输入的含噪复解析数字信号作FFT运算，并将最强谱线所对应的频率值作为正弦波信号频率的初估计值，记作f_e，具体如下：

记含噪复解析正弦波信号s₀(t)的采样频率、采样点数分别为f_s、N；对s₀(t)的采样信号序列s₀(n)作M点FFT运算得s₀(t)的幅度谱|S(k)|，n＝1,2,3,…,N-1，M≥N，k＝1,2,3,…,M-1；再计算幅度谱中最强谱线所对应的频率值记作f_e。

3.根据权利要求2所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤2所述利用正弦波信号频率的初估计值f_e，对输入的复解析数字信号进行下变频处理，具体如下：

利用正弦波信号频率的初估计值f_e对采样信号序列s₀(n)进行下变频处理，即

得下变频后的采样信号序列

n＝1,2,3,…,N-1。

4.根据权利要求3所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤3所述对下变频后的数据进行低通滤波处理，具体如下：

对下变频后的采样信号序列

进行低通滤波处理，得下变频及低通滤波后的采样数据序列

n＝0,1,2,…,N-1；低通滤波处理时，低通滤波器的带宽根据步骤1中f_e的最大估计误差及采样频率f_s共同决定，即低通滤波器的带宽大于f_s的0.04倍，并大于f_e的最大估计误差。

5.根据权利要求4所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤4所述数据抽取处理，具体如下：

对下变频及低通滤波后的采样数据序列

进行数据抽取运算，得到新的采样数据序列

n'＝0,1,2,…,Q-1，Q为进行数据抽取后所得采样数据序列的长度；

假设数据抽取因子为η，则进行数据抽取后所得采样数据序列的采样频率f_s'＝f_s/η，f_s'大于低通滤波器带宽的2倍；当不进行数据抽取操作时，η＝1。

6.根据权利要求5所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤5所述对所得数据进行解相位角运算，具体如下：

计算数据序列

的相位信息序列θ(n')；记数据序列

其中，s_I(n')、s_Q(n')分别为

的实部和虚部，则

-π≤θ(n')≤π，n'＝0,1,2,…,Q-1。

7.根据权利要求6所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤6所述对所得相位角数据进行预处理，具体如下：

首先，判断所提取的相位数据序列是否发生无周期性的相位折叠现象，当出现此种情况时，令

n'＝0,1,2,…,Q-1，并重新提取相位数据序列；

其次，判断是否发生了周期性的相位折叠现象：当发生周期性的相位折叠现象时，取最大没发生相位折叠的相位数据段，用于后续的相位直线拟合运算；

最后，从最大没发生相位折叠的相位数据段中剔除不符合要求的相位估计值，具体方法是：先计算Θ(n')＝θ(n'+1)-θ(n')，n'＝0,1,2,…,Q-1，再求Θ(n')的均值和标准差，当Θ(n')与其均值之差大于其标准差的3倍时，用Θ(n'-1)或Θ(n'+1)代替Θ(n')。

8.根据权利要求7所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤7所述对经预处理后的相位角数据进行直线拟合运算，具体如下：

假设所得拟合的直线为

则Δf_e即为频率估计误差；

为拟合后的相位信息序列，θ₀为拟合出的常数。

9.根据权利要求8所述的基于相位直线拟合的正弦波信号频率估计方法，其特征在于，步骤8所述判断频率估计误差Δf_e是否大于给定阈值：若Δf_e大于给定阈值，用Δf_e对经步骤4处理后的输入复解析数字信号进行再次下变频处理，然后转到步骤5；否则，停止正弦波频率迭代估计过程，并输出正弦波频率的最终估计值f_e，具体如下：

当Δf_e大于给定的阈值时，令f_e＝f_e+Δf_e，并利用Δf_e对数据序列

进行下变频处理，即

其中，n'＝0,1,2,…,Q-1，f_s'为进行数据抽取后的采样频率；当Δf_e小于给定的门限值时，输出频率估计值f_e。

Images