数字化谐波分析方法技术领域
本发明具体涉及一种数字化谐波分析方法。
背景技术
随着国家经济技术的发展和人们生活水平的提高，数字化技术已经深入人心，而
且也已经广泛应用于人们的生产和生活当中。
数字化变电站是智能电网发展的趋势，相对于传统变电站，数据采集数字化是数
字化变电站技术应用的标志之一。根据IEC61850-9-2的规定，合并单元的采样频率被固定
成了几个数据，按照国内工频50HZ，每周波常用的点数转换成了40，80，200，其中最高是200
点。由于每周波常用的点数越多，则数字化采样得到的波形越准确，因此200点已经成为了
每周波常用点数的趋势。按照电能质量标准《IEC61000-4-30》的要求，谐波计算窗口变成了
10周波，对应的采样点就成了2000点。由于2000点并不满足2的整数次幂，无法直接使用FFT
进行谐波分析。在这种情况下，可以采用的方法有DFT算法，插值成2的整数次幂的插值方
法，基于因数合成原理的混合基算法，素因子算法（简称PFA）、Winograd算法（简称WFTA算
法）等等，但单个这些方法的运算速度都不是很高效，达不到实际工程应用中快速运算的目
的，所以现有的谐波分析方法还存在着实际应用的一些问题。
发明内容
本发明的目的在于提供一种计算速度极高的数字化谐波分析方法。
本发明提供的这种数字化谐波分析方法，包括如下步骤：
S1. 对需要进行谐波分析的对象进行采样，得到n点的实数序列[x1,x2,x3,……xn]；
S2. 对得到的X点实数序列[x1,x2,x3,……xn]进行分割，并对分割后的序列分别采用
混合基算法、PFA算法和WFTA算法进行谐波分析，完成数字化谐波分析。
所述的数字化谐波分析方法，在步骤S1后，步骤S2之前还包括如下步骤：
S1.5 对得到的n点的实数序列[x1,x2,x3,……xn]，根据如下算式将序列[x1,x2,
x3,……xn]转换为复数序列Y[y1,y2,y3,……,ym]：
y1=x1+j*x2，y2=x3+j*x4，y3=x5+j*x6，……，yi=x2i-1+j*x2i，……，ym=xn-1+j*xn；其中m=n/
2；
在上述步骤下，步骤S2所述的对得到的n点实数序列[x1,x2,x3,……xn]进行分割，即为
对m个点的复数序列进行分割。
所述的对得到的m点复数序列进行分割，具体包括如下步骤：
Ⅰ. 根据WFTA算法的特性和m的值，选定WFTA算法最佳的计算点数A；
Ⅱ. 根据PFA算法的特性和m的值，选定PFA算法的最佳计算点数B；
Ⅲ. 根据混合基计算算法的特性和m的值，选定混合基计算算法的最佳计算点数C；
Ⅳ. 对选定的计算点数A、B和C，采用如下算式进行检验：
B=N1*A；
C=N2*B；
m=N3*C；
N2=N4*A；
式中N1、N2、N3和N4均为自然数；
Ⅴ. 若检验成功，则将m点复数序列分割成为N3组C个点数的序列，并将C个点的序列分
割成为N2个B个点数的序列，再将B个点的序列分割成为N1个A个点数的序列。
所述的对分割后的序列分别采用混合基算法、PFA算法和WFTA算法进行谐波分析，
具体包括如下步骤：
ⅰ. 将m个点数的序列分解为N3组C个点数的序列，进行N3次，每次C个点的混合基计算；
ⅱ. 将C个点数的序列分解成N2组B个点数的序列，进行N2次，每次B个点的PFA计算；
ⅲ. 将B个点数的序列分解为N1组A个点数的序列，进行N1次，每次A个点的WFTA计算；
ⅳ. 将步骤ⅲ计算后的B个点数序列进行转换，将转换后的B个点数的序列分解为A个
N1个点的序列，再进行A次，每次N1个点的WFTA计算；
ⅴ. 将步骤ⅳ计算后的B个点数序列进行译序，完成B个点数序列的PFA计算；再将步骤
ⅱ~步骤ⅴ得到的C个点数的序列乘以C个点的旋转因子；
ⅵ. 将步骤ⅴ得到的旋转后的C个点的序列分解为B组N2个点数的序列，再进行B次，每
次N2个点PFA计算；
ⅶ. 将N2个点的序列分解为N4组A个点数的序列，进行N4次，每次A个点的WFTA计算；
ⅷ. 将步骤ⅶ计算后的N2个点的序列转换，再把转换后的N2个点的序列分解为A组N4
个点数的序列，再进行A次，每次N4个点的WFTA计算；
ⅸ. 将步骤ⅷ计算得到的N2个点数的序列进行N2点PFA译序，完成N2个点序列的PFA计
算；
ⅹ. 将步骤ⅱ~步骤ⅸ计算得到的C个点数的序列进行译序，完成C个点数的混合基计
算；
xi. 将步骤ⅰ~步骤ⅹ计算得到的m个点数的序列进行旋转；
xⅱ. 将步骤xi得到的m个点数的序列分解为C组每组N3个点数的序列，再进行C次，每
次N3个点数的WFTA计算；
xⅲ. 对步骤xⅱ计算得到的m个点数的序列进行译序，完成m个点数的复数序列分析；
xⅳ. 将步骤xⅲ得到的m个点数的复数序列转换为n个点的实数序列，并根据该实数序
列进行谐波幅值和相位的计算，完成数字化谐波分析。
所述的对分割后的序列分别采用混合基算法、PFA算法和WFTA算法进行谐波分析，
在步骤ⅰ之前还包括如下步骤：
O. 根据m个点数复数序列的分割结果，将FFT算法、混合基算法、PFA算法和WFTA算法中
需要用到的常数和常量事先计算并存储，以提高本发明方法的计算速度。
所述的FFT算法、混合基算法、PFA算法和WFTA算法中需要用到的常数和常量包括m
点组合的FFT译序数组，C点混合基译序数组，B点PFA译序数组、N2点PFA译序数组，A点WFTA
三角函数常量、N4点WFTA三角函数常量、N2点WFTA三角函数常量，C个点的旋转因子，m个点
的旋转因子和m点复数转实数的三角函数常量。
所述的n取值为2000；m取值为1000；N3取值为5；C取值为200；N2取值为10；B取值为
20；N1取值为4；A取值为5；N4取值为2。
本发明提供的这种数字化谐波分析方法，算法层次清晰：组合算法层次清晰，顺序
执行，充分发挥实数分解、混合基、PFA、WFTA四个算法的优点，扬长避短；而且运算量小，计
算过程中需要的译序数组和旋转因子提前计算好，减少运算量；同时测量精度高，达到国标
谐波A级计算要求；此外，本方法运算速度快，按照工控机MCU的主频300MHZ，计算速度可以
达到10ms，远远低于其他单独算法的几十和几百ms，工程应用效果好，达到谐波实时快速分
析的目的，给数字化变电站的谐波分析提供强有力的技术保证。
附图说明
图1为本发明方法的具体实施流程示意图。
具体实施方式
如图1所示为本发明方法的具体实施流程示意图：本发明提供的这种数字化谐波
分析方法，首先对需要进行谐波分析的对象进行采样，得到n点的实数序列[x1,x2,x3,……
xn]；然后对得到的X点实数序列[x1,x2,x3,……xn]进行分割，并对分割后的序列分别采用混
合基算法、PFA算法和WFTA算法进行谐波分析，完成数字化谐波分析。
以下结合附图1对本发明方法进行详细说明：以按照电能质量标准《IEC61000-4-
30》的要求，谐波计算窗口变成了10周波，每个周波采样200个点，对应的采样点就成了2000
点；选定WFTA算法最佳的计算点数A为5点，选定PFA算法的最佳计算点数B为20点，选定混合
基计算算法的最佳计算点数C为200点；
（1）将FFT算法、混合基算法、PFA算法和WFTA算法中需要用到的常数和常量事先计算并
存储，以提高本发明方法的计算速度： FFT算法、混合基算法、PFA算法和WFTA算法中需要用
到的常数和常量，包括1000点组合的FFT译序数组，200点混合基译序数组，20点PFA译序数
组，10点PFA译序数组，5点WFTA三角函数常量，4点WFTA三角函数常量、2点WFTA三角函数常
量、200个点的旋转因子，1000个点的旋转因子和1000点复数转实数的三角函数常量，；
（2）将2000点的实数序列[x1,x2,x3,……x2000]变成1000点的复数序列[y1,y2,y3,……,
y1000]；
y1=x1+j*x2，y2=x3+j*x4，y3=x5+j*x6，……，yi=x2i-1+j*x2i，……，y1000=x1999+j*x2000；
（3）将1000点复数序列分解成5*200，M=5，L=200，M是循环次数，进行循环5次，每次200
点的混合基计算；
（4）将200点序列分解成10*20，P=10，Q=20，P是循环次数，进行循环10次，每次20点序列
的PFA计算；
（5）将20点序列分解成4*5，进行循环4次，每次5点序列的WFTA计算；
（6）将（5）计算后的20点序列行列转换后，再把得到的20点序列又分解成5*4，进行循环
5次，每次4点序列的WFTA计算；
（7）将（6）计算后的序列进行20点序列的译序，使用初始化好的译序数组，完成20点序
列的PFA计算，再将从步骤（4）到此得到的200点序列乘以200点序列初始化好的旋转因子；
（8）与步骤（4）相对应，将（7）计算出的200点序列分解成20*10，Q=20，P=10，Q为循环次
数，进行循环20次，每次10点序列的PFA计算；
（9）将10点序列分解成2*5，进行循环2次，每次5点的WFTA计算；
（10）将（9）计算后的10点序列行列转换后，再把得到的10点序列分解成5*2，进行循环5
次，每次2点的WFTA计算；
（11）将（10）计算后的序列进行10点PFA的译序，使用初始化好的译序数组，完成10点序
列的PFA计算；
（12）对从（4）到（11）计算出的200点序列进行译序，使用初始化好的译序数组，完成200
点的混合基计算；
（13）将从（3）到此得到的1000点序列乘以初始化好的1000点旋转因子；
（14）将（13）计算出的1000点序列又分解成200*5，L=200，M=5，L是循环次数，进行循环
200次，每次5点的WFTA计算；
（15）对（14）计算出的1000个序列使用初始化好的1000点译序数组进行译序，完成1000
点的复数序列分析；
（16）将（15）得到的1000点复数序列使用初始化好的三角函数进行2000点的实数序列
转换；
（17）根据分解出的2000点的实部和虚部，进行各次谐波幅值和相位的计算，完成2000
点序列的谐波分析。