CN113281566A - 一种基于组合复信号相位差的频率估计方法 - Google Patents

Abstract

本发明属于信号处理领域，具体涉及一种基于组合复信号相位差的频率估计方法；所述方法包括通过信号采集器对原始信号采样，获得长度为N+M的数字信号；从该段数字信号0点起取N点构造出第一时间序列，从该段数字信号第M点起取N点构造出第二时间序列；构造出长度为N的对称窗函数，分别对第一时间序列和第二时间序列进行对称加窗处理得到第一加窗函数和第二加窗函数；构造出第三时间序列，对第三时间序列进行离散傅里叶变换得到该序列的频谱；利用频谱的峰值谱线解算得到相位差，通过相位差估计得到原始信号的频率值。本发明只需要一次FFT运算即可对信号频率进行高精度估计，大大减少了算法的运算量，提高了实时性。

Description

一种基于组合复信号相位差的频率估计方法

技术领域

本发明属于信号处理领域，涉及一种用于信号中频率估计的方法，具体涉及一种基于组合复信号相位差的频率估计方法。

背景技术

精确而高效地测量谐波信号的幅值、频率和相位是信号处理中的研究重点。工程中通常采用快速傅里叶变换(Fast Transform,FFT)，并通过离散频谱分析对测量信号进行频谱参数估计。然而，时域信号的截断和离散化处理，不可避免会造成能量泄漏和栅栏效应，最终带来较大的估计误差。针对这一问题，国内外学者提出了离散频谱校正技术，主要包括插值校正法、能量重心校正法、相位差校正法和FFT+FT连续细化校正法。相位差校正法适用于各类窗函数，并且具有将强的抗噪性能。其中，基于时移的相位差校正法应用较为广泛。

目前，基于时移的相位差校正法存在如下缺点：

1、需要进行两次FFT运算，计算量较大，极大地影响了算法的实时性；

2、在随机噪声干扰下存在由于谱线定位错误导致的相位绕卷问题，最终造成较大的估计误差。

发明内容

本发明旨在解决上述技术问题。提出了一种基于组合复信号相位差的频率估计方法。本发明的技术方案如下：

一种基于组合复信号相位差的频率估计方法，包括：

通过信号采集器对原始信号采样，获得长度为N+M的数字信号；

从该段数字信号0点起取N点构造出第一时间序列，再从第M点起取N点构造出第二时间序列；

构造出长度为N的对称窗函数，分别对所述第一时间序列和所述第二时间序列进行对称加窗处理得到第一加窗函数和第二加窗函数；

将所述第一加窗函数和所述第二加窗函数构造出第三时间序列，对所述第三时间序列进行离散傅里叶变换得到该序列的频谱；

利用所述频谱的峰值谱线解算得到第一时间序列与第二时间序列的初始相位差，通过所述相位差估计得到原始信号的频率值。

进一步的，所述利用所述频谱的峰值谱线解算得到相位差包括利用频率幅值搜索得到所述频谱中峰值谱线对应的频点l，将其正负频点±l带入离散傅里叶变换后的第三时间序列中，得到第一参量和第二参量；根据欧拉公式对所述第一参量和所述第二参量化简，得到第一相位差和第二相位差：当正频点l处的幅值小于负频点-l处的幅值时，所述解算得到相位差取第一相位差，反之取第二相位差。

进一步的，将所述相位差估计得到原始信号的频率值表示为：

其中，

表示原始信号的频率估计值；l表示频谱中峰值谱线对应的频点；

表示调整后的归一化频率校正量；f_s为采样频率。

进一步的，所述调整后的归一化频率校正量的获取过程包括根据时移系数和取模后的相位差计算出归一化频率校正量

按照最高谱线幅值|Y(l)|、次高谱线幅值|Y(l±1)|、第三高谱线幅值

计算出第一调整参数和第二调整参数；根据所述第一调整参数和所述第二调整参数对所述归一化频率校正量的符号进行调整，调整后的归一化频率校正量的公式表示为

ε＝sign(ε₁+ε₂)；ε₁表示第一调整参数；ε₂表示第二调整参数；sign(·)为符号函数。

进一步的，所述归一化频率校正量

的计算过程包括基于相位绕卷构建出相位差等式，对所述相位差等式两边取模，对取模后可变相位进行调整，将调整后的可变相位与时移系数的比值作为归一化频率校正量。

本发明的有益效果：

本发明的基于组合复信号的相位差法只需要进行一次FFT运算即可对信号频率进行高精度估计，相较于基于时移的相位差校正法，大大减少了算法的运算量，提高了实时性。并且，本发明能有效降低相位绕卷问题导致的较大频率估计误差，具有良好的抗噪声高性能。

附图说明

图1是本发明实施例中基于组合复信号相位差的频率估计方法流程图；

图2是不同信噪比下现有技术在加汉宁窗时谱线定位错误的概率曲线图；

图3是不同信噪比下本发明中谱线定位错误概率曲线图；

图4是时移长度为L＝96时本发明与现有技术基于时移相位差的频率估计均方根误差(RMSE)曲线图；

图5是时移长度为L＝112时本发明与现有技术基于时移相位差的频率估计均方根误差(RMSE)曲线图；

图6是时移长度为L＝128时本发明与现有技术基于时移相位差的频率估计均方根误差(RMSE)曲线图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1是本发明实施例中一种基于组合复信号相位差的频率估计方法流程图，如图1所示，所述估计方法包括：

101、通过信号采集器对原始信号采样，获得长度为N+M的数字信号；

假设原始信号是一种单频实正弦信号，该原始信号表示为：

x(t)＝Acos(2πf₀t+θ₀) (1)

其中，A为原始信号的幅值，f₀为信号频率，θ₀为信号的初始相位。原始信号x(t)经采样后得到长度为N+M的时间序列x(n)：

其中，f_s为采样频率，n表示采样点序列，同时，为了满足奈奎斯特采样定理，采样频率满足f_s>2f₀。

其中，在本发明中，N表示采样点数，M表示第二段信号相对于第一段信号的间隔长度，N和M均为正整数。

当然，可以理解的是，原始信号除了采用上述单频实正弦信号以外，还可以采用其他序列信号，本发明对此不作限定。

102、从该段数字信号0点起取N点构造出第一时间序列，从该段数字信号第M点起取N点构造出第二时间序列；

本步骤中，第一时间序列是x(n)中前N点数字信号，即从序列x(n)中的0点到N-1点；第二时间序列是以序列x(n)的0点时间平移M个点作为起始点开始取N点，即第M点到第M+N-1点。

第一时间序列x₀(n)和第二时间序列x₁(n)依次表示如下：

其中，λ₀为归一化频率，其表示为λ₀＝f₀N/f_s。

103、构造出长度为N的对称窗函数，分别对所述第一时间序列和所述第二时间序列进行对称加窗处理得到第一加窗函数和第二加窗函数；

分别对第一时间序列x₀(n)和第二时间序列x₁(n)进行加窗处理，则对应的两个加窗函数可表示为：

x_w0(n)＝x₀(n)w(n),n＝0,1,2,…,N-1 (5)

x_w1(n)＝x₁(n)w(n),n＝0,1,2,…,N-1 (6)

其中，w(n)表示N点对称窗函数。

104、将所述第一加窗函数和所述第二加窗函数构造出第三时间序列，对所述第三时间序列进行离散傅里叶变换得到该序列的频谱；

利用第一加窗函数x_w0(n)和第二加窗函数x_w1(n)构造的新时间序列即第三时间序列表示为：

y(n)＝x_w0(n)+jx_w1(n),n＝0,1,2,…,N-1 (7)

对第三时间序列进行离散傅里叶变换得：

Y(k)＝X_w0(k)+jX_w1(k),k＝0,1,2,…,N-1 (8)

其中，X_w0(k)和X_w1(k)分别为第一加窗函数x_w0(n)和第二加窗函数x_w1(n)的离散傅里叶变换。根据式(5)和式(6)可得：

其中，θ₁＝θ₀+2πλ₀S为x₁(n)的初始相位，S＝M/N表示时移系数，W(·)和τ(·)分别表示对称窗函数w(n)的窗谱函数和相频函数。将式(9)和(10)代入式(8)可得：

其中，Y(k)表示第三时间序列的频谱。

105、利用所述频谱的峰值谱线解算得到第一时间序列与第二时间序列的初始相位差，通过所述相位差估计得到原始信号的频率值。

在实际离散频谱分析中，由于难以实现对信号的整周期截断，以及有限的观测时间，不可避免会造成频谱的泄露和栅栏效应，因此归一化频率λ₀往往不是整数，其可表示为：

λ₀＝l+δ (12)

其中，l为峰值频点λ₀附近抽样幅值最大所对应的谱线号，其为整数。δ为归一化频率校正量，其为小数，取值范围为[-0.5,0.5]。通过频域幅值搜索法得到l,并将其正负频点即k＝±l代入式(11)可得：

分别记第一参量r₁＝Y(l)/Y^*(-l)，第二参量r₂＝Y(-l)/Y^*(l)，Y^*(·)表示Y(·)的共轭复数，则有：

根据欧拉公式，对式(15)的第一参量和式(16)的第二参量可进一步化简为：

其中Δθ＝θ₁-θ₀为相位差。进一步地，设γ＝Im(r₁)，η＝Im(r₂)，Im(·)表示取复数的虚部，则：

根据式(19)和(20)可解算得相位差为：

或

当|Y(l)|＜|Y(-l)|时，相位差取前者，反之取后者。

当正频点l处的幅值小于负频点-l处的幅值时即|Y(l)|＜|Y(-l)|时，解算得到相位差取第一相位差，反之取第二相位差；其中：

所述第一相位差表示为：

所述第二相位差表示为：

其中，Y(l)表示频点l对应的频谱值；Y(-l)表示频点-l对应的频谱值；Y^*(·)表示Y(·)的共轭复数；

在本发明实施例中，通过相位差估计原始信号频率值，具体包括：

考虑到相位绕卷问题，相位差Δθ满足如下关系式：

Δθ+2πh＝2πλ₀S＝2π(l+δ)S, h∈Z (25)

其中，h为相位绕卷次数；则有：

其中，

对式(26)取模得：

{·}表示对1的取模运算；由于S＜1时等式左边取值范围-0.5～0.5，而相位绕卷问题可能导致等式右边超出该取值范围；所以令

并对

做如下调整：

举个例子，假设

为0.7，那么调整后的

则为-0.3，假设

为-0.7，那么调整后的

为0.3。

根据式(28)，可计算得归一化频率校正量表示为：

进一步的，通过第一调整参数ε₁和第二调整参数ε₂对

的符号进行调整：

其中|Y(l)|，|Y(l±1)|，

分别为最高谱线，次高谱线，第三高谱线的幅值。ε₀表示第零调整参数，

sign(·)为符号函数。则调整后的归一化频率校正量为：

式中，ε＝sign(ε₁+ε₂)。则归一化频率估计值为：

在一些实施例中，在实际的工程应用中，所采集到的信号中往往包含较强的随机噪声，这可能导致谱线定位错误，带来较大的校正误差。为了评估噪声对新算法中谱线定位的影响，考虑带有加性高斯白噪声的信号由式(1)给出。设置幅值A＝1，采样频率f_s＝128Hz，数据长度N＝128，故频率分辨率为1Hz。频率在区间[31.5,32.5]内以步距0.025变化，初相θ₀在[-π,π]内随机变化，时移长度M＝128。图2和图3所示为10000次模拟实验中，不同信噪比下加汉宁窗的基于时移相位差校正法和基于组合复信号相位差校正法的谱线定位错误概率曲线。如图2和图3所示，谱线定位错误的概率随着绝对频率偏差|δ|的增大而升高。对于相同的频率偏差，谱线定位错误的概率随信噪比的减小而下降。本发明相较于传统的基于时移的相位差校正法，在谱线定位的准确性上具有明显优势。

在另一些实施例中，为了评估噪声对新算法频率估计性能的影响，考虑带有加性高斯白噪声的信号由式(1)给出。设置幅值A＝1，采样频率f_s＝128Hz，数据长度N＝128，故频率分辨率为1Hz。频率在区间[31.5,32.5]内以步距0.025变化，初相θ₀在[-π,π]内随机变化，信噪比SNR＝0dB，时移长度M＝96,112,128。图4～图6所示为10000次模拟实验中基于时移的相位差校正法和基于组合复信号的相位差校正法的均方根误差(RMSE)曲线。其中，图4中是时移长度为M＝96时本发明与现有技术基于时移相位差的频率估计均方根误差(RMSE)曲线图；图5是时移长度为M＝112时本发明与现有技术基于时移相位差的频率估计均方根误差(RMSE)曲线；图6是时移长度为M＝128时本发明与现有技术基于时移相位差的频率估计均方根误差(RMSE)曲线图；如图4～图6所示，在高斯白噪声的干扰下，由于相位绕卷问题导致传统的基于时移相位差的校正方法在|δ|较大处存在极大频率估计误差。而本发明所提出的基于组合复信号的相位差校正方法具有更小的均方根误差。因此，本发明具有较好的抗噪声性能。

在本发明的描述中，需要理解的是，术语“同轴”、“底部”、“一端”、“顶部”、“中部”、“另一端”、“上”、“一侧”、“顶部”、“内”、“外”、“前部”、“中央”、“两端”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

在本发明中，除非另有明确的规定和限定，术语“安装”、“设置”、“连接”、“固定”、“旋转”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (10)

1.一种基于组合复信号相位差的频率估计方法，其特征在于，包括：

通过信号采集器对原始信号采样，获得长度为N+M的数字信号；

从该段数字信号0点起取N点构造出第一时间序列，从该段数字信号第M点起取N点构造出第二时间序列；

构造出长度为N的对称窗函数，分别对所述第一时间序列和所述第二时间序列进行对称加窗处理得到第一加窗函数和第二加窗函数；

将所述第一加窗函数和所述第二加窗函数构造出第三时间序列，对所述第三时间序列进行离散傅里叶变换得到该序列的频谱；

利用所述频谱的峰值谱线解算得到第一时间序列与第二时间序列的初始相位差，通过所述相位差估计得到原始信号的频率值。

2.根据权利要求1所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述长度为N+M的数字信号表示为：

其中，x(n)表示长度为N+M的数字信号；A为原始信号的幅值，f₀为信号频率，f_s为采样频率，采样频率满足f_s>2f₀；n表示采样点序列，θ₀为信号的初始相位；N表示采样点数，M表示第二时间序列相对于第一时间序列的间隔长度。

3.根据权利要求1所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述对称窗函数为任意对称窗函数。

4.根据权利要求1所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述第三时间序列表示为：

y(n)＝x_w0(n)+jx_w1(n),n＝0,1,2,…,N-1

其中，y(n)表示为第三时间序列；x_w0(n)表示第一加窗函数；x_w1(n)表示第二加窗函数。

5.根据权利要求1所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述利用所述频谱的峰值谱线解算得到相位差包括利用频率幅值搜索得到所述频谱中峰值谱线对应的频点l，将其正负频点±l带入离散傅里叶变换后的第三时间序列中，得到第一参量和第二参量；根据欧拉公式对所述第一参量和所述第二参量化简，得到第一相位差和第二相位差：当正频点l处的幅值小于负频点-l处的幅值时，所述解算得到相位差取第一相位差，反之取第二相位差。

6.根据权利要求5所述的一种基于组合复信号相位差的频率估计方法，其特征在于，

所述第一相位差表示为：

所述第二相位差表示为：

其中，Y(l)表示频点l对应的频谱值；Y(-l)表示频点-l对应的频谱值；Y^*(·)表示Y(·)的共轭复数；Im(·)表示取复数的虚部。

7.根据权利要求1所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述通过所述相位差估计得到原始信号的频率值表示为：

其中，

表示原始信号的频率估计值；l表示频谱中峰值谱线对应的频点；

表示调整后的归一化频率校正量；f_s为采样频率。

8.根据权利要求7所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述调整后的归一化频率校正量的获取过程包括根据时移系数和取模后的相位差计算出归一化频率校正量

按照最高谱线幅值|Y(l)|、次高谱线幅值|Y(l±1)|、第三高谱线幅值

计算出第一调整参数和第二调整参数；根据所述第一调整参数和所述第二调整参数对所述归一化频率校正量的符号进行调整，调整后的归一化频率校正量的公式表示为

ε＝sign(ε₁+ε₂)；ε₁表示第一调整参数；ε₂表示第二调整参数；sign(·)为符号函数。

9.根据权利要求8所述的一种基于组合复信号相位差的频率估计方法，其特征在于，所述归一化频率校正量

的计算过程包括基于相位绕卷构建出相位差等式，对所述相位差等式两边取模，对取模后可变相位进行调整，将调整后的可变相位与时移系数的比值作为归一化频率校正量。

10.根据权利要求7所述的一种基于组合复信号相位差的频率估计方法，其特征在于，

所述第一调整参数表示为：

所述第二调整参数表示为：

其中，ε₀表示第零调整参数，

Images