CN103399203A - High-precision harmonic parameter estimation method based on composite iterative algorithm - Google Patents

Abstract

The invention discloses a high-precision harmonic parameter estimation method based on a composite iterative algorithm. The method comprises the steps as follows: firstly, a harmonic amplitude value is directly calculated by utilizing a three-peak interpolation algorithm, and meanwhile, a harmonic phase is calculated to serve as an initial phase value for composite iteration, then the initial phase value is transmitted to a phase difference correction algorithm to obtain a frequency deviation value by utilizing the phase difference correction algorithm, and the phase difference correction algorithm is combined with a spectral leakage cancellation algorithm so as to obtain a phase value for the first-time iterative computation; circulating iteration is performed between the phase difference correction algorithm and the spectral leakage cancellation algorithm repeatedly until a phase iteration result meeting the error limit requirement is obtained. According to the invention, with the adoption of the three-peak interpolation algorithm, the estimation precision of the harmonic amplitude value is improved, and the phase obtained by utilizing the three-peak interpolation algorithm serves as an initial value, and the phase difference correction algorithm and the spectral leakage cancellation algorithm are combined to construct the composite iterative algorithm, so that the estimation precision of a harmonic phase and frequency is greatly improved.

Description

A kind of high precision of harmonic parameters based on complex iteration algorithm method of estimation

Technical field

The present invention relates to a kind of method of estimation of harmonic parameters for steady periodic signal, specifically refer to a kind of high precision of harmonic parameters based on complex iteration algorithm method of estimation, this method, for estimating harmonic amplitude, phase place and frequency parameter, belongs to power network signal frequency analysis technical field.

Background technology

Along with development and the industrial expansion of Power Electronic Technique, in electric system, nonlinear-load constantly increases.Nonlinear-load has also injected a large amount of higher hamonic waves in electrical network when bringing great economic benefit, strengthened the content of higher hamonic wave in the electrical network, has aggravated the distortion degree of electric signal.Therefore nonsinusoidal signal is carried out to frequency analysis electric energy metrical and power quality analysis and improvement are had to great Research Significance and practical value.

Fast Fourier Transform (FFT) (FFT) is in electric system, to carry out the most frequently used method of frequency analysis, but it exists spectrum leakage and fence effect, has affected the precision of frequency analysis, and window function and interpolation algorithm can address these problems to a certain extent.Unimodal interpolation and double peak interpolation algorithm utilize respectively a single spectral line and two spectral line information to estimate harmonic parameters, but these two kinds of algorithms are all ignored or part has been ignored long-range spectrum and revealed, and therefore utilize three peak interpolation algorithms of the Linear Combination Model structure of the maximum spectral line of amplitude and twice the large spectral lines in left and right thereof can further improve the estimated accuracy of harmonic amplitude parameter.

Although three peak interpolation algorithms estimate to have degree of precision to harmonic amplitude, the method only adopts the phase information of maximum spectral line to estimate harmonic phase, and precision is not high.

Summary of the invention

For the prior art above shortcomings, the purpose of this invention is to provide a kind of high precision of harmonic parameters based on complex iteration algorithm method of estimation, this method can improve the computational accuracy of frequency, amplitude and the phase place of power network signal greatly.

The present invention realizes that the technical solution of above-mentioned purpose is as follows:

A kind of high precision of harmonic parameters based on complex iteration algorithm method of estimation, at first utilize three peak interpolation algorithms directly to calculate harmonic amplitude, calculates simultaneously the phase place initial value of harmonic phase as complex iteration; Again the phase place initial value is delivered in the phase difference correction algorithm, utilizes the phase difference correction algorithm calculate to obtain the frequency departure amount, then leak and offset algorithm and combine with spectrum, obtain the phase value of iterative computation for the first time; So repeatedly in phase difference correction algorithm and spectrum, leak in offseting algorithm and carry out loop iteration, until finally obtain to meet the phase place iteration result that the limits of error requires.

Utilize three peak interpolation algorithms to calculate the harmonic amplitude step as follows:

1.1) signal x to be analyzed (t) is added to the hanning window, N+L the point (getting L=N/4) of sampling, top n point forms sequence 1, x _W1(n); Rear N point forms sequence 2, x _W2(n); Respectively two sequences are done to the FFT computing;

1.2) determine in two sequences peak value spectral line i under each harmonic frequency _lWith two large position of spectral line i _l+ 1, i _l-1, and obtain the amplitude X of its corresponding spectral line _W1(i _l-1), X _W1(i _l), X _W1(i _l+ 1); X _W2(i _l-1), X _W2(i _l), X _W2(i _l+ 1) and phase angle arg[X _W1(i _l-1)], arg[X _W1(i _l)], arg[X _W1(i _l+ 1)]; Arg[X _W2(i _l-1)], arg[X _W2(i _l)], arg[X _W2(i _l+ 1)];

1.3) make y ₁=X _W1(i _l-1), y ₂=X _W1(i _l), y ₃=X _W1(i _l+ 1), and by

Calculate β _l

1.4) by β _lThe substitution following formula calculates harmonic frequency deviation δ _l

${\mathrm{\δ}}_l=0.0011{\mathrm{\β}}_l^7-0.0079{\mathrm{\β}}_l^6+0.0181{\mathrm{\β}}_l^5\_0.0048{\mathrm{\β}}_l^4-0.0785{\mathrm{\β}}_l^3+0.0015{\mathrm{\β}}_l^2+0.6665{\mathrm{\β}}_l-0.5$

1.5) by δ _lThe substitution following formula calculates the amplitude A of each harmonic _l

$A_1=N(y_l+2y_2+y_3)(1.4726+0.5891{\mathrm{\&delta;}}_l+0.7221{\mathrm{\&delta;}}_l^2+0.288{\mathrm{\&delta;}}_l^3+0.2039{\mathrm{\&delta;}}_{\mathrm{<figure-callout\; id="4"\; label="l"\; filenames="HDA00003664945200011.png"\; state="\{\{state\}\}">l</figure-callout>}}^4+0.0895{\mathrm{\&delta;}}_l^5+0.049{\mathrm{\&delta;}}_l^6)$

Utilize the complex iteration algorithm to determine harmonic phase and frequency parameter, step is as follows:

2.1) by formula

Calculate the fundamental phase that obtains two sequences And as the initial value of complex iteration; Specification error limit ε, 0 → k;

2.2) by Calculate the fundamental frequency deviation, and by

Calculate other each harmonic frequency departures;

2.3) if 0≤δ _l<0.5, get μ=1; Otherwise get μ=-1, sequence 1 is calculated to the amplitude ratio

Sequence 2 is calculated to the amplitude ratio

2.4) by following formula, in conjunction with the amplitude that upper step obtains, compare the iterative value of calculating two sequence each harmonic phase angles again

2.5) if meet

Stop iteration; Otherwise,

K+1 → k, go to step 2.2);

2.6) phase angle of output each harmonic last iteration

And frequency departure

And by formula f _l=(i _l+ δ _l) △ f, l=1,2 ..., K, the Frequency Estimation of acquisition each harmonic is f as a result _l ^k.

Compared with prior art, the present invention has following beneficial effect:

1, the present invention has overcome the deficiency that in traditional frequency analysis, employing is unimodal or the double peak interpolation algorithm is lower to the harmonic amplitude estimated accuracy, by adopting three peak interpolation algorithms, has improved the estimated accuracy of harmonic amplitude.

2, utilize the phase difference correction method frequency departure to be had to degree of precision and compose the leakage opposition method and phase estimation is had to the characteristics of degree of precision, and utilize phase place that three peak interpolations obtain as initial value, both are combined and construct the complex iteration algorithm, greatly improved the estimated accuracy of harmonic phase and frequency departure.

3, the Simulation Example result of harmonic wave shows: the calculation of parameter precision of complex iteration algorithm is far above adding Hanning window method of interpolation, and computing time is suitable, thereby algorithm of the present invention has obvious advantage.

The accompanying drawing explanation

Fig. 1-harmonic parameters of the present invention is estimated process flow diagram.

Embodiment

The present invention utilizes three peak interpolations to have degree of precision to amplitude analysis, the phase difference correction algorithm has degree of precision to frequency analysis, spectrum is leaked opposition method and aspect phase analysis, is had the characteristics of better precision, three peak interpolation algorithms, phase difference correction algorithm and spectrum leakage are offseted to algorithm and combine, propose the algorithm of harmonics analysis of complex iteration.The superior function of utilizing three peak interpolations to calculate amplitude is directly asked for amplitude, calculate simultaneously phase place, phase result is delivered in the phase difference correction algorithm, utilize the phase difference correction Algorithm Analysis to obtain the frequency correction amount, with the spectrum leakage, offset algorithm again and combine, obtain the phase value of iterative computation for the first time.So repeatedly in phase difference correction algorithm and spectrum, leak in offseting algorithm and carry out loop iteration, until finally obtain the phase analysis result of degree of precision.

Below first introduce ultimate principle and the crest meter calculator body step of three peak interpolation algorithms, then introduce the complex iteration algorithm.

(1) three peak interpolation algorithm.

1) three peak interpolation method ultimate principles are as follows:

If power network signal is:

In formula, fk, Ak,

Represent respectively frequency, amplitude and phase place that the k subharmonic is corresponding, the higher harmonics number of times of K for asking for.

With cycle T _sThis signal is carried out to discrete sampling, obtains N point discrete series x (n):

ω in formula _k=2 π f _kT _s, T _sFor the sampling period.

Being provided with length is the discrete window function w (n) of N, and its spectrum expression formula is:

W(e ^jω)=W _o(ω)e ^-jCω (3)

In formula, W ₀(ω) be the real function of variable ω, C is the constant relevant to window function.

X (n) is carried out obtaining discrete signal x after windowing process _w(n):

x _w(n)=x(n)·w(n) (4)

Following formula is carried out to discrete Fourier transformation and can obtain x _w(n) frequency domain presentation is:

Non-integer-period sampled due in the actual samples process, make time window and between the signal period, do not meet the relation of integral multiple, establishes:

$\frac{{\mathrm{NT}}_s}{T_l}=i_l+{\mathrm{\&delta;}}_l---\left(6\right)$

In formula, i ₁For positive integer, δ ₁For frequency departure, T ₁For the signal period.

Any k subharmonic is had:

$\begin{array}{cc}\frac{{\mathrm{\&omega;}}_k}{\mathrm{\&Delta;\&omega;}}=\frac{2\mathrm{\&pi;k}{f}_1{T}_s}{2\mathrm{\&pi;}/N}=\frac{\mathrm{kN}{\mathrm{<figure-callout\; id="1"\; label="T\; s\; T"\; filenames="HDA00003664945200011.png"\; state="\{\{state\}\}">T</figure-callout>}}_s}{T_{}}=i_k+{\mathrm{\&delta;}}_k& k=\mathrm{1,2}...,K---\left(7\right)\end{array}$

In formula,

i _kBe position of spectral line corresponding to k subharmonic.

By (7), be easy to get:

ω _k=(i _k+δ _k)Δω (8)

For certain concrete harmonic wave k=l for example once, if the amplitude versus frequency characte W of window function ₀(ω) meet following formula:

$\begin{array}{ccc}{W}_0(i_l\mathrm{\&Delta;\&omega;}+{\mathrm{\&omega;}}_k)=0& k=\mathrm{1,2},...K& \\ W_0(i_l\mathrm{\&Delta;\&omega;}+{\mathrm{\&omega;}}_k)=0& k=\mathrm{1,2},...K& k\&NotEqual;l\end{array}---\left(9\right)$

Can think at ω=i this moment _l△ ω place, comprising first-harmonic is 0 in the positive and negative frequency component of other interior each harmonics and the negative frequency components of l subharmonic self, this spectral line can be regarded as and not affected by spectrum leakage.Composite type (5) and formula (8) can obtain:

I _lWith i _lThe amplitude expression formula of+1 spectral line is respectively:

$\left|X_w\left(i_l\right)\right|=\frac{A_l}{2}{W}_0(-{\mathrm{\&delta;}}_l\mathrm{\&Delta;\&omega;})---\left(11\right)$

$\left|X_w(i_l+1)\right|=\frac{A_l}{2}{W}_0(1-{\mathrm{\&delta;}}_l\mathrm{\&Delta;\&omega;})---\left(12\right)$

Order: ${\mathrm{\&beta;}}_l=\frac{y_3-y_1}{y_2}---\left(13\right)$

Wherein, y ₁, y ₂, y ₃Represent respectively successively X _w(i _l-1), X _w(i _l), X _w(i _l+ 1), namely harmonic wave frequency place comprises the amplitude of three nearest spectral lines of peak value spectral line.

Breadth of spectral line value expression shown in (11) formula is updated in (13) and has:

${\mathrm{\&beta;}}_l=\frac{|W_0((1-{\mathrm{\&delta;}}_l)2\mathrm{\&pi;}/N|-|W_0((-1-{\mathrm{\&delta;}}_l)2\mathrm{\&pi;}/N|}{|W_0((-{\mathrm{\&delta;}}_l)2\mathrm{\&pi;}/N|}---\left(14\right)$

Frequency departure δ while adopting polynomial fitting method to obtain to add Hanning window to following formula _lWith β _lBetween funtcional relationship as follows:

${\mathrm{\&delta;}}_1=0.0011{\mathrm{\&beta;}}_l^7-0.0079{\mathrm{\&beta;}}_l^6+0.0181{\mathrm{\&beta;}}_{l\_}^{5}0.0048{\mathrm{\&beta;}}_l^4-0.0785{\mathrm{\&beta;}}_l^3+0.0015{\mathrm{\&beta;}}_l^2+0.6665{\mathrm{\&beta;}}_1-0.5---\left(15\right)$

Thereby obtain comparatively accurate harmonic amplitude:

$A_l=\frac{2({\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="HDA00003664945200011.png"\; state="\{\{state\}\}">y</figure-callout>}}_1+2y_2+y_3)}{\left|W_0((-1-{\mathrm{\&delta;}}_l)2\mathrm{\&pi;}/N)\right|+2\left|W_0((-{\mathrm{\&delta;}}_l)2\mathrm{\&pi;}/N)\right|+|W_0((1-{\mathrm{\&delta;}}_1)2\mathrm{\&pi;}/N)}---\left(16\right)$

Adopt identical polynomial fitting method partly to carry out match to the denominator of formula (16), the harmonic amplitude estimation formulas while obtaining to add Hanning window is:

$A_1=N(y_l+2y_2+y_3)(1.4726+0.5891{\mathrm{\&delta;}}_l+0.7221{\mathrm{\&delta;}}_l^2+0.288{\mathrm{\&delta;}}_l^3+0.2039{\mathrm{\&delta;}}_{\mathrm{<figure-callout\; id="4"\; label="l"\; filenames="HDA00003664945200011.png"\; state="\{\{state\}\}">l</figure-callout>}}^4+0.0895{\mathrm{\&delta;}}_l^5+0.049{\mathrm{\&delta;}}_l^6)---\left(17\right)$

The computing formula of harmonic phase is:

2) three peak interpolations estimate that the implementation step of harmonic amplitude parameter is as follows:

The first step: choose sampled signal to be analyzed, it is added to Hanning window blocks and carry out Fast Fourier Transform (FFT);

Second step: find the peak value position of spectral line i under each harmonic frequency _lAnd the position i of these two large spectral lines in spectral line left and right _l+ 1 and i _l-1;

The 3rd step: the amplitude y that obtains these three spectral lines ₁=X _w(i _l-1), y ₂=X _w(i _l), y ₃=X _w(i _l+ 1), and according to formula (13) calculate β _l

The 4th step: by β _lSubstitution formula (15), thus frequency departure amount δ obtained _lValue;

The 5th step: by δ _lSubstitution formula (17) calculates the amplitude A that each harmonic is corresponding _l.

(2) the complex iteration algorithm for estimating of harmonic phase

Adopt three peak interpolation methods in step (1) to obtain harmonic amplitude and phase place, phase result is delivered in the phase difference correction algorithm, utilize the phase difference correction algorithm to obtain the frequency departure amount, then leak and offset algorithm and combine with spectrum, obtain the phase value of iterative computation for the first time.So repeatedly in phase difference correction algorithm and spectrum, leak in offseting algorithm and carry out loop iteration, until finally obtain the phase analysis result of degree of precision.

1) adopt the phase difference correction method to obtain frequency departure

Get time window T _w=mT, being located at sampling number in time window is N, sampling period T _s=mT/N, △ f=f _s/ N=1/NT _s.According to N+L point of time domain translation ratio juris sampling, the sequence of top n point is designated as to x ₁(n), a rear N point sequence is designated as x ₂(n).Sequence 2 is △ t=LT than 1 retardation time of sequence _s, the first-harmonic initial phase angle of sequence 2

First-harmonic initial phase angle with sequence 1

Between meet following relation:

In electrical network, the effect of regulating due to system makes frequency departure can not surpass 0.5Hz, i ₁The root spectral line is nearest apart from first-harmonic peak value frequency, and f is arranged this moment ₁=(i ₁+ δ ₁) △ f establishment, in substitution (19):

Thereby fundamental frequency deviation δ ₁For:

The frequency departure of other each harmonics is: δ _l=l δ ₁-round (l δ ₁) l=1,2 ..., K

The frequency parameter that therefore, can obtain each harmonic wave by formula (8) is:

f _l=(i _l+δ _l)△fl=1,2,…,K(21-2)

2) spectrum of harmonic phase is leaked and is offseted the high precision algorithm for estimating

Thereby according to the linear combination between many spectral lines, can eliminate long-range and compose the thought that the impact of leaking realizes the high accuracy analysis of harmonic amplitude, asking for of harmonic phase also can take identical mode to proofread and correct, thereby make up the impact of leaking due to the long-range spectrum when a single spectral line carries out the phase angle estimation, makes the defect that the angle values that obtains precision is lower.Considering the impact that long and short journey spectrum is leaked, and utilize maximum to be weighted and to offset with time large two amplitude spectral lines, is that the spectrum that realizes the analysis of high precision harmonic wave phase angle is leaked the core concept that offsets algorithm.

Cosine combination window function commonly used is:

$W\left(n\right)=\underset{h=0}{\overset{J-1}{\mathrm{\&Sigma;}}}{(-1)}^h{a}_h\cos \left(\frac{2\mathrm{\&pi;n}}{N}h\right)---\left(22\right)$

In formula, a _hFor the coefficient of cosine combination window function, J is the item number of cosine composite window.

To the l subharmonic, do not ignore the signal windowing expression formula that the long-range spectrum is leaked:

In formula, mainly comprise two parts, take minus sign as boundary.First half is negative frequency produces in frequency domain component, is called the departure vector, is designated as △ (i _l); Latter half is positive frequency produces in frequency domain component, is called short scope and leaks, and is designated as X (i _l) _s.

X _W(i _l)=X (i _l) _s± △ (i _l) (24)

Wherein: $\left|X{\left(i_l\right)}_s\right|=\frac{A_l{\mathrm{\&delta;}}_l\sin \left(\mathrm{\&pi;}{\mathrm{\&delta;}}_l\right)}{2\mathrm{\&pi;}}\underset{h=0}{\overset{J-1}{\mathrm{\&Sigma;}}}\frac{{(-1)}^h}{{\left({\mathrm{\&delta;}}_l\right)}^2-h^2}---\left(26\right)$

Suppose i _lAnd i _l+ μ root spectral line is respectively the maximum and inferior large spectral line of amplitude.As 0≤δ _l<0.5 o'clock, i _l+ 1 is amplitude time large spectral line, i.e. μ=1; As-0.5≤δ _l<0 o'clock, il-1 was amplitude time large spectral line, i.e. μ=-1.At this moment, the expression formula of inferior large spectral line is:

Wherein:

$|X_w(i_l+\mathrm{\&mu;}){|}_s=\frac{A_l(\mathrm{\&mu;}-{\mathrm{\&delta;}}_l)\sin \left(\mathrm{\&pi;}(\mathrm{\&mu;}-{\mathrm{\&delta;}}_l)\right)}{2\mathrm{\&pi;}}\underset{h=0}{\overset{J-1}{\mathrm{\&Sigma;}}}\frac{{(-1)}^h{a}_h}{{\mathrm{\&delta;}}_l^2-h^2}---\left(29\right)$

In phase place expression formula due to amplitude maximum, inferior large spectral line

With Value less, so work as X _w(i _l) while remaining unchanged, if the vectorial △ (i of departure _l) and short scope leakage rate X (i _l) _sBetween phase differential be 90 while spending,

Maximal value is arranged.This moment, angle deviation ratio maximum, inferior large spectral line can have following approximate:

The ratio that can be obtained maximum, inferior large spectral line amplitude deviation by formula (23) (24) is:

$\frac{\left|\mathrm{\&Delta;}\left(i_l\right)\right|}{\left|\mathrm{\&Delta;}(i_l+\mathrm{\&mu;})\right|}=\frac{(2i_l+{\mathrm{\&delta;}}_l)\underset{h=0}{\overset{J-1}{\mathrm{\&Sigma;}}}\frac{{(-1)}^h{a}_h}{{(2i_l+{\mathrm{\&delta;}}_l)}^2-h^2}}{(2i_l+{\mathrm{\&mu;}+\mathrm{\&delta;}}_l)\underset{h=0}{\overset{J-1}{\mathrm{\&Sigma;}}}\frac{{(-1)}^h{a}_h}{{(2i_l+\mathrm{\&mu;}+{\mathrm{\&delta;}}_l)}^2-h^2}}---\left(32\right)$

As the position of maximum spectral line i _lMuch larger than 1 o'clock, △ (i is arranged _l) ≈ △ (i _l+ μ) set up.Make X _w(i _l+ μ)=τ _h, X _w(i _l)=τ _s, formula (31) can be reduced to:

Convolution (23), (25) and (27) can obtain:

According to formula (33), formula (34) (35) is weighted on average, the phase place that can obtain harmonic wave is:

3) complex iteration method

As can be known by formula (21) and (36), frequency departure is relevant with the measuring accuracy of phase angle, and the precision of harmonic phase is relevant with frequency departure, therefore phase difference correction method and spectrum are leaked and offset the algorithm formation complex iteration algorithm that combines, can improve simultaneously the precision of frequency departure and phase angle measurement.

Complex iteration algorithm implementation step is as follows:

First step: choose N+L Window sampling signal, top n point forms sequence x ₁(n), rear N point forms sequence x ₂(n).

Second step: application three peak interpolation methods are analyzed these two sequences respectively, obtain the phase place of first-harmonic in these two sections sequences

With

0 → k.

Third step: with

With

As the phase place iterative initial value, in substitution formula (21), calculate the fundamental frequency deviation

And by

Calculate the frequency departure of l subharmonic.

The 4th step: two sequences are carried out to following identical operation: obtain harmonic frequency peak value breadth of spectral line value X in each sequence _w(i _l) and phase angle arg[X _w(i _l)] and the amplitude X of minor peaks spectral line _w(i _l+ μ) and phase angle arg[X _w(i _l+ μ)], and by $\frac{\left|X_w(i_l+\mathrm{\&mu;})\right|}{\left|X_w\left(i_l\right)\right|}=\frac{{\mathrm{\&tau;}}_h}{{\mathrm{\&tau;}}_s}$ Calculate the amplitude ratio;

The 5th step: calculate to obtain the each harmonic phase angle according to formula (36) With

The 6th step: establishing error precision is ε=10 ^-4If,

Stop iteration; Otherwise

K+1 → k, turn third step, until meet accuracy requirement.

The 7th step: the phase place of the last iteration of output And frequency departure

The 8th step: will

In substitution formula (21-2), obtain the Frequency Estimation result of each harmonic

Therefore, comprehensive above-mentioned introduction, the step that the present invention asks for harmonic parameters is as follows, asks simultaneously referring to Fig. 1:

1) at first utilize three peak interpolation algorithms to calculate harmonic amplitude, step is as follows:

1.1) signal x to be analyzed (t) is added to the hanning window, N+L the point of sampling, get L=N/4 usually, and top n point forms sequence 1, x _W1(n); Rear N point forms sequence 2, x _W2(n); Respectively two sequences are done to the FFT computing;

1.2) determine in two sequences peak value spectral line i under each harmonic frequency _lWith two large position of spectral line i _l+ 1, i _l-1, and obtain the amplitude X of its corresponding spectral line _W1(i _l-1), X _W1(i _l), X _W1(i _l+ 1); X _W2(i _l-1), X _W2(i _l), X _W2(i _l+ 1) and phase angle arg[X _W1(i _l-1)], arg[X _W1(i _l)], arg[X _W1(i _l+ 1)]; Arg[X _W2(i _l-1)], arg[X _W2(i _l)], arg[X _W2(i _l+ 1)];

1.3) make y ₁=X _W1(i _l-1), y ₂=X _W1(i _l), y ₃=X _W1(i _l+ 1), and by

Calculate β _l

1.4) by β _lThe substitution following formula calculates harmonic frequency deviation δ _l

${\mathrm{\&delta;}}_l=0.0011{\mathrm{\&beta;}}_l^7-0.0079{\mathrm{\&beta;}}_l^6+0.0181{\mathrm{\&beta;}}_l^5\_0.0048{\mathrm{\&beta;}}_l^4-0.0785{\mathrm{\&beta;}}_l^3+0.0015{\mathrm{\&beta;}}_l^2+0.6665{\mathrm{\&beta;}}_l-0.5$

1.5) by δ _lThe substitution following formula calculates the amplitude A of each harmonic _l

$A_1=N(y_l+2y_2+y_3)(1.4726+0.5891{\mathrm{\&delta;}}_l+0.7221{\mathrm{\&delta;}}_l^2+0.288{\mathrm{\&delta;}}_l^3+0.2039{\mathrm{\&delta;}}_{\mathrm{<figure-callout\; id="4"\; label="l"\; filenames="HDA00003664945200011.png"\; state="\{\{state\}\}">l</figure-callout>}}^4+0.0895{\mathrm{\&delta;}}_l^5+0.049{\mathrm{\&delta;}}_l^6)$

2) recycling complex iteration algorithm is determined harmonic phase and frequency parameter, and step is as follows:

2.1) by formula

Calculate the fundamental phase that obtains two sequences

And as the initial value of complex iteration; Specification error limit ε, the limits of error is set as ε=10 usually ^-4, 0 → k;

2.2) by Calculate the fundamental frequency deviation, and by

Calculate other each harmonic frequency departures;

2.3) if 0≤δ _l<0.5, get μ=1; Otherwise get μ=-1, sequence 1 is calculated to the amplitude ratio

Sequence 2 is calculated to the amplitude ratio

2.4) by following formula, in conjunction with the amplitude that upper step obtains, compare the iterative value of calculating two sequence each harmonic phase angles again

2.5) if meet Stop iteration; Otherwise,

Go to step 2.2);

2.6) phase angle of output each harmonic last iteration

And frequency departure

And by formula f _l=(i _l+ δ _l) △ f, l=1,2 ..., K, the Frequency Estimation result of acquisition each harmonic

The present invention is further elaborated below in conjunction with one, to implement example.

Getting harmonic wave power network signal model is:

$x\left(t\right)=100\sin (2\mathrm{\&pi;}{f}_{0}t+\frac{\mathrm{\&pi;}}{3})+2\sin (3\&times;\mathrm{\&pi;}{f}_{0}t+\frac{\mathrm{\&pi;}}{6})+10\sin (5\&times;2\mathrm{\&pi;}{f}_{0}t+\frac{\mathrm{\&pi;}}{2})+\sin (7\&times;2\mathrm{\&pi;}{f}_{0}t+\frac{\mathrm{\&pi;}}{2})$

Wherein, f ₀=49.8Hz, get sample frequency F _s=1kHz, sampling time 0.1s, sampling number N=100, L=25; To adding different window function and interpolation method, carry out simulation analysis, its amplitude relative error is as shown in table 1.

Table 1 amplitude relative error

Adopt respectively unimodal interpolation, double peak interpolation, three peak interpolations, phase difference correction method, spectrum leakage opposition method and complex iteration method to carry out frequency analysis to this power network signal, the error of establishing in the complex iteration algorithm is limited to ε=10 ^-4, its phase place relative error is as shown in table 2.

Table 2 phase place relative error

From Table 1 and Table 2, in frequency analysis, adopt three peak interpolations can greatly improve the harmonic amplitude accuracy of detection, and adopt the complex iteration algorithm, to the computational accuracy of harmonic phase, significantly raising is arranged.

Technical scheme of the present invention can be applicable to electric harmonic analysis, electric energy metrical and electric energy quality monitoring.

Finally explanation is, above embodiment is only unrestricted in order to technical scheme of the present invention to be described, although with reference to preferred embodiment, the present invention is had been described in detail, those of ordinary skill in the art is to be understood that, can modify or be equal to replacement technical scheme of the present invention, and not breaking away from aim and the scope of technical solution of the present invention, it all should be encompassed in the middle of claim scope of the present invention.

Claims (3)

1. the high precision of the harmonic parameters based on a complex iteration algorithm method of estimation, is characterized in that: at first utilize three peak interpolation algorithms directly to calculate harmonic amplitude, calculate simultaneously the phase place initial value of harmonic phase as complex iteration; Again the phase place initial value is delivered in the phase difference correction algorithm, utilizes the phase difference correction algorithm calculate to obtain the frequency departure amount, then leak and offset algorithm and combine with spectrum, obtain the phase value of iterative computation for the first time; So repeatedly in phase difference correction algorithm and spectrum, leak in offseting algorithm and carry out loop iteration, until finally obtain to meet the phase place iteration result that the limits of error requires.

2. harmonic parameters high precision method of estimation according to claim 1 is characterized in that: utilize three peak interpolation algorithms to calculate the harmonic amplitude steps as follows:

1.1) signal x to be analyzed (t) is added to the hanning window, N+L the point (getting L=N/4) of sampling, top n point forms sequence 1, x _W1(n); Rear N point forms sequence 2, x _W2(n); Respectively two sequences are done to the FFT computing;

1.2) determine in two sequences peak value spectral line i under each harmonic frequency _lWith two large position of spectral line i _l+ 1, i _l-1, and obtain the amplitude X of its corresponding spectral line _W1(i _l-1), X _W1(i _l), X _W1(i _l+ 1); X _W2(i _l-1), X _W2(i _l), X _W2(i _l+ 1) and phase angle arg[X _W1(i _l-1)], arg[X _W1(i _l)], arg[X _W1(i _l+ 1)]; Arg[X _W2(i _l-1)], arg[X _W2(i _l)], arg[X _W2(i _l+ 1)];

1.3) make y ₁=X _W1(i _l-1), y ₂=X _W1(i _l), y ₃=X _W1(i _l+ 1), and by Calculate β _l

1.4) by β _lThe substitution following formula calculates the harmonic frequency deviation

${\mathrm{\&delta;}}_l=0.0011{\mathrm{\&beta;}}_l^7-0.0079{\mathrm{\&beta;}}_l^6+0.0181{\mathrm{\&beta;}}_l^5\_0.0048{\mathrm{\&beta;}}_l^4-0.0785{\mathrm{\&beta;}}_l^3+0.0015{\mathrm{\&beta;}}_l^2+0.6665{\mathrm{\&beta;}}_l-0.5$

1.5) by δ _lThe substitution following formula calculates the amplitude A of each harmonic _l

$A_1=N(y_l+2y_2+y_3)(1.4726+0.5891{\mathrm{\&delta;}}_l+0.7221{\mathrm{\&delta;}}_l^2+0.288{\mathrm{\&delta;}}_l^3+0.2039{\mathrm{\&delta;}}_l^4+0.0895{\mathrm{\&delta;}}_l^5+0.049{\mathrm{\&delta;}}_l^6)$

Utilize the complex iteration algorithm to determine harmonic phase and frequency parameter, step is as follows:

2.1) by formula

Calculate the fundamental phase that obtains two sequences

And as the initial value of complex iteration; Specification error limit ε, 0 → k;

2.2) by

Calculate the fundamental frequency deviation, and by

Calculate other each harmonic frequency departures;

2.3) if 0≤δ _l<0.5, get μ=1; Otherwise get μ=-1, sequence 1 is calculated to the amplitude ratio

Sequence 2 is calculated to the amplitude ratio

2.4) by following formula, in conjunction with the amplitude that upper step obtains, compare the iterative value of calculating two sequence each harmonic phase angles again

2.5) if meet

Stop iteration; Otherwise, Go to step 2.2);

2.6) phase angle of output each harmonic last iteration

And frequency departure

And by formula f _l=(i _l+ δ _l) △ f, l=1,2 ..., K, the Frequency Estimation result of acquisition each harmonic

3. harmonic parameters high precision method of estimation according to claim 1 and 2, it is characterized in that: the described limits of error is set as ε=10 ^-4.

Images