本發(fā)明具體涉及一種數(shù)字化諧波分析方法。

背景技術(shù)：

隨著國家經(jīng)濟技術(shù)的發(fā)展和人們生活水平的提高，數(shù)字化技術(shù)已經(jīng)深入人心，而且也已經(jīng)廣泛應(yīng)用于人們的生產(chǎn)和生活當(dāng)中。

數(shù)字化變電站是智能電網(wǎng)發(fā)展的趨勢，相對于傳統(tǒng)變電站，數(shù)據(jù)采集數(shù)字化是數(shù)字化變電站技術(shù)應(yīng)用的標(biāo)志之一。根據(jù)IEC61850-9-2的規(guī)定，合并單元的采樣頻率被固定成了幾個數(shù)據(jù)，按照國內(nèi)工頻50HZ，每周波常用的點數(shù)轉(zhuǎn)換成了40，80，200，其中最高是200點。由于每周波常用的點數(shù)越多，則數(shù)字化采樣得到的波形越準(zhǔn)確，因此200點已經(jīng)成為了每周波常用點數(shù)的趨勢。按照電能質(zhì)量標(biāo)準(zhǔn)《IEC61000-4-30》的要求，諧波計算窗口變成了10周波，對應(yīng)的采樣點就成了2000點。由于2000點并不滿足2的整數(shù)次冪，無法直接使用FFT進行諧波分析。在這種情況下，可以采用的方法有DFT算法，插值成2的整數(shù)次冪的插值方法，基于因數(shù)合成原理的混合基算法，素因子算法（簡稱PFA）、Winograd算法（簡稱WFTA算法）等等，但單個這些方法的運算速度都不是很高效，達不到實際工程應(yīng)用中快速運算的目的，所以現(xiàn)有的諧波分析方法還存在著實際應(yīng)用的一些問題。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的在于提供一種計算速度極高的數(shù)字化諧波分析方法。

本發(fā)明提供的這種數(shù)字化諧波分析方法，包括如下步驟：

S1. 對需要進行諧波分析的對象進行采樣，得到n點的實數(shù)序列[x₁,x₂,x₃,……x_n]；

S2. 對得到的X點實數(shù)序列[x₁,x₂,x₃,……x_n]進行分割，并對分割后的序列分別采用混合基算法、PFA算法和WFTA算法進行諧波分析，完成數(shù)字化諧波分析。

所述的數(shù)字化諧波分析方法，在步驟S1后，步驟S2之前還包括如下步驟：

S1.5 對得到的n點的實數(shù)序列[x₁,x₂,x₃,……x_n]，根據(jù)如下算式將序列[x₁,x₂,x₃,……x_n]轉(zhuǎn)換為復(fù)數(shù)序列Y[y₁,y_2,y₃,……,y_m]：

y₁=x₁+j*x₂，y₂=x₃+j*x₄，y₃=x₅+j*x₆，……，y_i=x_2i-1+j*x_2i，……，y_m=x_n-1+j*x_n；其中m=n/2；

在上述步驟下，步驟S2所述的對得到的n點實數(shù)序列[x₁,x₂,x₃,……x_n]進行分割，即為對m個點的復(fù)數(shù)序列進行分割。

所述的對得到的m點復(fù)數(shù)序列進行分割，具體包括如下步驟：

Ⅰ. 根據(jù)WFTA算法的特性和m的值，選定WFTA算法最佳的計算點數(shù)A；

Ⅱ. 根據(jù)PFA算法的特性和m的值，選定PFA算法的最佳計算點數(shù)B；

Ⅲ. 根據(jù)混合基計算算法的特性和m的值，選定混合基計算算法的最佳計算點數(shù)C；

Ⅳ. 對選定的計算點數(shù)A、B和C，采用如下算式進行檢驗：

B=N1*A；

C=N2*B；

m=N3*C；

N2=N4*A；

式中N1、N2、N3和N4均為自然數(shù)；

Ⅴ. 若檢驗成功，則將m點復(fù)數(shù)序列分割成為N3組C個點數(shù)的序列，并將C個點的序列分割成為N2個B個點數(shù)的序列，再將B個點的序列分割成為N1個A個點數(shù)的序列。

所述的對分割后的序列分別采用混合基算法、PFA算法和WFTA算法進行諧波分析，具體包括如下步驟：

ⅰ. 將m個點數(shù)的序列分解為N3組C個點數(shù)的序列，進行N3次，每次C個點的混合基計算；

ⅱ. 將C個點數(shù)的序列分解成N2組B個點數(shù)的序列，進行N2次，每次B個點的PFA計算；

ⅲ. 將B個點數(shù)的序列分解為N1組A個點數(shù)的序列，進行N1次，每次A個點的WFTA計算；

ⅳ. 將步驟ⅲ計算后的B個點數(shù)序列進行轉(zhuǎn)換，將轉(zhuǎn)換后的B個點數(shù)的序列分解為A個N1個點的序列，再進行A次，每次N1個點的WFTA計算；

ⅴ. 將步驟ⅳ計算后的B個點數(shù)序列進行譯序，完成B個點數(shù)序列的PFA計算；再將步驟ⅱ~步驟ⅴ得到的C個點數(shù)的序列乘以C個點的旋轉(zhuǎn)因子；

ⅵ. 將步驟ⅴ得到的旋轉(zhuǎn)后的C個點的序列分解為B組N2個點數(shù)的序列，再進行B次，每次N2個點PFA計算；

ⅶ. 將N2個點的序列分解為N4組A個點數(shù)的序列，進行N4次，每次A個點的WFTA計算；

ⅷ. 將步驟ⅶ計算后的N2個點的序列轉(zhuǎn)換，再把轉(zhuǎn)換后的N2個點的序列分解為A組N4個點數(shù)的序列，再進行A次，每次N4個點的WFTA計算；

ⅸ. 將步驟ⅷ計算得到的N2個點數(shù)的序列進行N2點PFA譯序，完成N2個點序列的PFA計算；

ⅹ. 將步驟ⅱ~步驟ⅸ計算得到的C個點數(shù)的序列進行譯序，完成C個點數(shù)的混合基計算；

xi. 將步驟ⅰ~步驟ⅹ計算得到的m個點數(shù)的序列進行旋轉(zhuǎn)；

xⅱ. 將步驟xi得到的m個點數(shù)的序列分解為C組每組N3個點數(shù)的序列，再進行C次，每次N3個點數(shù)的WFTA計算；

xⅲ. 對步驟xⅱ計算得到的m個點數(shù)的序列進行譯序，完成m個點數(shù)的復(fù)數(shù)序列分析；

xⅳ. 將步驟xⅲ得到的m個點數(shù)的復(fù)數(shù)序列轉(zhuǎn)換為n個點的實數(shù)序列，并根據(jù)該實數(shù)序列進行諧波幅值和相位的計算，完成數(shù)字化諧波分析。

所述的對分割后的序列分別采用混合基算法、PFA算法和WFTA算法進行諧波分析，在步驟ⅰ之前還包括如下步驟：

O. 根據(jù)m個點數(shù)復(fù)數(shù)序列的分割結(jié)果，將FFT算法、混合基算法、PFA算法和WFTA算法中需要用到的常數(shù)和常量事先計算并存儲，以提高本發(fā)明方法的計算速度。

所述的FFT算法、混合基算法、PFA算法和WFTA算法中需要用到的常數(shù)和常量包括m點組合的FFT譯序數(shù)組，C點混合基譯序數(shù)組，B點PFA譯序數(shù)組、N2點PFA譯序數(shù)組，A點WFTA三角函數(shù)常量、N4點WFTA三角函數(shù)常量、N2點WFTA三角函數(shù)常量，C個點的旋轉(zhuǎn)因子，m個點的旋轉(zhuǎn)因子和m點復(fù)數(shù)轉(zhuǎn)實數(shù)的三角函數(shù)常量。

所述的n取值為2000；m取值為1000；N3取值為5；C取值為200；N2取值為10；B取值為20；N1取值為4；A取值為5；N4取值為2。

本發(fā)明提供的這種數(shù)字化諧波分析方法，算法層次清晰：組合算法層次清晰，順序執(zhí)行，充分發(fā)揮實數(shù)分解、混合基、PFA、WFTA四個算法的優(yōu)點，揚長避短；而且運算量小，計算過程中需要的譯序數(shù)組和旋轉(zhuǎn)因子提前計算好，減少運算量；同時測量精度高，達到國標(biāo)諧波A級計算要求；此外，本方法運算速度快，按照工控機MCU的主頻300MHZ，計算速度可以達到10ms，遠遠低于其他單獨算法的幾十和幾百ms，工程應(yīng)用效果好，達到諧波實時快速分析的目的，給數(shù)字化變電站的諧波分析提供強有力的技術(shù)保證。

附圖說明

圖1為本發(fā)明方法的具體實施流程示意圖。

具體實施方式

如圖1所示為本發(fā)明方法的具體實施流程示意圖：本發(fā)明提供的這種數(shù)字化諧波分析方法，首先對需要進行諧波分析的對象進行采樣，得到n點的實數(shù)序列[x₁,x₂,x₃,……x_n]；然后對得到的X點實數(shù)序列[x₁,x₂,x₃,……x_n]進行分割，并對分割后的序列分別采用混合基算法、PFA算法和WFTA算法進行諧波分析，完成數(shù)字化諧波分析。

以下結(jié)合附圖1對本發(fā)明方法進行詳細說明：以按照電能質(zhì)量標(biāo)準(zhǔn)《IEC61000-4-30》的要求，諧波計算窗口變成了10周波，每個周波采樣200個點，對應(yīng)的采樣點就成了2000點；選定WFTA算法最佳的計算點數(shù)A為5點，選定PFA算法的最佳計算點數(shù)B為20點，選定混合基計算算法的最佳計算點數(shù)C為200點；

（1）將FFT算法、混合基算法、PFA算法和WFTA算法中需要用到的常數(shù)和常量事先計算并存儲，以提高本發(fā)明方法的計算速度： FFT算法、混合基算法、PFA算法和WFTA算法中需要用到的常數(shù)和常量，包括1000點組合的FFT譯序數(shù)組，200點混合基譯序數(shù)組，20點PFA譯序數(shù)組，10點PFA譯序數(shù)組，5點WFTA三角函數(shù)常量，4點WFTA三角函數(shù)常量、2點WFTA三角函數(shù)常量、200個點的旋轉(zhuǎn)因子，1000個點的旋轉(zhuǎn)因子和1000點復(fù)數(shù)轉(zhuǎn)實數(shù)的三角函數(shù)常量，；

（2）將2000點的實數(shù)序列[x₁,x₂,x₃,……x₂₀₀₀]變成1000點的復(fù)數(shù)序列[y₁,y_2,y₃,……,y₁₀₀₀]；

y₁=x₁+j*x₂，y₂=x₃+j*x₄，y₃=x₅+j*x₆，……，y_i=x_2i-1+j*x_2i，……，y₁₀₀₀=x1₉₉₉+j*x₂₀₀₀；

（3）將1000點復(fù)數(shù)序列分解成5*200，M=5，L=200，M是循環(huán)次數(shù)，進行循環(huán)5次，每次200點的混合基計算；

（4）將200點序列分解成10*20，P=10，Q=20，P是循環(huán)次數(shù)，進行循環(huán)10次，每次20點序列的PFA計算；

（5）將20點序列分解成4*5，進行循環(huán)4次，每次5點序列的WFTA計算；

（6）將（5）計算后的20點序列行列轉(zhuǎn)換后，再把得到的20點序列又分解成5*4，進行循環(huán)5次，每次4點序列的WFTA計算；

（7）將（6）計算后的序列進行20點序列的譯序，使用初始化好的譯序數(shù)組，完成20點序列的PFA計算，再將從步驟（4）到此得到的200點序列乘以200點序列初始化好的旋轉(zhuǎn)因子；

（8）與步驟（4）相對應(yīng)，將（7）計算出的200點序列分解成20*10，Q=20，P=10，Q為循環(huán)次數(shù)，進行循環(huán)20次，每次10點序列的PFA計算；

（9）將10點序列分解成2*5，進行循環(huán)2次，每次5點的WFTA計算；

（10）將（9）計算后的10點序列行列轉(zhuǎn)換后，再把得到的10點序列分解成5*2，進行循環(huán)5次，每次2點的WFTA計算；

（11）將（10）計算后的序列進行10點PFA的譯序，使用初始化好的譯序數(shù)組，完成10點序列的PFA計算；

（12）對從（4）到（11）計算出的200點序列進行譯序，使用初始化好的譯序數(shù)組，完成200點的混合基計算；

（13）將從（3）到此得到的1000點序列乘以初始化好的1000點旋轉(zhuǎn)因子；

（14）將（13）計算出的1000點序列又分解成200*5，L=200，M=5，L是循環(huán)次數(shù)，進行循環(huán)200次，每次5點的WFTA計算；

（15）對（14）計算出的1000個序列使用初始化好的1000點譯序數(shù)組進行譯序，完成1000點的復(fù)數(shù)序列分析；

（16）將（15）得到的1000點復(fù)數(shù)序列使用初始化好的三角函數(shù)進行2000點的實數(shù)序列轉(zhuǎn)換；

（17）根據(jù)分解出的2000點的實部和虛部，進行各次諧波幅值和相位的計算，完成2000點序列的諧波分析。