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

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


背景技术：

2.谐波实信号指含有多个不同频率的正弦分量的实信号，谐波实信号的参数估计是从含有噪声的谐波采样信号中检测出谐波各分量的频率、幅度、初相位等参数，这种技术在电力系统监测、雷达阵列处理、生物工程、地球物理等领域有着广泛应用，具有重要的理论研究意义和实际应用价值。
3.近年来，许多学者针对谐波实信号的参数估计方法展开了许多研究，这些方法主要可以分为基于矩阵变换的参数估计方法和基于离散傅里叶变换(dft)的参数估计方法两大类。基于矩阵变换的参数估计方法构造信号子空间和噪声子空间矩阵，通过矩阵的特征值分解或奇异值分解来求解谐波参数，具有超分辨率、精度高等特点，但实现困难且计算量非常大，主要有多重信号分类法、线性预测法等。基于dft的参数估计方法根据傅立叶变换获得信号频谱，并利用快速傅里叶变换(fft)技术实现快速处理，具有实现简单，实时性好等特点，硬件实现相对容易，因而得到了广泛研究，目前主要有加窗插值法、频谱搬移法等。
4.(1)多重信号分类法(参考文献[1]：schmidt r，schmidt r o.multiple emitter location and signal parameter estimation[j].ieee transactions on antennas&propagation，1986，34(3)：276
‑
280.)，该方法通过计算信号的协方差矩阵，构建信号的噪声子空间，通过求解矩阵的特征值和特征向量来计算采样信号的参数。该方法为高分辨率的现代谱估计方法，但矩阵的特征值分解运算量大，计算复杂，实时性差。
[0005]
(2)线性预测法(参考文献[2]：so h c，chan k w，chan y t，et al.linear prediction approach for efficient frequency estimation of multiple real sinusoids：algorithms and analyses[j].ieee transactions on signal processing，2005，53(7)：2290
‑
2305.)，该方法利用谐波实信号的线性预测性质估计信号参数，并利用迭代计算进一步提高精度。但该方法中需要矩阵的求逆运算，时间复杂度较大，同时迭代计算增加了算法复杂度。
[0006]
(3)加窗插值法(参考文献[3]：zeng b，teng z，cai y l，guo s y and qing b y.harmonic phasor analysis based on improved fft algorithm[j].ieee transactions on smart grid，2011，2(1)：51
‑
59.)，该方法首先对谐波采样信号进行加窗处理以抑制频谱泄漏，然后通过多项式拟合得到插值公式，利用dft频谱谱线插值方法快速计算频率。该方法实现速度快，原理简单，但窗函数的选择受限，无法完全消除频谱泄漏带来的影响，估计精度难以提高。
[0007]
(4)频谱搬移法(参考文献[4]：djukanovic s，popovic
‑
bugarin v.efficient and accurate detection and frequency estimation of multiple sinusoids[j].ieee access，2019，7：1118
‑
1125.)，该方法首先获得谐波实信号频率的粗估计，再利用频谱搬移技术，针对某一谐波分量的估计，以滤除直流分量的形式消除其他分量的影响，然后利用三
点抛物线插值获得频率估计。但是频谱搬移难以完全滤除粗估计得到的谐波分量，使得算法性能与估计精度难以提高。


技术实现要素：

[0008]
本发明旨在提出一种估计精度更高、抗噪性更好、实时性强、易于硬件实现的信号参数估计方法，适用于含噪声的谐波实信号参数估计，解决现有频域法受谐波实信号各分量的正频率与负频率频谱泄露相互影响的问题。
[0009]
本发明所提具体方法说明如下：
[0010]
方法的基本思想：在估计某一谐波实信号分量时，将其他分量视为干扰分量并去除，通过计算正交信号来生成此分量的解析信号，并利用迭代插值算法估计当前分量的各项参数，对谐波中所含的每个分量均重复以上步骤，直到谐波实信号的每个分量的参数均得到计算。
[0011]
主要包括以下步骤：首先，通过fft粗估计谐波实信号各分量中的显著分量，即能量最大的信号频率，生成此分量的正交信号并生成解析信号，对解析信号使用迭代插值算法粗估计频率，在迭代过程中更新解析信号；接着，利用解析信号与频率粗估计计算当前分量的幅度、初相位，根据参数粗估计将当前能量最大的分量从原始信号中去除，重复以上操作直到所有谐波分量的参数粗估计计算完成；然后，针对某一谐波分量，从原始信号中利用参数粗估计将其他分量去除，生成此分量的解析信号并利用迭代插值算法计算参数精估计，重复此操作直到所有谐波分量的参数精估计计算完成；最后，重复整个精估计操作，进一步提高谐波实信号的参数估计精度。
[0012]
设谐波实信号的采样信号模型如式(1)所示
[0013][0014]
式中，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，方差为σ2的高斯白噪声。对采样谐波实信号的第k个分量，其信噪比定义为snr
k
＝a2/2σ2。
[0015]
对采样信号参数估计的频域分析而言，频率估计是最重要的步骤，信号的幅度、相位参数均可在频率估计的基础上得到。针对采样谐波信号，第k阶谐波频率ω
k
可表示为
[0016][0017]
式中，m
k
为采样信号频谱中第k个能量极大值点的索引值，m
k
＝round(ω
k
n/2π)，round(
·
)表示对
·
进行四舍五入取整；δ
k
为频谱偏移量，δ
k
∈(
‑
0.5，0.5)。
[0018]
为抑制谐波实信号各分量的正频率与负频率的相互影响，提高谐波实信号参数估计精度，基于上述思想和信号模型，提出一种解析信号生成的谐波实信号参数估计方法，主要包括参数粗估计和参数精估计两部分，原理图如图1所示。
[0019]
(1)参数粗估计，原理图如图2所示，主要包括以下步骤：
[0020]
第一步：初始化变量。
[0021]
对含有k个谐波分量的实信号x(n)，为了不破坏原始信号，设置临时变量s(n)＝x
(n)，并初始化频率残差初值：δ
1，2，..，k
＝0，设置粗估计迭代次数q1。
[0022]
第二步：估计能量最大分量的频谱索引值与对应频率。
[0023]
利用式(3)对s(n)进行快速傅里叶变换，估计当前临时变量s(n)中能量最大的频率分量索引值m
k
，k＝1，2，...，k，并利用式(2)计算频率粗估计ω
k(c)
；
[0024]
x(m)＝fft[s(n)]m＝0，1，
…
，n
‑1ꢀꢀꢀꢀꢀꢀꢀ
(3)
[0025][0026]
第三步：生成当前索引值m
k
对应谐波分量的解析信号。
[0027]
首先，利用式(5)计算正交信号生成参数b：
[0028]
b＝round(π/2ω
k(c)
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0029]
需要注意的是，若计算得到b＝0，则令b＝1。
[0030]
然后，利用式(6)计算正交信号
[0031][0032]
最后，利用式(7)计算解析信号s
(j)
(n)
[0033][0034]
第四步：根据解析信号利用迭代插值算法进行频率估计。
[0035]
首先，利用式(8)计算解析信号的频谱插值。
[0036][0037]
然后，利用式(9)更新频率残差
[0038][0039]
最后，利用式(2)更新频率粗估计值。
[0040]
第五步：迭代计算第二步至第四步q1次，更新解析信号，提高频率估计精度。
[0041]
第六步：消除当前计算的第k阶谐波分量的影响。
[0042]
首先，利用式(10)计算复幅值的粗估计值
[0043][0044]
然后，利用式(11)计算幅值和相角
[0045]
a
k(c)
＝|a
k(c)
|，φ
k(c)
＝∠a
k(c)
ꢀꢀꢀꢀꢀꢀꢀ
(11)
[0046]
式中，|
·
|表示取复数
·
的绝对值，∠
·
表示取复数
·
的相角。
[0047]
最后，利用式(12)更新临时变量s(n)，去除第k阶分量的影响。
[0048]
s(n)＝s(n)
‑
a
k(c)
cos(ω
k(c)
n+φ
k(c)
)
ꢀꢀꢀꢀꢀꢀꢀ
(12)
[0049]
第七步：循环计算第五步至第六步k次，最终得到所有k阶谐波分量的参数粗估计值
[0050]
(2)参数精估计，原理图如图3所示，主要包括以下步骤：
[0051]
第一步：参数初始化。
[0052]
令ω
k
＝ω
k(c)
，a
k
＝a
k(c)
，φ
k
＝φ
k(c)
，k＝1，2，...，k，其中，ω
k
，a
k
，φ
k
分别为第k阶谐波分量的参数精估计值，设置精估计迭代次数q2。
[0053]
第二步：去除干扰分量，生成解析信号。
[0054]
首先，对第k阶谐波，利用式(13)计算仅含第k阶谐波分量的序列s
k
(n)。
[0055][0056]
接着，利用式(5)，根据ω
k
计算解析信号生成参数b。
[0057]
然后，利用式(6)，将s(n)替换为s
k
(n)，计算s
k
(n)的正交信号。
[0058]
最后，利用式(7)生成s
k
(n)的解析信号
[0059]
第三步：根据解析信号进行参数估计。
[0060]
首先，利用式(14)计算解析信号的频谱插值。
[0061][0062]
接着，利用式(15)更新频率精估计值ω
k
。
[0063][0064]
然后，利用式(16)计算复幅值a
k
。
[0065][0066]
最后，利用式(17)计算幅值和相角的精估计值。
[0067]
a
k
＝|a
k
|，φ
k
＝∠a
k
ꢀꢀꢀꢀꢀꢀꢀ
(17)
[0068]
第四步：循环执行第二步至第四步k次，得到所有k阶谐波分量的参数精估计值；
[0069]
第五步：循环执行第四步q2次，更新参数精估计值，提高参数估计精度。
[0070]
相对于现有的技术，本发明的有益效果如下：本发明基于快速傅里叶变换实现，通过解析信号生成与迭代插值算法，在估计某一谐波分量时去除其他分量，具有实时性高，易于嵌入式系统的实现的特点，且本发明抗干扰性强，参数估计误差接近克拉美劳下限，在中高信噪比下估计性能好，比现有方法的精度更高。
附图说明
[0071]
图1是解析信号生成的谐波实信号参数估计方法的基本原理。
[0072]
图2是本方法参数粗估计部分的原理图。
[0073]
图3是本方法参数精估计部分的原理图。
具体实施方式
[0074]
本发明的具体实施方式如下：
[0075]
(1)参数粗估计
[0076]
第一步：对含有k个谐波分量的实信号x(n)，初始化临时变量s(n)＝x(n)，初始化频率残差初值：设置粗估计迭代次数q1；
[0077]
第二步：利用式x(m)＝fft[s(n)]m＝0，1，
…
，n
‑
1对s(n)进行快速傅里叶变换，并由式得到实信号频谱能量最大值的索引值m
k
，利用计算频率粗估计ω
k(c)
；
[0078]
第三步：利用b＝round(π/2ω
k(c)
)计算正交信号生成参数b，利用式：
[0079][0080]
计算正交信号利用式计算解析信号s
(1)
(n)；
[0081]
第四步：利用式p＝
±
0.5计算解析信号的频谱插值，利用更新频率残差，利用更新频率粗估计值；
[0082]
第五步：迭代计算第二步至第四步q1次，并利用更新频率粗估计值；
[0083]
第六步：利用式计算复幅值的粗估计值，利用式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)；
[0084]
第七步：循环计算第五步和第六步k次，最终得到所有k个谐波分量的参数粗估计值。
[0085]
(2)参数精估计：
[0086]
第一步：令ω
k
＝ω
k(c)
，a
k
＝a
k(c)
，φ
k
＝φ
k(c)
，k＝1，2，...，k，设置迭代次数q2；
[0087]
第二步：对第k阶谐波，利用式计算仅含第k阶谐波分量的序列s
k
(n)，利用式b＝round(π/2ω
k
)计算解析信号生成参数b，利用式：
[0088][0089]
计算s
k
(n)的正交信号，利用生成s
k
(n)的解析信号；
[0090]
第三步：利用式p＝
±
0.5计算解析信号的频谱插
值，利用式更新频率精估计值，利用式计算复幅值的精估计值，利用a
k
＝|a
k
|，φ
k
＝∠a
k
计算幅值和相角的精估计值；
[0091]
第四步：循环执行第三步和第四步k次，得到所有k阶谐波分量的参数精估计值；
[0092]
第五步：循环执行第四步q2次，更新参数精估计值，提高参数估计精度。