一种基于方波傅里叶级数分解的新型谐波分析方法
【专利摘要】本发明公开一种基于方波傅里叶级数分解的新型谐波分析方法，利用计算机，通过程序，先输入待检测谐波分析次数N和待分析信号，根据待分析次数产生相应频率的矩形波周期函数，将各矩形波与待分析信号相乘所得的函数作为被积函数，再将其积分，由于矩形波的值只有1或者?1，被积函数大大简化且在采样时相比正余玄函数能极大的节省时间。积分得到N个与各次谐波系数对应的参数，该N个参数构成的矩阵与系数矩阵相乘即得到各次谐波系数，其中系数矩阵由数学分析得到并由N决定且与待分析信号无关，一旦N确定，则系数矩阵不变。
【专利说明】
一种基于方波傅里叶级数分解的新型谐波分析方法
技术领域
[0001] 本发明属于电力系统谐波分析建模技术领域，具体涉及电力系统中谐波系数实时 计算建模方法。
【背景技术】
[0002] 随着电子技术与现代工业的发展，电力系统中非线性负荷不断增加，电力电子器 件的使用不断增加，使得电力系统谐波问题越来越严重，而且许多精密电子工业对电能质 量要求偏高，故及时准确掌握电网中谐波状况对现代工业的发展已经对电力系统的稳定运 行具有极大的意义。
[0003] 电力系统中谐波来源通常来讲有三个方面。一是发电机绕组不是绝对对称，铁芯 也不是绝对均匀，这导致在发电过程中产生少量谐波。二是来自电力变压器等输配电设备 产生的谐波。三是来自于系统中的非线性负荷。
[0004] 电力系统的谐波分析有多种算法，例如小波分析算法、神经网络算法、自适应谐波 算法等，但通常采用快速傅里叶变换(FFT)算法。其中快速傅里叶虽然实现简单使用方便， 但计算量大，实时性欠佳。小波分析法是近几年的热点领域，它是基于傅里叶变换的处理方 法，其优点在于分析暂态谐波时实时性比较好，但小波基的选取难度较大。神经网络法应用 于谐波检测的谐波源识别、谐波预测、谐波测量，但其依赖于算法的准确性与实现难易程 度，有的算法还存在收敛速度及局部最小点等问题，这些都会对检测的精度及实时性产生 影响。自适应谐波检测方法动态响应较慢。
【发明内容】
[0005] 本发明的目的是针对现有谐波检测方法的不足，提供一种简便新颖且具有高实时 性的谐波检测办法，并且该方法具有普遍适用性。
[0006] 为实现本发明目的而采用的技术方案是这样的，一种基于方波傅里叶级数分解的 新型谐波分析方法，其特征在于，包括以下步骤：
[0007] 1)输入待检信号：
[0008]
;其中，Ai和Bi为基波 振幅，Aq和Bq为谐波振幅；cot为相位；
[0009] 所述待检信号含有基波和(1个谐波，基波次数为m，谐波的次数分别为n2、n3…… n q；
[0010] 令系数 ki = ni，i = l、2......q;
[0011 ]根据公式 ki= (2mi-l)ru，计算出系数:mi、m2......mq;
[0012] 生成矩阵I:
[0013]首先，生成一个qXq的全0矩阵，此全0矩阵中，坐标为(ni，m* (2nu-l))的元素，若 Ili X ( 2lHi-1 ) Climax j
[0014] 生成矩阵II:
[0015]首先，生成一个qXq的全0矩阵，此全0矩阵中，坐标为(ni，m* (2nu-l))的元素，若 Ili X ( 2lHi-1 ) Climax
[0016] 2)生成矩形波Sanl( ?t)、sanq( c〇t)、Sbni( ?t)和sbnq( ?t):
[0021] 3)采样：
[0022] 设定步长St = 2Vrw/200，在信号u( ? t)上设定nmax X 200个采样点，每个采样点 的值为u( ? tf)，f = 1，2，3 ? ? ? nmaxX 200 ;nmax为次数最高的谐波的次数；
[0023] 设定步长 St = 23l/nmax/200,在矩形波 Sanl( 〇 t)、Sanq( 〇 t)、Sbnl( W t)和 Sbnq( W t)上 分别设定nmaxX200个采样点，每个采样点的值为Sam( 〇tf)、Sanq( 〇tf)、Sbni( 〇tf)和Sbnq (? tf)，f = 1，2，3 ? ? ? nmaxX 200 ;nmax为次数最高的谐波的次数；
[0024] 4)计算参数 asni、bsni、asnq 和 bsnq:
[0025] asnl = j。丨)，先求解被积函数；
[0026] 将对应采样点u(〇tf)和Sanl(〇tf)的值相乘，得到a nl(f)，其中f=l，2,3 ? ? ? 200Xnmax，令hani(f)=ani(f) ? St，asni = hani(l)+hani(2)+hani(3)+ ? ? ? +hani(200X Hmax)；
[0027] =^：u{c〇nShn[(a^)d(a)i)
[0028] 将对应采样点u(〇tf)和Sbnl(〇tf)的值相乘，得到b nl(f)，其中f=l，2,3 ? ? ? 200Xnmax，令hbni(f)=bni(f) ? St，bsni = hbni(l)+hbni(2)+hbni(3)+ ? ? ? +hbni(200X Hmax)；
[0029] = _fn "("斗(⑴〇6/(欣），先求解被积函数；
[0030] 将对应采样点u( cotf)和Sanq( cotf)的值相乘，得到anq(f)，其中f=l，2,3 ? ? ? 200Xnmax，令hanq(f)=anq(f) ? St，asnq = hanq(l)+hanq(2)+hanq(3)+ ? ? ? +hanq(200X Hmax)；
[0031] h两.=￡ w(故(祕)^/(如），先求解被积函数；
[0032] 将对应采样点u( cotf)和Sbnq( cotf)的值相乘，得到bnq(f)，其中f=l，2,3 ? ? ? 200Xnmax，令hbnq(f)=bnq(f) .St，bsnq = hbnq(l)+hbnq(2)+hbnq(3)+. ? .+hbnq(200X Ilmax)；
[0033] 5)分析得到基波和各次谐波正、余弦分量系数：
[0036] 其中，ai为基波余弦分量系数、a2、a3......aq为谐波余弦分量系数、bi为基波正弦分 量系数、b2、b3......bq为谐波正弦分量系数。
[0037] 进一步，根据步骤5)获得的分析结果，进行谐波的实时检测。
[0038] 值得说明的是，本发明可以利用计算机，通过程序，先输入待检测谐波分析次数N 和待分析信号，将各矩形波与待分析信号相乘所得的函数作为被积函数，再将其积分，由于 矩形波的值只有1或者-1，被积函数大大简化且在采样时相比正余玄函数能极大的节省时 间。积分得到N个与各次谐波系数对应的参数，该N个参数构成的矩阵与系数矩阵相乘即得 到各次谐波系数，其中系数矩阵由数学分析得到并由N决定且与待分析信号无关，一旦N确 定，则系数矩阵不变。所述具体数学分析如下：
[0039]傅里叶级数的展开表达式为：
[0044] 本发明是用特定的矩形波代替上述被积函数u( ? t)cos n? t中的余弦函数cos n ?t。矩形波San如图1所示：
[0045] 其数学表达式为：
[0047]现推导其傅里叶级数：
[0048] 已知在区间[_1，1]上有：
[0052]依据该式子可以推导矩形波傅里叶级数：
[0066]现将谐波系数 中的cosnot替换为San的表达式，SP式 (2)，并做推导：
[0077]其中矩阵I为系数矩阵，各个元素的值由式(7)中&"的系数决定，且该矩阵与输入 信号无关，只与需要计算的谐波次数N有关，N-旦确定，该系数矩阵则确定不变。同理由式 (8)可求得bJ^bsn相对应的系数矩阵矩阵II。对式(9)稍作改变得：
[0079]
，可以用程序直接计算结果，则&"的求解过 程简化为两个数字矩阵的乘积，且式中San( ? t)为矩形波函数，其值只有1或者-1，使得被积 函数大大简化，节省了计算量。流程图如图3。
[0080] 本发明采用上述技术方案后，主要有以下效果：
[0081 ] 1.本发明方法能够准确快速计算任意次谐波幅值，奇次偶次谐波能独立计算，能 够很容易的提取某次谐波参数进行分析。
[0082] 2.本发明方法具有通用性，只要提出需要检测的谐波次数，则可利用计算机程序 算出系数矩阵1、11，并结合输入信号即可快速计算出谐波信息。方法简单，实用性强，便于 推广应用。
[0083] 3.可以高效快速实现输入信号谐波信息的实时检测。
[0084] 本发明计算结果精确，计算速度快，效率高，可计算到任意高次谐波，广泛应用于 电力系统中电能质量分析，为谐波治理打下可靠的基础。
【附图说明】
[0085] 图1为矩形波San;
[0086]图2为矩形波Sbn;
[0087]图3为本发明的流程图；
[0088]图4为积分计算示意图；
[0089]图5为实时谐波检测示意图。
[0090]图中：。
【具体实施方式】
[0091] 下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅 限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯 用手段，做出各种替换和变更，均应包括在本发明的保护范围内。
[0092] 实施例1:
[0093] (1)输入待测函数
[0094] 指定一个信号，只含基波，3、5、7、9次谐波，如下
[0095] u( ? t) = (cos ? t+0.6sin ? t) +
[0096] (0?lcos3 ? t+0?9sin3 ? t) + (0?2cos5 ? t+0?5sin5 ? t) +
[0097] (0?3cos7 ? t+0?3sin7 ? t) + (0?6cos9 ? t+0?2sin9 ? t)
[0098] (11)
[0099] 其中〇^=[0,231]，设定步长5丨=231/9/200，（注:最高为9次谐波，即步长设定为9 次谐波周期除以200)，计算出输入信号在每个采样点的值。
[0100]根据式(7)，有k= (2m-l )n，由于最高算到9次谐波，k的取值为k = 1，3，5，7，9，可得 m、n的取值分别为11 = 1，3,5,7,9;111=1，2,3,4,5。根据式(7)，即可计算出系数矩阵如下：
[0102] (2)周期矩形波的计算
[01 03]根据第（1)步计算得出的n的取值，在Mat lab中生成n = 1，3，5，7，9次的矩形波San (?〇、^(〇^)(注:11=1时矩形波周期为231，11 = 3时周期为231/3，依次类推）。同样根据步骤 (1)中设定的步长，采样出5个不同周期矩形波在每个采样点的值，这些值均为1或-1。
[0104] (3)计算 asn
[0105] 第⑵步完成后，开始计算asn的值。根据式(7)中(叫<;/(w/)，积分 计算具体步骤如下：
[0106] 以求asn为例，先求解被积函数。将u( ?t)和San( ?t)对应采样点的值相乘，令a(k) =u( 〇t(k)) ? San( 〇t(k))(其中k为采样点编号，n = l，3,5,7,9)。令h(k)=a(k) ? St，即h (k)为如下图4阴影部分所示。再将h(k)从k=l，2,3* ? ? 1800求和即可算出asn的值，SPasn = h(l)+h(2)+h(3)+ ? ? ? +11(1800)。&如和 San(?t(k))--对应，从 n = l，3,5,7,9，SalSA 计算求得asi，Sa2带入计算求得a S2，依次类推。
[0107] (4)各次谐波余弦分量系数&"的计算
[0108]第(3)步计算出asn的值后与第（1)步系数矩阵相乘即可得出各次谐波中正弦量与 余弦量的系数。即：
[0110] 用相同的步骤可以计算出bsn及其对应的系数矩阵jbsn然后计算出bn的值，计算结 果如下表：
[0111] 表1算例信号中各次谐波正余弦分量系数
[0112]
[0113] 实施例2:
[0114] 本实施例涉及20次谐波
[0115] (1)输入需计算谐波的最高次数
[0116]
和其注解知k=(2m-l)n，k最 大值为20,则n、m的取值为n=l~20、m=l~10。计算机生成矩阵的步骤:先生成20X20的0 矩阵，再往里面填充元素，填充规则是在n X ( 2m-l )彡20的条件下，其元素
[0117] 2) jbsn:同样生成20X20的0矩阵，在nX(2m-l)彡20的条件下，
[0118] (2)计算 asn，bsn
[0119]给定的计算信号依然为：
[0120] u( ? t) = (cos ? t+0.6sin ? t) +
[0121] (0.lcos3 ? t+0.9sin3 ? t) + (0.2cos5 ? t+0.5sin5 ? t) +
[0122] (0?3cos7 ? t+0?3sin7 ? t) + (0?6cos9 ? t+0?2sin9 ? t)
[0123] 设分析长度为T0 = 2jt，步长St = T0/20/100(保证20次谐波的采样点为100)。计算u (cot)在每个采样点的值
[0124] u( ?t(k))(k = l，2,3 ? ? ? 2000)。
[0125] 周期矩形波San( co t)的matlab实现需要用到heaviside函数，具体实现代码如下：
[0126] San( co t) ：m=heaviside(wt)
[0127] for k = l :n
[0128] m=m
[0129] -2*heaviside(wt_(4*k_3)*pi/2/n)+2*heaviside(wt_
[0130] (4*k_l)*pi/2/n);
[0131] end
[0132] Sbn( co t) ： sbn = heaviside(wt)；
[0133] for k = l :n
[0134] sbn = sbn
[0135] -2*heaviside(wt-(2*k_l)*pi/n)+2*heaviside(wt_2*k*pi/n);
[0136] end
[0137] 计算每个采样点 San( 〇t)、Sbn( cot)的值 San( 〇t(k))、Sbn( c〇t(k))其中（k=l，2， 3 ? ? ? 2000) 〇
[0138] l)asn 的计算
[0139] 根据式(7)中知=J^w(w/)S"丨，先求解被积函数。将 u( ?t(k)#PSai( ? t(k))对应采样点的值相乘，其中k=l，2,3 ? ? ? 2000,令a(k)=u( ?t(k)) ? Sai( ?t(k))。 令h(k)=a(k) *St，即h(k)为如下图4阴影部分所示。再将h(k)从k=l，2,3* ? ?SOOO求和 即可算出&31的值，即&31 = 11(1)+11(2)+11(3)+...+11(2000)。同理可计算&；52，&；53，&； 54...。
[0140] 根据式（8)中=〔《，（如)6'",,|>/)6/(^)以及上述步骤，同样可计算出1331，1332， bS3...的值。
[0141] 2)bsn 的计算
[0142] 根据式（8)中，(⑴/)\,,(敁)^/(如）以及上述步骤，同样可计算出^1，匕2， bs3...的值。
[0143] (3)计算各次谐波正弦余弦分量系数
[0144] 有了（1)步中系数矩阵，第(2)步中asn、bsn，代入(10)中即可算得a n、bn。
[0147] 计算结果如下表2所示
[0148] 表2 20次谐波计算结果对比
[0149]
[0150] 在上述计算的基础上可以进一步实现谐波的实时检测。基于上述分析，谐波检测 次数N确定后，系数矩阵jasn、j bsn&已确定不会跟随信号的变换而变换，式（12)中唯一变化 的值是asn、b sn，即不断的跟随信号改变计算asn、bsn的实时数值即可实现谐波幅值的实时检 测 。
[0151] 由式(7)可得丨〇，(如/)5"丨（《/W(⑴/)，asl即为图5中曲线在0到2JI的积分。本方 法谐波实时检测的基本思路是在一个周期的计算基础上，继续往下采样。以20次谐波通用 计算【具体实施方式】中的采样方式为例:一个周期内采样点k = 2000,一个周期之后的第一个 采样点（k = 2001)为2ii+St，该处as冲被积函数的值为u( ?t(2001)) ? Sal( ?t(2001))， St ? u(?t(2001)) ? Sal(?t(2001))即为图中2ii+St处阴影部分面积。图中St处阴影部分面 积为St ? u(〇t(l)) ? SaKotQnoasi-St ? u(〇t(l)) ? Sai(〇t(l))+St ? u(〇t (2001)) ? Sal(c〇t(2001))即为被积函数在[St，23i+St]区间上的积分，令其为asl'。同理可 计算& 32'，&33'，&34'***。所：得的数组[& 31'，&32'，&33'，&34'***& 320']中包含了输入信号 在[St，23I+St]区间上的信息，将其带入式（12)即可求得输入信号在[St，23I+St]上的谐波余 弦分量系数信息。通过不断地采样反复用此方法计算不断变换的a sn即可得到不断变换的谐 波信息，实现谐波的实时检测。
[0152]同理可以计算正弦分量系数bn。
【主权项】
1. 一种基于方波傅里叶级数分解的新型谐波分析方法，其特征在于，包括以下步骤： 1) 输入待检信号：其中，AdPBi为所述基波 振幅，Aq和Bq为谐波振幅。t为相位； 所述待检信号含有基波和q个谐波，基波次数为m，谐波的次数分别为n2、n3……nq; 令系数 ki = m，i = l、2......q; 根据公式ki = (2mi_l )m，计算出系数:mi、m2......mq; 生成矩阵I: 首先，生成一个qXq的全0矩阵，此全0矩阵中，坐标为(m，m ? (2nu-l))的元素，若mX (2lHi- 1 ) ^ Ilmax j生成矩阵II: 首先，生成一个qXq的全0矩阵，此全0矩阵中，坐标为(m，m ? (2mi-l))的元素，若mX (2lHi- 1 ) ^ Ilmax2) 生成矩形波 Sanl( W t)、Sanq( W t)、 Sbnl( W t)和 Sbnq( W t):3) 米样： 设定步长St = 2Vnmax/200,在信号u( cot)上设定nmaxX200个采样点，每个采样点的值 为11(〇^)，€=1，2,3***]1111￡?\200;11 111￡?为次数最高的谐波的次数； 设定步长 St = 23l/nmax/200,在矩形波 Sanl( 〇 t)、Sanq( 〇 t)、Sbnl( W t)和 Sbnq( W t)上分别 设定nmaxX200个采样点，每个采样点的值为Sani( 〇tf)、Sanq( 〇tf)、Sbni( 〇tf)和Sbnq( 〇 tf)，f = 1，2，3 ? ? ? nmaxX 200 ;nmax为次数最高的谐波的次数； 4) 计算参数asni、bsni > asnq^FBbsnq ：将对应采样点u(〇tf)和Sanl(〇tf)的值相乘，得到a nl(f)，其中f=l，2,3* ? WOOX nmax，令hani(f)=ani(f) ? St，asni = hani(l)+hani(2)+hani(3)+ ? ? ? +hani(200Xnmax);将对应采样点u(〇tf)和Sbnl(〇tf)的值相乘，得到b nl(f)，其中f=l，2,3* ? WOOX nmax，令hbni(f)=bni(f) ? St，bsni = hbni(l)+hbni(2)+hbni(3)+? ? ?+hbni(200Xnmax);将对应采样点u(〇tf)和Sanq(〇tf)的值相乘，得到a nq(f)，其中f=l，2,3* ? WOOX Umax，令hanq(f )-&nq(f) * St?&snq - hanq( 1 )+hanq( 2 )+hanq( 3 ) + ? ? ? +hanq( 200 X Umax);将对应采样点u(〇tf)和Sbnq(〇tf)的值相乘，得到b nq(f)，其中f=l，2,3* ? WOOX nmax，令hbnq(f)=bnq(f) ? St，bsnq = hbnq(l)+hbnq(2)+hbnq(3)+ ? ? ? +hbnq(200Xnmax); 5) 分析得到基波和各次谐波正、余弦分量系数：其中，&1为基波余弦分量系数、aq为谐波余弦分量系数、匕为基波正弦分量系数、b q为谐 波正弦分量系数。2.根据权利要求1所述的一种基于方波傅里叶级数分解的新型谐波分析方法，其特征 在于，根据步骤5)获得的分析结果，进行谐波的实时检测。
【文档编号】G06F17/14GK106053940SQ201610647819
【公开日】2016年10月26日
【申请日】2016年8月9日 公开号201610647819.6, CN 106053940 A, CN 106053940A, CN 201610647819, CN-A-106053940, CN106053940 A, CN106053940A, CN201610647819, CN201610647819.6
【发明人】陆治国, 黎越, 王友, 罗邱银
【申请人】重庆大学