技术领域本发明属于电力系统谐波信号检测技术领域，尤其涉及一种基于中心频移的谐波参数估计算法。
背景技术：
随着非线性设备使用量的增加，由其引起的谐波污染问题也越来越严重。这些问题不仅威胁到电力系统的安全运行，也缩短了电力系统中设备的使用寿命，所以高精度的谐波分析一直是人们研究的热点。快速傅里叶变化(fastFouriertransform,FFT)是一种应用极为广泛的信号处理方法。在长期的实际应用和理论分析中发现，由FFT计算得到的频谱在估算谐波参数时，会产生极大的误差。减小这种误差使用最为广泛的校正方法是插值法，即寻找频率归一后最接近主瓣峰值的两条谱线，并根据它们的窗谱函数的比值建立一个以频率参数为变量的方程，然后求解校正频率，同时估算出幅值和相角。为了提高插值法的精确度，会选用旁瓣性能更好的窗函数，比如Rife-Vincent窗、Nuttall窗和自卷积窗；或者是参考更多谱线信息。这些改进算法在提高精确度的同时会出现以下问题：(1)算法复杂；无论是改用更高性能的窗函数，还是增加参考谱线的个数，都会加大算法的复杂程度，使得谐波计算时间变长，影响算法的实时性。(2)抗干扰性弱；在实际环境中，以上的改进算法会使算法中受干扰的因素更多，难以保证谐波检测在干扰环境中的可靠性。所以现有的谐波检测环境更需要一种计算简便、抗干扰性强的谐波计算方法。
技术实现要素：
针对现有的谐波分析算法过程复杂和抗干扰性弱的缺点，本发明提出了一种基于中心频移的谐波参数估计算法。本发明提出的算法与传统的加窗插值算法相比，省去了拟合插值的过程，同时利用中心频移的方法提高了算法的抗干扰性，更适用于电力系统中的谐波检测。本发明所采用的技术方案是：一种基于中心频移的谐波参数估计算法，包括以下步骤：步骤1：待测谐波信号s(t)进行采样得到离散信号s(n)，然后对离散信号s(n)加汉宁窗函数，即sw(n)＝s(n)·w(n)，最后sw(n)进行FFT计算得到频谱Sw(k)；步骤2：搜寻基波谱线，判断其位置是否位于量化频点的二分之一；若是，则直接进行后续步骤；若不是，则根据距离差D对整体信号进行时域频移，得到sD(n)，然后经过FFT计算得到新的频谱SD(k)；步骤3：搜寻各次谐波谱峰的位置khm，然后依据各次谐波谱线位于量化频点二分之一处的条件，和同时考虑其对应的偏移因数Eh，对信号的幅值、频率和相角进行估算。计算基波谱线与量化频点二分之一处距离D的公式为：D=12-|Sw(k1m+r)||Sw(k1m+r)|+|Sw(k1m)|]]>其中，|Sw()|为信号谱线的幅值；k1m为基波谱峰对应的位置；当|Sw(k1m+1)|>|Sw(k1m-1)|时，r＝1；当|Sw(k1m+1)|<|Sw(k1m-1)|时，r＝-1。上式中的谱线幅值易于得到，用来估算谐波谱线位置简单可靠。所述的判断基波谱线位置是否位于量化频点的二分之一处的方法是：比较距离D与阈值ε的大小，当D>ε时，可认为基波谱线不位于量化频点的二分之一处，否则条件成立；阈值ε大小可视实际精度要求而定，通过阈值判断频移是否，可以减轻算法的运算量。所述各次谐波对应的偏移因素Eh是：Eh=12-|SD(khm+l)||SD(khm)|+|SD(khm+l)|]]>其中，|SD()|为频移后信号谱线的幅值；khm为第h次谐波谱峰对应的位置；当|SD(khm+1)|>|SD(khm-1)|时，l＝1；当|SD(khm+1)|<|SD(khm-1)|时，l＝-1；估算谐波参数时考虑到偏移因素，可协助校正各次谐波的谱线位置，增加算法的估算精度。所述信号的幅值估算公式为：Ah=2·[|SD(khm)|+|SD(khm+l)|]|W(-12l+l·Eh)|+|W(12l+l·Eh)|]]>其中，|W()|为窗函数频域幅值表达式，|SD()|为频移后信号谱线的幅值；Eh为各次谐波对应的偏移因素；khm为第h次谐波谱峰对应的位置；当|SD(khm+1)|>|SD(khm-1)|时，l＝1；当|SD(khm+1)|<|SD(khm-1)|时，l＝-1。所述信号的频率估算公式为：fh=(khm+l2-l·Eh-r·D)·Δf]]>其中，fs为信号采样频率，N为采样个数，Eh为各次谐波对应的偏移因素。所述信号的相角估算公式为：其中，|SD()|为频移后信号谱线的幅值；Eh为各次谐波对应的偏移因素；khm为第h次谐波谱峰对应的位置。本发明一种基于中心频移的谐波参数估计算法，有益效果如下：1)、实时性高：本发明提出的算法与传统加窗插值的算法相比，省去了拟合插值的步骤，缩短了运算时间，提高了算法在实际应用中的实时性。2)、抗干扰性强：本方发明加入了中心频移的步骤，该步骤可保证双谱线幅值近乎相等，不会存在次谱线幅值过弱造成精度受环境干扰影响较大的情况。3)、精度高：本发明在估算谐波参数的过程中加入了偏移因素的概念，用偏移因素来修正谐波谱线的位置，可以抑制在频移过程中算法受到的影响，从来提高参数估算的精确度。4)、易于嵌入式系统的运用：本发明算法步骤简便，易于编程与应用于嵌入式系统中。附图说明图1为本发明算法的流程图。图2为本发明信号中心频移示意图。图3为本发明白噪声下频率相对误差对比图。图4为本发明白噪声下幅值相对误差对比图。图5为本发明白噪声下相角相对误差对比图。具体实施方式一种基于中心频移的谐波参数估计算法，包括以下步骤：步骤1：待测谐波信号s(t)进行采样得到离散信号s(n)，然后对离散信号s(n)加汉宁窗函数，即sw(n)＝s(n)·w(n)，最后sw(n)进行FFT计算得到频谱Sw(k)；步骤2：搜寻基波谱线，判断其位置是否位于量化频点的二分之一；若是，则直接进行后续步骤；若不是，则根据距离差D对整体信号进行时域频移，得到sD(n)，然后经过FFT计算得到新的频谱SD(k)；步骤3：搜寻各次谐波谱峰的位置khm，然后依据各次谐波谱线位于量化频点二分之一处的条件，和同时考虑其对应的偏移因数Eh，对信号的幅值、频率和相角进行估算。计算基波谱线与量化频点二分之一处距离D的公式为：D=12-|Sw(k1m+r)||Sw(k1m+r)|+|Sw(k1m)|]]>其中，|Sw()|为信号谱线的幅值；k1m为基波谱峰对应的位置；当|Sw(k1m+1)|>|Sw(k1m-1)|时，r＝1；当|Sw(k1m+1)|<|Sw(k1m-1)|时，r＝-1。上式中的谱线幅值易于得到，用来估算谐波谱线位置简单可靠。所述的判断基波谱线位置是否位于量化频点的二分之一处的方法是：比较距离D与阈值ε的大小，当D>ε时，可认为基波谱线不位于量化频点的二分之一处，否则条件成立；阈值ε大小可视实际精度要求而定，通过阈值判断频移是否，可以减轻算法的运算量。所述各次谐波对应的偏移因素Eh是：Eh=12-|SD(khm+l)||SD(khm)|+|SD(khm+l)|]]>其中，|SD()|为频移后信号谱线的幅值；khm为第h次谐波谱峰对应的位置；当|SD(khm+1)|>|SD(khm-1)|时，l＝1；当|SD(khm+1)|<|SD(khm-1)|时，l＝-1；估算谐波参数时考虑到偏移因素，可协助校正各次谐波的谱线位置，增加算法的估算精度。所述信号的幅值估算公式为：Ah=2·[|SD(khm)|+|SD(khm+l)|]|W(-12l+l·Eh)|+|W(12l+l·Eh)|]]>其中，|W()|为窗函数频域幅值表达式，|SD()|为频移后信号谱线的幅值；Eh为各次谐波对应的偏移因素；khm为第h次谐波谱峰对应的位置；当|SD(khm+1)|>|SD(khm-1)|时，l＝1；当|SD(khm+1)|<|SD(khm-1)|时，l＝-1。所述信号的频率估算公式为：fh=(khm+l2-l·Eh-r·D)·Δf]]>其中，fs为信号采样频率，N为采样个数，Eh为各次谐波对应的偏移因素。所述信号的相角估算公式为：其中，|SD()|为频移后信号谱线的幅值；Eh为各次谐波对应的偏移因素；khm为第h次谐波谱峰对应的位置。下面结合附图，对优选实例进行详细说明。本发明估算谐波参数的步骤为：1、谐波信号的采样：假设检测的信号共含H次谐波，其信号表达为：式中：Ah为第h次谐波的幅值；fh为第h次谐波的频率；为第h次谐波的相角。信号经采样后得到离散信号s(n),采样频率为fs，采样点共N个。2、采样信号加窗FFT计算：将离散信号加汉宁窗函数w(n)，即sw(n)＝s(n)w(n),再进行FFT变换得到信号频谱Sw(k)。上式中：为频谱中频点的间隔；W()为对应窗函数离散傅里叶变换后的表达式。式(3)可简化为：3、搜寻基波谱线位置，判断是否频移：搜寻基波谱线的位置k1m，求出基波谱线与量化频点二分之一处的距离D。D=12-|Sw(k1m+r)||Sw(k1m+r)|+|Sw(k1m)|---(4)]]>式中：当|Sw(k1m+1)|>|Sw(k1m-1)|时，r＝1；当|Sw(k1m+1)|<|Sw(k1m-1)|时，r＝-1。阈值ε为判断中心频移的依据，ε的大小可视实际的精度要求而定。若D>ε，对信号进行频移；反之，跳过下一步进行后续计算。4、信号中心频移：对信号进行中心频移：sD(n)=sw(n)·ej2πnNrD---(5)]]>频移过程如图2所示，频移后再通过FFT计算得到新的频谱SD(k)。5、计算偏移因素和信号参数：在频谱中，计算每一次谐波的偏移因素Eh，即各次谐波谱线与量化频点二分之一处的距离差。Eh=12-|SD(khm+l)||SD(khm)|+|SD(khm+l)|---(6)]]>上式中：khm为第h次谐波对应的谱峰位置；l的大小与每次谐波谱线的相对位置有关，当|SD(khm+1)|>|SD(khm-1)|时，l＝1；当|SD(khm+1)|<|SD(khm-1)|时，l＝-1。频率估算公式为：fh=(khm+l2-l·Eh-r·D)·Δf---(7)]]>相角估算公式为：幅值参数可以根据谱峰与次谱峰的幅值平均化求得：Ah=2·[|SD(khm)|+|SD(khm+l)|]|W(-12l+l·Eh)|+|W(12l+l·Eh)|---(9)]]>下面将通过实施例进一步说明本发明技术效果：以采样频率1500Hz采样含10次谐波的电力系统信号s(t),采样点数N＝1024，基波频率f0＝50Hz，各次谐波的幅值和相位见表1。选取窗函数长度N＝1024，对比本发明算法与双谱线插值的分析结果。它们的频率相对误差结果如表2所示，幅值相对误差如表3所示，相角相对误差如表4所示。表1信号模型参数表2频率相对误差对比表3幅值相对误差对比表4相角相对误差对比从仿真结果可知，本发明算法可以有效估算出谐波参数，且整体分析精度要高于双谱线插值算法。为了验证本发明在干扰环境下的分析效果，在上信号模型中加入白噪声干扰，仿真观察本发明算法的分析效果，并与双谱线插值进行对比。仿真结果如图3、图4和图5所示。由仿真结果可知在白噪声环境中，本发明依然可以准确地估算出各次谐波信号的参数，即使信噪比(SNR)在30dB以下，本发明的估算误差也满足于标准要求。从相对误差曲线的光滑程度可知，本发明在噪声环境中稳定性更强。相对于双谱线插值算法而言，本发明受噪声影响小且在精确度上占有优势，更适用于噪声干扰下的谐波检测。