一种高精度的简谐信号参数估计方法与流程

本发明属于信号处理领域，更具体地，涉及一种高精度的简谐信号参数估计方法。

背景技术：

信号的频率参数在工程实际中占有重要的地位，精确的参数估计不仅可以给理论分析提供依据，还可以验证理论分析，因而精确的参数估计得到众多科学家的重视。比如在基于频率的故障诊断领域，裂纹扩展会导致频率发生改变，在某些情况下，由裂纹扩展导致的频率改变量比较小，特别是第一阶固有频率的改变量。为了准确反映裂纹的状况，就需要准确的频率测量，因而频率精度的高低直接决定着诊断的可靠性。众所周知，对加矩形窗的信号进行离散傅里叶变换，其频率误差可达频率分辨率的二分之一，其幅值误差可达36.3％，相位误差为±90°。为了解决上述问题，前人做了不懈的努力，提出多种频率估计方法。这些方法大致可以分为四类，比值法(插值法)，相位差法，能量重心法和fft+ft连续细化分析傅里叶变换法。

传统理论认为在不考虑噪声的情况下，比值法和相位差法是精确的频率估计方法，而能量重心法和fft+ft连续细化分析傅里叶变换法是精度较高的频率估计方法。根据我们的分析，比值法和相位差法仅仅考虑正频率，忽略了负频率，因而得到的结果不是频率的精确估计。严格意义上来说，这四种方法均不是精确的离散频谱矫正方法，很重要的一个原因是离散窗函数的相位计算公式存在着误差，现有技术大都存着这上述问题，这是误差来源之一；离散傅里叶变换得到的是双边谱，而很多研究者仅仅把离散傅里叶变换得到的频谱作为单边谱，这是误差来源之二；对于fft+ft连续细化分析傅里叶变换法来说，必须考虑直流对细化谱的影响，该方法忽略了这一影响，这是误差来源之三。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种高精度的简谐信号参数估计方法，通过矫正关键常数c，同时在计算过程中考虑负频率分量和直流分量，避免简化计算带来的误差，由此解决现有单一简谐信号频率估计值、相位估计值和幅值估计值精度低等问题。

为实现上述目的，按照本发明，提供了一种高精度的简谐信号参数估计方法，其特征在于，该方法包括下列步骤：

(a)按照采样频率fs采集待求解简谐信号的离散信号，在该离散信号上加矩形窗函数获得长度为n的信号x(n)，然后对该信号x(n)进行傅里叶变换获得该信号的频谱x(f)，其中，f是频率；

(b)计算所述频谱x(f)的幅值a(f)和相位p(f)，然后按照所述赋值将所述频谱x(f)分为三个部分，分别为：第一部分为a(1)，第二部分为a(2)～a(n/2)，第三部分为a(n/2+1)～a(n)，在所述第二部分选取两个点分别作为第一频率点f1和第二频率点f2，获取第一频率点的相位和幅值a1，获取第二频率点的相位和幅值a2；

(c)利用所述第一频率点和第二频率点，按照下列关系式(一)，构造所述简谐信号的频率、相位和幅值与两个频率点对应的方程，

其中，a0是所述简谐信号的幅值，f0是所述简谐信号的频率，是所述简谐信号的相位，i是虚数单位，c是常数，取值为-(n-1)π/n。

进一步优选地，在步骤(b)中，所述第一频率点优选为所述第二部分中幅值最大的点，所述第二频率点优选分别为所述第二部分中幅值第二大的点。

进一步优选地，在步骤(c)中，按照所述表达式(一)计算所述简谐信号的频率时，化简后，所述简谐信号的频率优选按照下列表达式(二)或(三)进行：

sinc(f1-f0)+sinc(f1+f0)＝ri1(sinc(f2-f0)-sinc(f2+f0))(二)

sinc(f1-f0)-sinc(f1+f0)＝ri2(sinc(f2-f0)+sinc(f2+f0))(三)

其中，ri1和ri2均为中间变量，f1和f2分别是所述第一频率点和第二频率点的频率，a1和a2分别是所述第一频率点和第二频率点的幅值，和分别是所述第一频率点和第二频率点的相位。

进一步优选地，利用所述表达式(二)或(三)进计算所述简谐信号的频率时，优选在区间[f1-0.5，f1+0.5]的解作为所述简谐信号的频率f0的解。

进一步优选地，在步骤(c)中，按照所述表达式(一)计算所述简谐信号的相位时，化简后，所述简谐信号的相位优选按照下列表达式(四)或(五)进行：

其中，和分别是所述第一频率点和第二频率点的相位。

进一步优选地，在步骤(c)中，按照所述表达式(一)计算所述简谐信号的幅值时，化简后，所述简谐信号的幅值优选按照下列表达式(六)、(七)、(八)或(九)进行：

进一步优选地，在步骤(c)中，当获得的所述简谐信号的幅值为负值时，优选将其相位和幅值进行下列变换：

a’0＝-a0

其中，所述是变换后的所述简谐信号的相位，a’0是变换后的所述简谐信号的幅值。

进一步优选地，在步骤(c)后，利用所述简谐信号计算其直流分量dc优选按照下列表达式：

其中，x(1)是所述第一部分中a(1)对应的频谱。

进一步优选地，在步骤(c)后，优选按照下列表达式将所述简谐信号的频率的单位进行转换：

其中，f’0是所述简谐信号的单位进行转换后的频率，其单位为hz，f0的单位为bin。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1、本发明通过将窗函数的离散傅里叶变换中的常数c的值设定为-(n-1)π/n，频率估计更加准确，任意一个信号的傅里叶变换都是信号与窗函数的卷积，因而傅里叶变换与窗函数密切相关，而窗函数的离散傅里叶变换(dft)和连续傅里叶变换(cft)计算公式不相同，在cft中，常数c的值为-π，在dft中，常数c的值为-(n-1)π/n，在实际傅里叶变换中，傅里叶变换都是dft变换，因而将常数设定为-(n-1)π/n；

2、本发明通过考虑负频率分量，使得频率估计更加准确。该负频率分量的引入是由于：任意一个简谐实信号，在频率域都有两个频率分量，一个正频率分量，一个负频率分量，采用离散傅里叶变换得到的是双边谱，有些窗函数的旁瓣迅速衰减，由负频率导致的谱间干涉非常小，因而通常将负频率分量导致的谱间干涉忽略，其本质是近似计算，达到简化计算的目的，但是，近似计算虽然可以简化计算，但是引入了误差，导致计算精度降低了，本发明没有忽略负频率分量，提高了计算精度。

附图说明

图1是按照本发明的优选实施例所构建的高精度的简谐信号参数估计方法的流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

图1是按照本发明的优选实施例所构建的高精度的简谐信号参数估计方法的流程图，如图1所示，一种高精度的简谐信号参数估计方法，具体包括以下步骤：

步骤一，离散傅里叶变换

将以采样频率为fs获取得到的离散信号加矩形窗函数获得信号x(n)，然后进行离散傅里叶变换，得到其傅里叶变换为x(f)，假设信号x(n)长度为n，那么其离散傅里叶变换后得到的x(f)长度也为n，f的取值为整数，1，2，3，…，n。

步骤二，选择两个频率点

计算步骤(1)求解得到离散傅里叶变换x(f)的模，并乘以2/n得到a(f)，计算步骤一求解得到离散傅里叶变换x(f)的相位p(f)。按照先后顺序将a(f)分成三个部分，第一部分为a(1)，第二部分为a(2)～a(n/2)，第三部分为a(n/2+1)～a(n)，从第二部分选择幅值最大的点作为第一个频率点f1，那么a1＝a(f1)，然后选择幅值第二大的点作为第二频率点f2，通常该点在f1附近，那么a2＝a(f2)，

步骤三，

i、计算参数

常数c的取值为-(n-1)π/n。其中上式的n为离散信号x(n)的长度，计算参数ri1和ri2，ri1和ri2的计算公式如下：

ii、求解频率

按照下列关系式之一求解频率：

sinc(f1-f0)+sinc(f1+f0)＝ri1(sinc(f2-f0)-sinc(f2+f0))

sinc(f1-f0)-sinc(f1+f0)＝ri2(sinc(f2-f0)+sinc(f2+f0))

理论上f0有多个解，选择在区间[f1-0.5，f1+0.5]的解作为频率f0的解，即求解得到频率f0。

iii、求解相位

相位的计算公式有两个，如下两式所示，任意选择一个公式，即可得到相位

最后将求解得到的相位转换到[-π，π]之间，具体转换计算公式为下式，其中mod表示计算余数，具体计算为用第一个数除以第二个数，计算其余数，

iv、求解幅值

幅值a0的计算公式有四个，如下四式所示，任意选择一个公式，即可得到幅值a0，

v、相位和幅值a0的调整

判断所求的幅值a0是否为负值，如果所求得到的幅值为负值，则更新相位和幅值a0，相位和幅值a0的计算公式为：

a’0＝-a0

其中，mod表示计算余数，具体计算为用第一个数除以第二个数，计算其余数，更新后，相位和幅值分别为和a0’。

vi、计算直流分量

直流分量dc的具体计算公式为：

vii、计算频率

实际频率的计算值为：

上式中为fs采样频率，n为采样点数，至此整个计算过程完成，计算得到的频率f’0、幅值a’0、相位和直流分量dc均为数值解。

本研究在前人的基础上，给出了窗函数新的离散傅里叶变换计算公式，并在此基础上提出了一种数值法，用于求解简谐信号的参数。经过验证，在不考虑噪声的情况下，频率、幅值和相位的误差非常小。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。