CN112816779B - 一种解析信号生成的谐波实信号参数估计方法 - Google Patents

Abstract

本发明涉及信号处理领域，特别是含噪谐波实信号的参数估计方法。本发明的适用对象为含噪谐波实信号的参数估计，主要包括以下步骤：首先，通过FFT粗估计采样信号的显著分量，生成此分量的解析信号并使用迭代插值算法粗估计频率；接着，计算当前分量的幅度、初相位粗估计并将当前显著分量从原始信号中去除，重复以上操作直到得到所有分量的参数粗估计；然后，对每个谐波分量，将其他分量去除后生成此分量的解析信号并计算参数精估计，最终得到所有谐波分量的参数精估计；最后，重复整个精估计操作，进一步提高参数估计精度。本发明涉及的谐波实信号的参数估计方法原理清晰、抗噪性好、实时性强，提高了参数估计精度。

Description

一种解析信号生成的谐波实信号参数估计方法

技术领域

本发明涉及信号处理领域，特别是含噪谐波实信号的参数估计方法。

背景技术

谐波实信号指含有多个不同频率的正弦分量的实信号，谐波实信号的参数估计是从含有噪声的谐波采样信号中检测出谐波各分量的频率、幅度、初相位等参数，这种技术在电力系统监测、雷达阵列处理、生物工程、地球物理等领域有着广泛应用，具有重要的理论研究意义和实际应用价值。

近年来，许多学者针对谐波实信号的参数估计方法展开了许多研究，这些方法主要可以分为基于矩阵变换的参数估计方法和基于离散傅里叶变换(DFT)的参数估计方法两大类。基于矩阵变换的参数估计方法构造信号子空间和噪声子空间矩阵，通过矩阵的特征值分解或奇异值分解来求解谐波参数，具有超分辨率、精度高等特点，但实现困难且计算量非常大，主要有多重信号分类法、线性预测法等。基于DFT的参数估计方法根据傅立叶变换获得信号频谱，并利用快速傅里叶变换(FFT)技术实现快速处理，具有实现简单，实时性好等特点，硬件实现相对容易，因而得到了广泛研究，目前主要有加窗插值法、频谱搬移法等。

(1)多重信号分类法(参考文献[1]：Schmidt R，Schmidt R O.Multiple emitterlocation and signal parameter estimation[J].IEEE Transactions on Antennas&Propagation，1986，34(3)：276-280.)，该方法通过计算信号的协方差矩阵，构建信号的噪声子空间，通过求解矩阵的特征值和特征向量来计算采样信号的参数。该方法为高分辨率的现代谱估计方法，但矩阵的特征值分解运算量大，计算复杂，实时性差。

(2)线性预测法(参考文献[2]：So H C，Chan K W，Chan Y T，et al.Linearprediction approach for efficient frequency estimation of multiple realsinusoids：algorithms and analyses[J].IEEE Transactions on Signal Processing，2005，53(7)：2290-2305.)，该方法利用谐波实信号的线性预测性质估计信号参数，并利用迭代计算进一步提高精度。但该方法中需要矩阵的求逆运算，时间复杂度较大，同时迭代计算增加了算法复杂度。

(3)加窗插值法(参考文献[3]：Zeng B，Teng Z，Cai Y L，Guo S Y and Qing BY.Harmonic Phasor Analysis Based on Improved FFT Algorithm[J].IEEETransactions on Smart Grid，2011，2(1)：51-59.)，该方法首先对谐波采样信号进行加窗处理以抑制频谱泄漏，然后通过多项式拟合得到插值公式，利用DFT频谱谱线插值方法快速计算频率。该方法实现速度快，原理简单，但窗函数的选择受限，无法完全消除频谱泄漏带来的影响，估计精度难以提高。

(4)频谱搬移法(参考文献[4]：Djukanovic S，Popovic-Bugarin V.Efficientand Accurate Detection and Frequency Estimation of Multiple Sinusoids[J].IEEEAccess，2019，7：1118-1125.)，该方法首先获得谐波实信号频率的粗估计，再利用频谱搬移技术，针对某一谐波分量的估计，以滤除直流分量的形式消除其他分量的影响，然后利用三点抛物线插值获得频率估计。但是频谱搬移难以完全滤除粗估计得到的谐波分量，使得算法性能与估计精度难以提高。

发明内容

本发明旨在提出一种估计精度更高、抗噪性更好、实时性强、易于硬件实现的信号参数估计方法，适用于含噪声的谐波实信号参数估计，解决现有频域法受谐波实信号各分量的正频率与负频率频谱泄露相互影响的问题。

本发明所提具体方法说明如下：

方法的基本思想：在估计某一谐波实信号分量时，将其他分量视为干扰分量并去除，通过计算正交信号来生成此分量的解析信号，并利用迭代插值算法估计当前分量的各项参数，对谐波中所含的每个分量均重复以上步骤，直到谐波实信号的每个分量的参数均得到计算。

主要包括以下步骤：首先，通过FFT粗估计谐波实信号各分量中的显著分量，即能量最大的信号频率，生成此分量的正交信号并生成解析信号，对解析信号使用迭代插值算法粗估计频率，在迭代过程中更新解析信号；接着，利用解析信号与频率粗估计计算当前分量的幅度、初相位，根据参数粗估计将当前能量最大的分量从原始信号中去除，重复以上操作直到所有谐波分量的参数粗估计计算完成；然后，针对某一谐波分量，从原始信号中利用参数粗估计将其他分量去除，生成此分量的解析信号并利用迭代插值算法计算参数精估计，重复此操作直到所有谐波分量的参数精估计计算完成；最后，重复整个精估计操作，进一步提高谐波实信号的参数估计精度。

设谐波实信号的采样信号模型如式(1)所示

式中，K为谐波阶数，a_k＞0为第k阶谐波的幅度，θ_k∈(-π，π)为第k阶谐波的初相位，ω_k＝2πf_k/f_s为第k阶谐波的角频率，f_k为第k阶谐波的模拟频率，f_s为信号采样频率，由f_s＞2f_k得ω_k∈(0，π)，N为信号长度，n表示信号的索引值；w(n)是均值为0，方差为σ²的高斯白噪声。对采样谐波实信号的第k个分量，其信噪比定义为SNR_k＝a²/2σ²。

对采样信号参数估计的频域分析而言，频率估计是最重要的步骤，信号的幅度、相位参数均可在频率估计的基础上得到。针对采样谐波信号，第k阶谐波频率ω_k可表示为

式中，m_k为采样信号频谱中第k个能量极大值点的索引值，m_k＝round(ω_kN/2π)，round(·)表示对·进行四舍五入取整；δ_k为频谱偏移量，δ_k∈(-0.5，0.5)。

为抑制谐波实信号各分量的正频率与负频率的相互影响，提高谐波实信号参数估计精度，基于上述思想和信号模型，提出一种解析信号生成的谐波实信号参数估计方法，主要包括参数粗估计和参数精估计两部分，原理图如图1所示。

(1)参数粗估计，原理图如图2所示，主要包括以下步骤：

第一步：初始化变量。

对含有K个谐波分量的实信号x(n)，为了不破坏原始信号，设置临时变量s(n)＝x(n)，并初始化频率残差初值：δ_{1，2，..，K}＝0，设置粗估计迭代次数Q₁。

第二步：估计能量最大分量的频谱索引值与对应频率。

利用式(3)对s(n)进行快速傅里叶变换，估计当前临时变量s(n)中能量最大的频率分量索引值m_k，k＝1，2，...，K，并利用式(2)计算频率粗估计ω_k ^(c)；

X(m)＝FFT[s(n)]m＝0，1，…，N-1 (3)

第三步：生成当前索引值m_k对应谐波分量的解析信号。

首先，利用式(5)计算正交信号生成参数b：

b＝round(π/2ω_k ^(c)) (5)

需要注意的是，若计算得到b＝0，则令b＝1。

然后，利用式(6)计算正交信号

最后，利用式(7)计算解析信号s^(j)(n)

第四步：根据解析信号利用迭代插值算法进行频率估计。

首先，利用式(8)计算解析信号的频谱插值。

然后，利用式(9)更新频率残差

最后，利用式(2)更新频率粗估计值。

第五步：迭代计算第二步至第四步Q₁次，更新解析信号，提高频率估计精度。

第六步：消除当前计算的第k阶谐波分量的影响。

首先，利用式(10)计算复幅值的粗估计值

然后，利用式(11)计算幅值和相角

a_k ^(c)＝|A_k ^(c)|，φ_k ^(c)＝∠A_k ^(c) (11)

式中，|·|表示取复数·的绝对值，∠·表示取复数·的相角。

最后，利用式(12)更新临时变量s(n)，去除第k阶分量的影响。

s(n)＝s(n)-a_k ^(c)cos(ω_k ^(c)n+φ_k ^(c)) (12)

第七步：循环计算第五步至第六步K次，最终得到所有K阶谐波分量的参数粗估计值

(2)参数精估计，原理图如图3所示，主要包括以下步骤：

第一步：参数初始化。

令ω_k＝ω_k ^(c)，a_k＝a_k ^(c)，φ_k＝φ_k ^(c)，k＝1，2，...，K，其中，ω_k，a_k，φ_k分别为第k阶谐波分量的参数精估计值，设置精估计迭代次数Q₂。

第二步：去除干扰分量，生成解析信号。

首先，对第k阶谐波，利用式(13)计算仅含第k阶谐波分量的序列s_k(n)。

接着，利用式(5)，根据ω_k计算解析信号生成参数b。

然后，利用式(6)，将s(n)替换为s_k(n)，计算s_k(n)的正交信号。

最后，利用式(7)生成s_k(n)的解析信号

第三步：根据解析信号进行参数估计。

首先，利用式(14)计算解析信号的频谱插值。

接着，利用式(15)更新频率精估计值ω_k。

然后，利用式(16)计算复幅值A_k。

最后，利用式(17)计算幅值和相角的精估计值。

a_k＝|A_k|，φ_k＝∠A_k (17)

第四步：循环执行第二步至第四步K次，得到所有K阶谐波分量的参数精估计值；

第五步：循环执行第四步Q₂次，更新参数精估计值，提高参数估计精度。

相对于现有的技术，本发明的有益效果如下：本发明基于快速傅里叶变换实现，通过解析信号生成与迭代插值算法，在估计某一谐波分量时去除其他分量，具有实时性高，易于嵌入式系统的实现的特点，且本发明抗干扰性强，参数估计误差接近克拉美劳下限，在中高信噪比下估计性能好，比现有方法的精度更高。

附图说明

图1是解析信号生成的谐波实信号参数估计方法的基本原理。

图2是本方法参数粗估计部分的原理图。

图3是本方法参数精估计部分的原理图。

具体实施方式

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

(1)参数粗估计

第一步：对含有K个谐波分量的实信号x(n)，初始化临时变量s(n)＝x(n)，初始化频率残差初值：设置粗估计迭代次数Q₁；

第二步：利用式X(m)＝FFT[s(n)]m＝0，1，…，N-1对s(n)进行快速傅里叶变换，并由式得到实信号频谱能量最大值的索引值m_k，利用/>计算频率粗估计ω_k ^(c)；

第三步：利用b＝round(π/2ω_k ^(c))计算正交信号生成参数b，利用式：

计算正交信号利用式/>计算解析信号s⁽¹⁾(n)；

第四步：利用式p＝±0.5计算解析信号的频谱插值，利用/>更新频率残差，利用/>更新频率粗估计值；

第五步：迭代计算第二步至第四步Q₁次，并利用更新频率粗估计值；

第六步：利用式计算复幅值的粗估计值，利用式a_k ^(c)＝|A_k ^(c)|，φ_k ^(c)＝∠A_k ^(c)计算幅值和相角，利用式s(n)＝s(n)-a_k ^(c)cos(ω_k ^(c)n+φ_k ^(c))更新临时变量s(n)；

第七步：循环计算第五步和第六步K次，最终得到所有K个谐波分量的参数粗估计值。

(2)参数精估计：

第一步：令ω_k＝ω_k ^(c)，a_k＝a_k ^(c)，φ_k＝φ_k ^(c)，k＝1，2，...，K，设置迭代次数Q₂；

第二步：对第k阶谐波，利用式计算仅含第k阶谐波分量的序列s_k(n)，利用式b＝round(π/2ω_k)计算解析信号生成参数b，利用式：

计算s_k(n)的正交信号，利用生成s_k(n)的解析信号；

第三步：利用式p＝±0.5计算解析信号/>的频谱插值，利用式/>更新频率精估计值，利用式/>计算复幅值的精估计值，利用a_k＝|A_k|，φ_k＝∠A_k计算幅值和相角的精估计值；

第四步：循环执行第三步和第四步K次，得到所有K阶谐波分量的参数精估计值；

第五步：循环执行第四步Q₂次，更新参数精估计值，提高参数估计精度。

Claims (1)

1.一种解析信号生成的谐波实信号参数估计方法，其特征在于：适用对象为含噪谐波实信号的参数估计，该方法包括以下步骤：

(1)参数粗估计：

第一步：对含有K个谐波分量的实信号序列x(n)，设置临时变量s(n)＝x(n)，初始化频率残差初值：δ_k＝0，k＝1，2，...，K，设置粗估计迭代次数Q₁；

式中：s(n)为临时变量，n表示s(n)时刻点，n＝0，1，…，N-1，N表示信号长度，K为谐波所含分量个数；

第二步：利用式X(m)＝FFT[s(n)] m＝0，1，…，N-1对s(n)进行快速傅里叶变换，并由式得到实信号频谱能量最大值的索引值，利用/>计算频率粗估计ω_k ^(c)；

式中：FFT(·)表示对·进行快速傅里叶变换，arg max X(m)表示X(m)取最大值时k的取值，m_k为频谱能量最大值的索引值，ω_k ^(c)表示第k阶谐波的频率粗估计；

第三步：利用b＝round(π/2ω_k ^(c))计算解析信号生成参数b，利用式：

计算正交信号利用式/>计算解析信号s⁽¹⁾(n)；

式中：round(·)表示对·进行四舍五入取整；

第四步：利用式计算解析信号的频谱插值，利用更新频率残差，利用/>更新频率粗估计值；

式中：S_p表示信号的频谱值，p＝±0.5表示插值点间隔；

第五步：迭代计算第二步至第四步Q₁次，并利用更新频率粗估计值；

第六步：利用式计算复幅值的粗估计值，利用式a_k ^(c)＝|A_k ^(c)|，φ_k ^(c)＝∠A_k ^(c)计算幅值和相角，利用式s(n)＝s(n)-a_k ^(c)cos(ω_k ^(c)n+φ_k ^(c))更新临时变量s(n)；

式中：A_k ^(c)表示复幅值的粗估计值，|·|表示取复数·的绝对值，∠·表示取复数·的相角；

第七步：循环计算第五步和第六步K次，最终得到所有K阶谐波分量的参数粗估计值；

(2)参数精估计：

第一步：令ω_k＝ω_k ^(c)，a_k＝a_k ^(c)，φ_k＝φ_k ^(c)，k＝1，2，...，K，设置精估计迭代次数Q₂；

式中，ω_k，a_k，φ_k分别为第k阶谐波分量频率、幅度、相位的精估计值；

第二步：对第k阶谐波，利用式计算仅含第k阶谐波分量的序列s_k(n)，利用式b＝round(π/2ω_k)计算解析信号生成参数b，利用式：

计算s_k(n)的正交信号，利用生成s_k(n)的解析信号；

式中：s_k(n)表示仅含第k阶谐波分量的序列，表示s_k(n)的正交信号，/>表示s_k(n)的解析信号；

第三步：利用式计算解析信号/>的频谱插值，利用式/>更新频率精估计值，利用式/>计算复幅值的精估计值，利用a_k＝|A_k|，φ_k＝∠A_k计算幅值和相角的精估计值；

第四步：循环执行第三步和第四步K次，得到所有K阶谐波分量的参数精估计值；

第五步：循环执行第四步Q₂次，更新参数精估计值，提高参数估计精度。