CN106053940A - 一种基于方波傅里叶级数分解的新型谐波分析方法 - Google Patents

Abstract

本发明公开一种基于方波傅里叶级数分解的新型谐波分析方法，利用计算机，通过程序，先输入待检测谐波分析次数N和待分析信号，根据待分析次数产生相应频率的矩形波周期函数，将各矩形波与待分析信号相乘所得的函数作为被积函数，再将其积分，由于矩形波的值只有1或者‑1，被积函数大大简化且在采样时相比正余玄函数能极大的节省时间。积分得到N个与各次谐波系数对应的参数，该N个参数构成的矩阵与系数矩阵相乘即得到各次谐波系数，其中系数矩阵由数学分析得到并由N决定且与待分析信号无关，一旦N确定，则系数矩阵不变。

Description

一种基于方波傅里叶级数分解的新型谐波分析方法

技术领域

本发明属于电力系统谐波分析建模技术领域，具体涉及电力系统中谐波系数实时计算建模方法。

背景技术

随着电子技术与现代工业的发展，电力系统中非线性负荷不断增加，电力电子器件的使用不断增加，使得电力系统谐波问题越来越严重，而且许多精密电子工业对电能质量要求偏高，故及时准确掌握电网中谐波状况对现代工业的发展已经对电力系统的稳定运行具有极大的意义。

电力系统中谐波来源通常来讲有三个方面。一是发电机绕组不是绝对对称，铁芯也不是绝对均匀，这导致在发电过程中产生少量谐波。二是来自电力变压器等输配电设备产生的谐波。三是来自于系统中的非线性负荷。

电力系统的谐波分析有多种算法，例如小波分析算法、神经网络算法、自适应谐波算法等，但通常采用快速傅里叶变换(FFT)算法。其中快速傅里叶虽然实现简单使用方便，但计算量大，实时性欠佳。小波分析法是近几年的热点领域，它是基于傅里叶变换的处理方法，其优点在于分析暂态谐波时实时性比较好，但小波基的选取难度较大。神经网络法应用于谐波检测的谐波源识别、谐波预测、谐波测量，但其依赖于算法的准确性与实现难易程度，有的算法还存在收敛速度及局部最小点等问题，这些都会对检测的精度及实时性产生影响。自适应谐波检测方法动态响应较慢。

发明内容

本发明的目的是针对现有谐波检测方法的不足，提供一种简便新颖且具有高实时性的谐波检测办法，并且该方法具有普遍适用性。

为实现本发明目的而采用的技术方案是这样的，一种基于方波傅里叶级数分解的新型谐波分析方法，其特征在于，包括以下步骤:

1)输入待检信号：

其中，A₁和B₁为基波振幅，A_q和B_q为谐波振幅；ωt为相位；

所述待检信号含有基波和q个谐波，基波次数为n₁，谐波的次数分别为n₂、n₃……n_q；

令系数k_i＝n_i，i＝1、2……q；

根据公式k_i＝(2m_i-1)n_i，计算出系数：m₁、m₂……m_q；

生成矩阵I：

首先，生成一个q×q的全0矩阵，此全0矩阵中，坐标为(n_i,n_i·(2m_i-1))的元素，若n_i×(2m_i-1)≤n_max，则填充为

生成矩阵II：

首先，生成一个q×q的全0矩阵，此全0矩阵中，坐标为(n_i,n_i·(2m_i-1))的元素，若n_i×(2m_i-1)≤n_max，则填充为

2)生成矩形波S_an1(ωt)、S_anq(ωt)、S_bn1(ωt)和S_bnq(ωt)：

$S_{{\mathrm{an}}_1}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{2}\\ -1& \frac{\mathrm{\&pi;}}{2}\&le;\mathrm{\&omega;}t<\frac{3\mathrm{\&pi;}}{2}\\ 1& \frac{3\mathrm{\&pi;}}{2}\&le;\mathrm{\&omega;}t\&le;2\mathrm{\&pi;}\end{array}\right.$

$S_{{\mathrm{an}}_q}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{2n_q}\\ -1& \frac{\mathrm{\&pi;}}{2n_q}\&le;\mathrm{\&omega;}t<\frac{3\mathrm{\&pi;}}{2n_q}\\ 1& \frac{3\mathrm{\&pi;}}{2n_q}\&le;\mathrm{\&omega;}t\&le;\frac{2\mathrm{\&pi;}}{n_q}\end{array}\right.,$

$S_{{\mathrm{bn}}_1}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\mathrm{\&pi;}\\ -1& \mathrm{\&pi;}\&le;\mathrm{\&omega;}t\&le;2\mathrm{\&pi;}\end{array}\right.$

$S_{{\mathrm{bn}}_q}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{n_q}\\ -1& \frac{\mathrm{\&pi;}}{n_q}\&le;\mathrm{\&omega;}t\&le;\frac{2\mathrm{\&pi;}}{n_q}\end{array}\right.;$

3)采样：

设定步长St＝2π/n_max/200，在信号u(ωt)上设定n_max×200个采样点，每个采样点的值为u(ωt_f)，f＝1,2,3···n_max×200；n_max为次数最高的谐波的次数；

设定步长St＝2π/n_max/200，在矩形波S_an1(ωt)、S_anq(ωt)、S_bn1(ωt)和S_bnq(ωt)上分别设定n_max×200个采样点，每个采样点的值为San₁(ωt_f)、San_q(ωt_f)、Sbn₁(ωt_f)和Sbn_q(ωt_f)，f＝1,2,3···n_max×200；n_max为次数最高的谐波的次数；

4)计算参数a_sn1、b_sn1、a_snq和b_snq：

先求解被积函数；

将对应采样点u(ωt_f)和S_an1(ωt_f)的值相乘，得到a_n1(f)，其中f＝1,2,3···200×n_max，令ha_n1(f)＝a_n1(f)·St，a_sn1＝ha_n1(1)+ha_n1(2)+ha_n1(3)+···+ha_n1(200×n_max)；

先求解被积函数；

将对应采样点u(ωt_f)和S_bn1(ωt_f)的值相乘，得到b_n1(f)，其中f＝1,2,3···200×n_max，令hb_n1(f)＝b_n1(f)·St，b_sn1＝hb_n1(1)+hb_n1(2)+hb_n1(3)+···+hb_n1(200×n_max)；

先求解被积函数；

将对应采样点u(ωt_f)和S_anq(ωt_f)的值相乘，得到a_nq(f)，其中f＝1,2,3···200×n_max，令ha_nq(f)＝a_nq(f)·St，a_snq＝ha_nq(1)+ha_nq(2)+ha_nq(3)+···+ha_nq(200×n_max)；

先求解被积函数；

将对应采样点u(ωt_f)和S_bnq(ωt_f)的值相乘，得到b_nq(f)，其中f＝1,2,3···200×n_max，令hb_nq(f)＝b_nq(f)·St，b_snq＝hb_nq(1)+hb_nq(2)+hb_nq(3)+···+hb_nq(200×n_max)；

5)分析得到基波和各次谐波正、余弦分量系数：

其中，a₁为基波余弦分量系数、a₂、a₃……a_q为谐波余弦分量系数、b₁为基波正弦分量系数、b₂、b₃……b_q为谐波正弦分量系数。

进一步，根据步骤5)获得的分析结果，进行谐波的实时检测。

值得说明的是，本发明可以利用计算机，通过程序，先输入待检测谐波分析次数N和待分析信号，将各矩形波与待分析信号相乘所得的函数作为被积函数，再将其积分，由于矩形波的值只有1或者-1，被积函数大大简化且在采样时相比正余玄函数能极大的节省时间。积分得到N个与各次谐波系数对应的参数，该N个参数构成的矩阵与系数矩阵相乘即得到各次谐波系数，其中系数矩阵由数学分析得到并由N决定且与待分析信号无关，一旦N确定，则系数矩阵不变。所述具体数学分析如下：

傅里叶级数的展开表达式为：

$u\left(\mathrm{\&omega;}t\right)=a_0+\underset{n=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}(a_n\cos n\mathrm{\&omega;}t+b_n\sin n\mathrm{\&omega;}t)$

式中：

$a_n=\frac{1}{\mathrm{\&pi;}}{\&Integral;}_0^{2\mathrm{\&pi;}}u\left(\mathrm{\&omega;}t\right)\cos n\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)$

$\begin{array}{c}{b}_n=\frac{1}{\mathrm{\&pi;}}{\&Integral;}_0^{2\mathrm{\&pi;}}u\left(\mathrm{\&omega;}t\right)\sin n\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)\\ n=1,2,3,\mathrm{...}\end{array}$

本发明是用特定的矩形波代替上述被积函数u(ωt)cos nωt中的余弦函数cos nωt。矩形波S_an如图1所示：

其数学表达式为：

$S_{an}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{2n}\\ -1& \frac{\mathrm{\&pi;}}{2n}\&le;\mathrm{\&omega;}t<\frac{3\mathrm{\&pi;}}{2n}\\ 1& \frac{3\mathrm{\&pi;}}{2n}\&le;\mathrm{\&omega;}t\&le;\frac{2\mathrm{\&pi;}}{n}\end{array}\right.---\left(1\right)$

现推导其傅里叶级数：

已知在区间[-l，l]上有：

$f\left(x\right)~\frac{a_0}{2}+\underset{n=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}(a_{n}cos\frac{n\mathrm{\&pi;}x}{l}+b_{n}sin\frac{n\mathrm{\&pi;}x}{l})$

$a_n=\frac{1}{l}{\&Integral;}_{-l}^{l}f\left(t\right)cos\frac{n\mathrm{\&pi;}t}{l}dt$

$b_n=\frac{1}{l}{\&Integral;}_{-l}^{l}f\left(t\right)sin\frac{n\mathrm{\&pi;}t}{l}dt$

依据该式子可以推导矩形波傅里叶级数：

$\begin{array}{c}{a}_{sbk}={\&Integral;}_0^{\frac{2\mathrm{\&pi;}}{n}}{S}_{an}\left(\mathrm{\&omega;}t\right)\sin \frac{k\mathrm{\&pi;}\mathrm{\&omega;}t}{\frac{\mathrm{\&pi;}}{n}}d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_0^{\frac{\mathrm{\&pi;}}{2n}}\sin kn\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)-{\&Integral;}_{\frac{\mathrm{\&pi;}}{2n}}^{\frac{3\mathrm{\&pi;}}{2n}}\sin kn\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)+{\&Integral;}_{\frac{3\mathrm{\&pi;}}{2n}}^{\frac{2\mathrm{\&pi;}}{n}}\sin kn\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)\\ =-\frac{1}{kn}\cos kn\mathrm{\&omega;}t{|}_0^{\frac{\mathrm{\&pi;}}{2n}}+\frac{1}{kn}\cos kn\mathrm{\&omega;}t{|}_{\frac{\mathrm{\&pi;}}{2n}}^{\frac{3\mathrm{\&pi;}}{2n}}-\frac{1}{kn}\cos kn\mathrm{\&omega;}t{|}_{\frac{3\mathrm{\&pi;}}{2n}}^{\frac{2\mathrm{\&pi;}}{n}}\\ =-\frac{2}{ku}\cos k\frac{\mathrm{\&pi;}}{2}+\frac{2}{kn}\cos k\frac{3\mathrm{\&pi;}}{2}\\ =0\end{array}$

所以可得到：

令图2矩形波为S_bn。

其表达式为：

$S_{bn}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{n}\\ -1& \frac{\mathrm{\&pi;}}{n}\&le;\mathrm{\&omega;}t\&le;\frac{2\mathrm{\&pi;}}{n}\end{array}\right.$

$\left(3\right)$

同理可推导其傅里叶级数：

$\begin{array}{c}{a}_{sbk}=\frac{1}{\frac{\mathrm{\&pi;}}{n}}{\&Integral;}_0^{\frac{2\mathrm{\&pi;}}{n}}{S}_{bn}\left(\mathrm{\&omega;}t\right)\cos \frac{k\mathrm{\&pi;}\mathrm{\&omega;}t}{\frac{\mathrm{\&pi;}}{n}}d\left(\mathrm{\&omega;}t\right)\\ =\frac{n}{\mathrm{\&pi;}}\&lsqb;{\&Integral;}_0^{\frac{\mathrm{\&pi;}}{n}}\cos kn\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)-{\&Integral;}_{\frac{\mathrm{\&pi;}}{n}}^{\frac{2\mathrm{\&pi;}}{n}}\cos kn\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)\&rsqb;\\ =\frac{n}{\mathrm{\&pi;}}\&lsqb;\frac{1}{kn}\sin kn\mathrm{\&omega;}t{|}_0^{\frac{\mathrm{\&pi;}}{n}}-\frac{1}{kn}\sin kn\mathrm{\&omega;}t{|}_{\frac{\mathrm{\&pi;}}{n}}^{\frac{2\mathrm{\&pi;}}{n}}\&rsqb;\\ =0\end{array}$

所以

现将谐波系数中的cosnωt替换为S_an的表达式，即式(2)，并做推导：

$\begin{array}{c}{a}_{s1}={\&Integral;}_0^{2\mathrm{\&pi;}}u\left(\mathrm{\&omega;}t\right)S_{a1}\left(\mathrm{\&omega;}t\right)d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_0^{2\mathrm{\&pi;}}\&lsqb;a_0+\underset{n=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}(a_n\cos n\mathrm{\&omega;}t+b_n\sin n\mathrm{\&omega;}t)\&rsqb;\&lsqb;\frac{4}{\mathrm{\&pi;}}\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{\cos (2m-1)\mathrm{\&omega;}t}{2m-1}\&rsqb;d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_2^{2\mathrm{\&pi;}}\frac{4}{\mathrm{\&pi;}}\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{\cos (2m-1)\mathrm{\&omega;}t}{2m-1}{a}_m\cos m\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)\\ =4\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{a_m}{2m-1}\end{array}---\left(5\right)$

同理将中sin nωt替换为S_bn：

$\begin{array}{c}{b}_{s1}={\&Integral;}_0^{2\mathrm{\&pi;}}u\left(\mathrm{\&omega;}t\right)S_{b1}\left(\mathrm{\&omega;}t\right)d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_0^{2\mathrm{\&pi;}}\&lsqb;a_0+\underset{n=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}(a_n\cos n\mathrm{\&omega;}t+b_n\sin n\mathrm{\&omega;}t)\&rsqb;\&lsqb;\frac{4}{\mathrm{\&pi;}}\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{\sin (2m-1)\mathrm{\&omega;}t}{2m-1}\&rsqb;d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_0^{2\mathrm{\&pi;}}\frac{4}{\mathrm{\&pi;}}\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{\sin (2m-1)\mathrm{\&omega;}t}{2m-1}{b}_m\sin m\mathrm{\&omega;}td\left(\mathrm{\&omega;}t\right)\\ =4\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{b_m}{2m-1}\end{array}---\left(6\right)$

由式(5)a_s1的表达式可以推导出a_sn的表达式：

$\begin{array}{c}{a}_{sn}={\&Integral;}_0^{2\mathrm{\&pi;}}u\left(\mathrm{\&omega;}t\right)S_{an}\left(\mathrm{\&omega;}t\right)d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_0^{2\mathrm{\&pi;}}\&lsqb;a_0+\underset{k=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}(a_k\cos k\mathrm{\&omega;}t+b_k\sin k\mathrm{\&omega;}t)\&rsqb;\&lsqb;-\frac{4}{\mathrm{\&pi;}}\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{1}{2m-1}{(-1)}^m\cos (2m-1)m\mathrm{\&omega;}t\&rsqb;d\left(\mathrm{\&omega;}t\right)\\ =-4\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{1}{2m-1}{(-1)}^m{a}_{(2m-1)n}\end{array}---\left(7\right)$

注：上述求和项中只有满足k＝(2m-1)n时才不为0。

同理可得b_sn表达式：

$\begin{array}{c}{b}_{sn}={\&Integral;}_0^{2\mathrm{\&pi;}}u\left(\mathrm{\&omega;}t\right)S_{bn}\left(\mathrm{\&omega;}t\right)d\left(\mathrm{\&omega;}t\right)\\ ={\&Integral;}_0^{2\mathrm{\&pi;}}\&lsqb;a_0+\underset{k=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}(a_k\cos k\mathrm{\&omega;}t+b_k\sin k\mathrm{\&omega;}t)\&rsqb;\&lsqb;\frac{4}{\mathrm{\&pi;}}\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{1}{2m-1}\sin (2m-1)n\mathrm{\&omega;}t\&rsqb;d\left(\mathrm{\&omega;}t\right)\\ =4\underset{m=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{1}{2m-1}{b}_{(2m-1)n}\end{array}---\left(8\right)$

观察式(7)可发现a_sn是a_n各项的求和(n＝1,2,3,···)。用矩阵表达如下：

其中矩阵I为系数矩阵，各个元素的值由式(7)中a_n的系数决定，且该矩阵与输入信号无关，只与需要计算的谐波次数N有关，N一旦确定，该系数矩阵则确定不变。同理由式(8)可求得b_n与b_sn相对应的系数矩阵矩阵II。对式(9)稍作改变得：

等式右边的可以用程序直接计算结果，则a_n的求解过程简化为两个数字矩阵的乘积，且式中S_an(ωt)为矩形波函数，其值只有1或者-1，使得被积函数大大简化，节省了计算量。流程图如图3。

本发明采用上述技术方案后，主要有以下效果：

1.本发明方法能够准确快速计算任意次谐波幅值，奇次偶次谐波能独立计算，能够很容易的提取某次谐波参数进行分析。

2.本发明方法具有通用性，只要提出需要检测的谐波次数，则可利用计算机程序算出系数矩阵I、II，并结合输入信号即可快速计算出谐波信息。方法简单，实用性强，便于推广应用。

3.可以高效快速实现输入信号谐波信息的实时检测。

本发明计算结果精确，计算速度快，效率高，可计算到任意高次谐波，广泛应用于电力系统中电能质量分析，为谐波治理打下可靠的基础。

附图说明

图1为矩形波S_an；

图2为矩形波S_bn；

图3为本发明的流程图；

图4为积分计算示意图；

图5为实时谐波检测示意图。

图中：。

具体实施方式

下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。

实施例1：

(1)输入待测函数

指定一个信号，只含基波，3、5、7、9次谐波，如下

u(ωt)＝(cosωt+0.6sinωt)+

(0.1cos3ωt+0.9sin3ωt)+(0.2cos5ωt+0.5sin5ωt)+

(0.3cos7ωt+0.3sin7ωt)+(0.6cos9ωt+0.2sin9ωt)

(11)

其中ωt＝[0,2π]，设定步长St＝2π/9/200，(注：最高为9次谐波，即步长设定为9次谐波周期除以200)，计算出输入信号在每个采样点的值。

根据式(7)，有k＝(2m-1)n，由于最高算到9次谐波，k的取值为k＝1，3，5，7，9，可得m、n的取值分别为n＝1，3，5，7，9；m＝1，2，3，4，5。根据式(7)，即可计算出系数矩阵如下：

$\begin{array}{cc}\begin{array}{c}n=1\\ 3\\ 5\\ 7\\ 9\end{array}& \begin{array}{c}m=\begin{array}{ccccc}1& 2& 3& 4& 5\end{array}\\ \left[\begin{array}{ccccc}-1& \frac{1}{3}& \frac{1}{5}& -\frac{1}{7}& -\frac{1}{9}\\ 0& -1& 0& 0& \frac{1}{3}\\ 0& 0& -1& 0& 0\\ 0& 0& 0& -1& 0\\ 0& 0& 0& 0& -1\end{array}\right]\end{array}\end{array}$

(2)周期矩形波的计算

根据第(1)步计算得出的n的取值，在Matlab中生成n＝1，3，5，7，9次的矩形波S_an(ωt)、S_bn(ωt)(注：n＝1时矩形波周期为2π，n＝3时周期为2π/3，依次类推)。同样根据步骤(1)中设定的步长，采样出5个不同周期矩形波在每个采样点的值，这些值均为1或-1。

(3)计算a_sn

第(2)步完成后，开始计算a_sn的值。根据式(7)中积分计算具体步骤如下：

以求a_sn为例，先求解被积函数。将u(ωt)和S_an(ωt)对应采样点的值相乘，令a(k)＝u(ωt(k))·S_an(ωt(k))(其中k为采样点编号，n＝1，3，5，7，9)。令h(k)＝a(k)·St，即h(k)为如下图4阴影部分所示。再将h(k)从k＝1,2,3···1800求和即可算出a_sn的值，即a_sn＝h(1)+h(2)+h(3)+···+h(1800)。a_sn和S_an(ωt(k))一一对应，从n＝1，3，5，7，9，S_a1带入计算求得a_s1，S_a2带入计算求得a_s2，依次类推。

(4)各次谐波余弦分量系数a_n的计算

第(3)步计算出a_sn的值后与第(1)步系数矩阵相乘即可得出各次谐波中正弦量与余弦量的系数。即：

$\left[\begin{array}{c}{a}_1\\ a_3\\ a_5\\ a_7\\ a_9\end{array}\right]=-\frac{1}{4}{\left[\begin{array}{ccccc}-1& \frac{1}{3}& \frac{1}{5}& -\frac{1}{7}& -\frac{1}{9}\\ 0& -1& 0& 0& \frac{1}{3}\\ 0& 0& -1& 0& 0\\ 0& 0& 0& -1& 0\\ 0& 0& 0& 0& -1\end{array}\right]}^{-1}\left[\begin{array}{c}{a}_{s1}\\ a_{s3}\\ a_{s5}\\ a_{s7}\\ a_{s9}\end{array}\right]$

用相同的步骤可以计算出b_sn及其对应的系数矩阵jbsn然后计算出b_n的值，计算结果如下表：

表1算例信号中各次谐波正余弦分量系数

实施例2：

本实施例涉及20次谐波

(1)输入需计算谐波的最高次数

1)jasn：N＝20，由式(7)和其注解知k＝(2m-1)n，k最大值为20，则n、m的取值为n＝1～20、m＝1～10。计算机生成矩阵的步骤：先生成20×20的0矩阵，再往里面填充元素，填充规则是在n×(2m-1)≤20的条件下，其元素

2)jbsn：同样生成20×20的0矩阵，在n×(2m-1)≤20的条件下，(依据式(8)中)。

(2)计算a_sn，b_sn

给定的计算信号依然为：

u(ωt)＝(cosωt+0.6sinωt)+

(0.1cos3ωt+0.9sin3ωt)+(0.2cos5ωt+0.5sin5ωt)+

(0.3cos7ωt+0.3sin7ωt)+(0.6cos9ωt+0.2sin9ωt)

设分析长度为T0＝2π，步长St＝T0/20/100(保证20次谐波的采样点为100)。计算u(ωt)在每个采样点的值

u(ωt(k))(k＝1,2,3···2000)。

周期矩形波S_an(ωt)的matlab实现需要用到heaviside函数，具体实现代码如下：

S_an(ωt)：m＝heaviside(wt)

for k＝1:n

m＝m

-2*heaviside(wt-(4*k-3)*pi/2/n)+2*heaviside(wt-

(4*k-1)*pi/2/n)；

end

S_bn(ωt)：sbn＝heaviside(wt)；

for k＝1:n

sbn＝sbn

-2*heaviside(wt-(2*k-1)*pi/n)+2*heaviside(wt-2*k*pi/n)；

end

计算每个采样点S_an(ωt)、S_bn(ωt)的值S_an(ωt(k))、S_bn(ωt(k))其中(k＝1,2,3···2000)。

1)a_sn的计算

根据式(7)中先求解被积函数。将u(ωt(k))和S_a1(ωt(k))对应采样点的值相乘，其中k＝1,2,3···2000，令a(k)＝u(ωt(k))·S_a1(ωt(k))。令h(k)＝a(k)·St，即h(k)为如下图4阴影部分所示。再将h(k)从k＝1,2,3···2000求和即可算出a_s1的值，即a_s1＝h(1)+h(2)+h(3)+···+h(2000)。同理可计算a_s2，a_s3，a_s4...。

根据式(8)中以及上述步骤，同样可计算出b_s1,b_s2,b_s3...的值。

2)b_sn的计算

根据式(8)中以及上述步骤，同样可计算出b_s1,b_s2,b_s3...的值。

(3)计算各次谐波正弦余弦分量系数

有了(1)步中系数矩阵，第(2)步中a_sn、b_sn，代入(10)中即可算得a_n、b_n。

$\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ .\\ .\\ .\\ a_{20}\end{array}\right]={\&lsqb;jasn\&rsqb;}^{-1}\left[\begin{array}{c}{a}_{s1}\\ a_{s2}\\ a_{s3}\\ .\\ .\\ .\\ a_{s20}\end{array}\right]$

$\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ .\\ .\\ .\\ b_{20}\end{array}\right]={\&lsqb;jbsn\&rsqb;}^{-1}\left[\begin{array}{c}{b}_{s1}\\ b_{s2}\\ b_{s3}\\ .\\ .\\ .\\ b_{s20}\end{array}\right]---\left(12\right)$

计算结果如下表2所示

表2 20次谐波计算结果对比

在上述计算的基础上可以进一步实现谐波的实时检测。基于上述分析，谐波检测次数N确定后，系数矩阵j_asn、j_bsn也已确定不会跟随信号的变换而变换，式(12)中唯一变化的值是a_sn、b_sn，即不断的跟随信号改变计算a_sn、b_sn的实时数值即可实现谐波幅值的实时检测。

由式(7)可得a_s1即为图5中曲线在0到2π的积分。本方法谐波实时检测的基本思路是在一个周期的计算基础上，继续往下采样。以20次谐波通用计算具体实施方式中的采样方式为例：一个周期内采样点k＝2000，一个周期之后的第一个采样点(k＝2001)为2π+St，该处a_s1中被积函数的值为u(ωt(2001))·S_a1(ωt(2001))，St·u(ωt(2001))·S_a1(ωt(2001))即为图中2π+St处阴影部分面积。图中St处阴影部分面积为St·u(ωt(1))·S_a1(ωt(1))。a_s1-St·u(ωt(1))·S_a1(ωt(1))+St·u(ωt(2001))·S_a1(ωt(2001))即为被积函数在[St,2π+St]区间上的积分，令其为a_s1'。同理可计算a_s2'，a_s3'，a_s4'···。所得的数组[a_s1',a_s2',a_s3',a_s4'···a_s20']中包含了输入信号在[St,2π+St]区间上的信息，将其带入式(12)即可求得输入信号在[St,2π+St]上的谐波余弦分量系数信息。通过不断地采样反复用此方法计算不断变换的a_sn即可得到不断变换的谐波信息，实现谐波的实时检测。

同理可以计算正弦分量系数b_n。

Claims (2)

1.一种基于方波傅里叶级数分解的新型谐波分析方法，其特征在于，包括以下步骤:

1)输入待检信号：

其中，A₁和B₁为所述基波振幅，A_q和B_q为谐波振幅。ωt为相位；

所述待检信号含有基波和q个谐波，基波次数为n₁，谐波的次数分别为n₂、n₃……n_q；

令系数k_i＝n_i，i＝1、2……q；

根据公式k_i＝(2m_i-1)n_i，计算出系数：m₁、m₂……m_q；

生成矩阵I：

首先，生成一个q×q的全0矩阵，此全0矩阵中，坐标为(n_i,n_i·(2m_i-1))的元素，若n_i×(2m_i-1)≤n_max，则填充为

生成矩阵II：

首先，生成一个q×q的全0矩阵，此全0矩阵中，坐标为(n_i,n_i·(2m_i-1))的元素，若n_i×(2m_i-1)≤n_max，则填充为

2)生成矩形波S_an1(ωt)、S_anq(ωt)、

S_bn1(ωt)和S_bnq(ωt)：

$S_{{\mathrm{an}}_1}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{2}\\ -1& \frac{\mathrm{\&pi;}}{2}\&le;\mathrm{\&omega;}t<\frac{3\mathrm{\&pi;}}{2}\\ 1& \frac{3\mathrm{\&pi;}}{2}\&le;\mathrm{\&omega;}t\&le;2\mathrm{\&pi;}\end{array}\right.$

$S_{{\mathrm{an}}_q}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{2n_q}\\ -1& \frac{\mathrm{\&pi;}}{2n_q}\&le;\mathrm{\&omega;}t<\frac{3\mathrm{\&pi;}}{2n_q}\\ 1& \frac{3\mathrm{\&pi;}}{2n_q}\&le;\mathrm{\&omega;}t\&le;\frac{2\mathrm{\&pi;}}{n_q}\end{array}\right.,$

$S_{{\mathrm{bn}}_1}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\mathrm{\&pi;}\\ -1& \mathrm{\&pi;}\&le;\mathrm{\&omega;}t\&le;2\mathrm{\&pi;}\end{array}\right.$

$S_{{\mathrm{bn}}_q}\left(\mathrm{\&omega;}t\right)=\left\{\begin{array}{cc}1& 0\&le;\mathrm{\&omega;}t<\frac{\mathrm{\&pi;}}{n_q}\\ -1& \frac{\mathrm{\&pi;}}{n_q}\&le;\mathrm{\&omega;}t\&le;\frac{2\mathrm{\&pi;}}{n_q}\end{array}\right.;$

3)采样：

设定步长St＝2π/n_max/200，在信号u(ωt)上设定n_max×200个采样点，每个采样点的值为u(ωt_f)，f＝1,2,3···n_max×200；n_max为次数最高的谐波的次数；

设定步长St＝2π/n_max/200，在矩形波S_an1(ωt)、S_anq(ωt)、S_bn1(ωt)和S_bnq(ωt)上分别设定n_max×200个采样点，每个采样点的值为San₁(ωt_f)、San_q(ωt_f)、Sbn₁(ωt_f)和Sbn_q(ωt_f)，f＝1,2,3···n_max×200；n_max为次数最高的谐波的次数；

4)计算参数a_sn1、b_sn1、a_snq和b_snq：

先求解被积函数；

将对应采样点u(ωt_f)和S_an1(ωt_f)的值相乘，得到a_n1(f)，其中f＝1,2,3···200×n_max，令ha_n1(f)＝a_n1(f)·St，a_sn1＝ha_n1(1)+ha_n1(2)+ha_n1(3)+···+ha_n1(200×n_max)；

先求解被积函数；

将对应采样点u(ωt_f)和S_bn1(ωt_f)的值相乘，得到b_n1(f)，其中f＝1,2,3···200×n_max，令hb_n1(f)＝b_n1(f)·St，b_sn1＝hb_n1(1)+hb_n1(2)+hb_n1(3)+···+hb_n1(200×n_max)；

先求解被积函数；

将对应采样点u(ωt_f)和S_anq(ωt_f)的值相乘，得到a_nq(f)，其中f＝1,2,3···200×n_max，令ha_nq(f)＝a_nq(f)·St，a_snq＝ha_nq(1)+ha_nq(2)+ha_nq(3)+···+ha_nq(200×n_max)；

先求解被积函数；

将对应采样点u(ωt_f)和S_bnq(ωt_f)的值相乘，得到b_nq(f)，其中f＝1,2,3···200×n_max，令hb_nq(f)＝b_nq(f)·St，b_snq＝hb_nq(1)+hb_nq(2)+hb_nq(3)+···+hb_nq(200×n_max)；

5)分析得到基波和各次谐波正、余弦分量系数：

其中，a₁为基波余弦分量系数、a_q为谐波余弦分量系数、b₁为基波正弦分量系数、b_q为谐波正弦分量系数。

2.根据权利要求1所述的一种基于方波傅里叶级数分解的新型谐波分析方法，其特征在于，根据步骤5)获得的分析结果，进行谐波的实时检测。