CN103399203B - A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm - Google Patents

A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm Download PDF

Info

Publication number
CN103399203B
CN103399203B CN201310354109.0A CN201310354109A CN103399203B CN 103399203 B CN103399203 B CN 103399203B CN 201310354109 A CN201310354109 A CN 201310354109A CN 103399203 B CN103399203 B CN 103399203B
Authority
CN
China
Prior art keywords
harmonic
phase
algorithm
delta
calculate
Prior art date
Application number
CN201310354109.0A
Other languages
Chinese (zh)
Other versions
CN103399203A (en
Inventor
李春燕
魏蔚
张淮清
许中
付志红
张谦
Original Assignee
重庆大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 重庆大学 filed Critical 重庆大学
Priority to CN201310354109.0A priority Critical patent/CN103399203B/en
Publication of CN103399203A publication Critical patent/CN103399203A/en
Application granted granted Critical
Publication of CN103399203B publication Critical patent/CN103399203B/en

Links

Abstract

The invention discloses a kind of High-precision harmonic parameter estimation method based on composite iterative algorithm, first utilize three peak interpolation algorithms directly to calculate harmonic amplitude, calculate the phase place initial value of harmonic phase as complex iteration simultaneously; Again phase place initial value is delivered in phase difference correction algorithm, utilizes phase difference correction algorithm to obtain frequency departure amount, then combine with spectrum feed-through nulling algorithm, obtain the phase value of first time iterative computation; So repeatedly in phase difference correction algorithm and spectrum feed-through nulling algorithm, carry out loop iteration, until obtain the phase place iteration result meeting the limits of error and require.The present invention adopts three peak interpolation algorithms to improve the estimated accuracy of harmonic amplitude, utilize the phase place of three peak interpolations acquisitions as initial value, combine phase difference correction method and spectrum feed-through nulling method structure composite iterative algorithm, substantially increases harmonic phase and frequency estimation accuracy.

Description

A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm

Technical field

The present invention relates to a kind of harmonic parameters method of estimation for steady periodic signal, a kind of High-precision harmonic parameter estimation method based on composite iterative algorithm of concrete finger, 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 is also filled with a large amount of higher hamonic waves while bringing great economic benefit in electrical network, increases the content of higher hamonic wave in electrical network, exacerbates the distortion degree of electric signal.Therefore frequency analysis is carried out to nonsinusoidal signal, to electric energy metrical and power quality analysis and improvement, there is great Research Significance and practical value.

The most frequently used method of frequency analysis is carried out in Fast Fourier Transform (FFT) (FFT) in electric system, but it exists spectrum leakage and fence effect, have impact on 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 a single spectral line and two spectral line information to estimate harmonic parameters respectively, but these two kinds of algorithms are all ignored or part have ignored the leakage of long-range spectrum, therefore utilize three peak interpolation algorithms of the Linear Combination Model structure of amplitude maximum spectral line and two, left and right time large spectral line thereof can improve the estimated accuracy of harmonic amplitude parameter further.

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 prior art above shortcomings, the object of this invention is to provide a kind of High-precision harmonic parameter estimation method based on composite iterative algorithm, this method can improve the computational accuracy of the frequency of power network signal, amplitude and phase place greatly.

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

Based on a High-precision harmonic parameter estimation method for composite iterative algorithm, first utilize three peak interpolation algorithms directly to calculate harmonic amplitude, calculate the phase place initial value of harmonic phase as complex iteration simultaneously; Again phase place initial value is delivered in phase difference correction algorithm, utilizes phase difference correction algorithm to calculate and obtain frequency departure amount, then combine with spectrum feed-through nulling algorithm, obtain the phase value of first time iterative computation; So repeatedly in phase difference correction algorithm and spectrum feed-through nulling algorithm, carry out loop iteration, obtain the phase place iteration result meeting the limits of error and require until final.

Three peak interpolation algorithms are utilized to calculate harmonic amplitude step as follows:

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

1.2) peak value spectral line i under each harmonic frequency is determined in two sequences lwith two large position of spectral line i l+ 1, i l-1, and the amplitude X obtaining 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) y is made 1=X w1(i l-1), y 2=X w1(i l), y 3=X w1(i l+ 1), and by calculate β l;

1.4) by β lsubstitute into following formula and calculate harmonic frequency deviation δ l;

δ l = 0.0011 β l 7 - 0.0079 β l 6 + 0.0181 β l 5 _ 0.0048 β l 4 - 0.0785 β l 3 + 0.0015 β l 2 + 0.6665 β l - 0.5

1.5) by δ lsubstitute into the amplitude A that following formula calculates each harmonic l;

A 1 = N ( y l + 2 y 2 + y 3 ) ( 1.4726 + 0.5891 δ l + 0.7221 δ l 2 + 0.288 δ l 3 + 0.2039 δ l 4 + 0.0895 δ l 5 + 0.049 δ l 6 )

Utilize composite iterative algorithm determination harmonic phase and frequency parameter, step is as follows:

2.1) by formula calculate the fundamental phase of acquisition two sequences and as the initial value of complex iteration; Specification error limit ε, 0 → k;

2.2) by calculate fundamental frequency deviation, and by calculate other each harmonic frequency departures;

2.3) if 0≤δ l<0.5, then get μ=1; Otherwise get μ=-1, Amplitude Ration is calculated to sequence 1 amplitude Ration is calculated to sequence 2

2.4) Amplitude Ration obtained in conjunction with upper step again by following formula calculates the iterative value of two sequence each harmonic phase angles

2.5) if meet then stop iteration; Otherwise, k+1 → k, goes to step 2.2);

2.6) phase angle of the last iteration of each harmonic is exported and frequency departure and by formula f l=(i l+ δ l) △ f, l=1,2 ..., K, obtains the Frequency Estimation result f of each harmonic l k.

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

1, instant invention overcomes in conventional harmonic analysis and adopt unimodal or that double peak interpolation algorithm is lower to harmonic amplitude estimated accuracy deficiency, by the estimated accuracy adopting three peak interpolation algorithms to improve harmonic amplitude.

2, utilize phase difference correction method to have degree of precision to frequency departure and compose feed-through nulling method has degree of precision feature to phase estimation, and utilize the phase place of three peak interpolations acquisitions as initial value, combine both structure composite iterative algorithm, substantially increases the estimated accuracy of harmonic phase and frequency departure.

3, the Simulation Example result of harmonic wave shows: the parameter computational accuracy of composite iterative algorithm is far above adding Hanning window method of interpolation, and computing time is suitable, and thus, algorithm of the present invention has obvious advantage.

Accompanying drawing explanation

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

Embodiment

The present invention utilizes three peak interpolations to have degree of precision to amplitude analysis, phase difference correction algorithm has degree of precision to frequency analysis, spectrum feed-through nulling method has the feature of better precision in phase analysis, three peak interpolation algorithms, phase difference correction algorithm and spectrum feed-through nulling algorithm are combined, proposes the algorithm of harmonics analysis of complex iteration.Three peak interpolations are utilized directly to ask for amplitude to the superior function that amplitude calculates, calculate phase place simultaneously, phase result is delivered in phase difference correction algorithm, phase difference correction Algorithm Analysis is utilized to obtain frequency correction amount, combine with spectrum feed-through nulling algorithm again, obtain the phase value of first time iterative computation.So repeatedly in phase difference correction algorithm and spectrum feed-through nulling algorithm, carry out loop iteration, until the final phase analysis result obtaining degree of precision.

Below ultimate principle and the crest meter of first introducing three peak interpolation algorithms calculate concrete steps, then introduce composite iterative 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 frequency, amplitude and phase place that kth subharmonic is corresponding respectively, K for required by the most higher harmonics number of times got.

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

ω in formula k=2 π f kt s, T sfor the sampling period.

Be provided with discrete window function w (n) that length is N, its spectrum expression formula is:

W(e )=W o(ω)e -jCω(3)

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

Discrete signal x is obtained after windowing process is carried out to x (n) w(n):

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

Discrete Fourier transformation is carried out to above formula and can obtain x wn the frequency domain presentation of () is:

Non-integer-period sampled due in actual samples process, makes the relation not meeting integral multiple between time window and signal period, if:

NT s T l = i l + &delta; l - - - ( 6 )

In formula, i 1for positive integer, δ 1for frequency departure, T 1for the signal period.

Have any kth subharmonic:

&omega; k &Delta;&omega; = 2 &pi;k f 1 T s 2 &pi; / N = kN T s T 1 = i k + &delta; k k = 1,2 . . . , K - - - ( 7 )

In formula, i kfor the position of spectral line that kth subharmonic is corresponding.

Be easy to get by (7):

ω k=(i kk)Δω (8)

For certain once concrete harmonic wave such as k=l, if the amplitude versus frequency characte W of window function 0(ω) following formula is met:

W 0 ( i l &Delta;&omega; + &omega; k ) = 0 k = 1,2 , . . . K W 0 ( i l &Delta;&omega; + &omega; k ) = 0 k = 1,2 , . . . K k &NotEqual; l - - - ( 9 )

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

Then i-th lwith i-th lthe amplitude expression of+1 spectral line is respectively:

| X w ( i l ) | = A l 2 W 0 ( - &delta; l &Delta;&omega; ) - - - ( 11 )

| X w ( i l + 1 ) | = A l 2 W 0 ( 1 - &delta; l &Delta;&omega; ) - - - ( 12 )

Order: &beta; l = y 3 - y 1 y 2 - - - ( 13 )

Wherein, y 1, y 2, y 3represent X respectively successively w(i l-1), X w(i l), X w(i l+ 1), namely harmonic wave frequency place comprises the amplitude of the nearest spectral line of three of peak value spectral line.

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

&beta; l = | W 0 ( ( 1 - &delta; l ) 2 &pi; / N | - | W 0 ( ( - 1 - &delta; l ) 2 &pi; / N | | W 0 ( ( - &delta; l ) 2 &pi; / N | - - - ( 14 )

Frequency departure δ when adopting polynomial fitting method to obtain to add Hanning window to above formula lwith β lbetween funtcional relationship as follows:

&delta; 1 = 0.0011 &beta; l 7 - 0.0079 &beta; l 6 + 0.0181 &beta; l _ 5 0.0048 &beta; l 4 - 0.0785 &beta; l 3 + 0.0015 &beta; l 2 + 0.6665 &beta; 1 - 0.5 - - - ( 15 )

Thus obtain comparatively accurate harmonic amplitude:

A l = 2 ( y 1 + 2 y 2 + y 3 ) | W 0 ( ( - 1 - &delta; l ) 2 &pi; / N ) | + 2 | W 0 ( ( - &delta; l ) 2 &pi; / N ) | + | W 0 ( ( 1 - &delta; 1 ) 2 &pi; / N ) - - - ( 16 )

Adopt the denominator part of identical polynomial fitting method to formula (16) to carry out matching, the harmonic amplitude estimation formulas obtained when adding Hanning window is:

A 1 = N ( y l + 2 y 2 + y 3 ) ( 1.4726 + 0.5891 &delta; l + 0.7221 &delta; l 2 + 0.288 &delta; l 3 + 0.2039 &delta; l 4 + 0.0895 &delta; l 5 + 0.049 &delta; l 6 ) - - - ( 17 )

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, adds Hanning window to it and 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 two large spectral lines about this spectral line l+ 1 and i l-1;

3rd step: the amplitude y obtaining these three spectral lines 1=X w(i l-1), y 2=X w(i l), y 3=X w(i l+ 1), and according to formula (13) calculate β l;

4th step: by β lsubstitution formula (15), thus obtain frequency departure amount δ lvalue;

5th step: by δ lsubstitution formula (17) calculates amplitude A corresponding to each harmonic l.

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

Three peak interpolation methods in step (1) are adopted to obtain harmonic amplitude and phase place, phase result is delivered in phase difference correction algorithm, utilize phase difference correction algorithm to obtain frequency departure amount, then combine with spectrum feed-through nulling algorithm, obtain the phase value of first time iterative computation.So repeatedly in phase difference correction algorithm and spectrum feed-through nulling algorithm, carry out loop iteration, until the final phase analysis result obtaining degree of precision.

1) phase difference correction method is adopted to obtain frequency departure

Get time window T w=mT, being located at sampling number in time window is N, then sampling period T s=mT/N, △ f=f s/ N=1/NT s.According to time domain translation ratio juris sampling N+L point, the sequence of top n point is designated as x 1n (), rear N number of point sequence is designated as x 2(n).Sequence 2 is △ t=LT than sequence 1 retardation time s, the first-harmonic initial phase angle of sequence 2 with the first-harmonic initial phase angle of sequence 1 between meet following relation:

In electrical network, the effect due to system fading margin makes frequency departure can not more than 0.5Hz, and i-th 1root spectral line distance first-harmonic peak value frequency is nearest, now has f 1=(i 1+ δ 1) △ f sets up, to substitute in (19):

Thus fundamental frequency deviation δ 1for:

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

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

f l=(i ll)△fl=1,2,…,K(21-2)

2) the spectrum feed-through nulling high precision algorithm for estimating of harmonic phase

Can eliminate according to the linear combination between many spectral lines impact that long-range spectrum leaks thus realize the thought of the high accuracy analysis of harmonic amplitude, asking for of harmonic phase also can be taked identical mode to correct, thus makes up the impact of leaking due to long-range spectrum when a single spectral line carries out phase angular estimation and make the defect that the angle values precision of acquisition is lower.Considering the impact that long and short journey spectrum is leaked, and utilize maximum and secondary large two amplitude spectral lines to be weighted to offset, is the core concept realizing the spectrum feed-through nulling algorithm that high precision harmonic wave phase angle is analyzed.

Conventional cosine combination window function is:

W ( n ) = &Sigma; h = 0 J - 1 ( - 1 ) h a h cos ( 2 &pi;n N h ) - - - ( 22 )

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

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

Mainly comprising two parts in formula, take minus sign as boundary.First half is the component that negative frequency produces in frequency domain, is called departure vector, is designated as △ (i l); Latter half is the component that positive frequency produces in frequency domain, is called that short scope is leaked, is designated as X (i l) s.

Then X w(i l)=X (i l) s± △ (i l) (24)

Wherein: | X ( i l ) s | = A l &delta; l sin ( &pi; &delta; l ) 2 &pi; &Sigma; h = 0 J - 1 ( - 1 ) h ( &delta; l ) 2 - h 2 - - - ( 26 )

Suppose i-th land i l+ μ root spectral line is respectively amplitude maximum and time large spectral line.As 0≤δ lduring <0.5, then i l+ 1 is amplitude time large spectral line, i.e. μ=1; As-0.5≤δ lduring <0, then il-1 is amplitude time large spectral line, i.e. μ=-1.Now, the expression formula of secondary large spectral line is:

Wherein:

| X w ( i l + &mu; ) | s = A l ( &mu; - &delta; l ) sin ( &pi; ( &mu; - &delta; l ) ) 2 &pi; &Sigma; h = 0 J - 1 ( - 1 ) h a h &delta; l 2 - h 2 - - - ( 29 )

In phase place expression formula due to amplitude maximum, secondary large spectral line with value less, so work as X w(i l) when remaining unchanged, if the vectorial △ (i of departure l) and short scope leakage rate X (i l) sbetween phase differential when being 90 degree, there is maximal value.Now the angle deviation ratio of maximum, secondary large spectral line can have with lower aprons:

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

| &Delta; ( i l ) | | &Delta; ( i l + &mu; ) | = ( 2 i l + &delta; l ) &Sigma; h = 0 J - 1 ( - 1 ) h a h ( 2 i l + &delta; l ) 2 - h 2 ( 2 i l + &mu; + &delta; l ) &Sigma; h = 0 J - 1 ( - 1 ) h a h ( 2 i l + &mu; + &delta; l ) 2 - h 2 - - - ( 32 )

As the position i of maximum spectral line lmuch larger than 1 time, have △ (i l) ≈ △ (i l+ μ) set up.Make X w(i l+ μ)=τ h, X w(i l)=τ s, then formula (31) can be reduced to:

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

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

3) complex iteration method

From 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 combine phase difference correction method with spectrum feed-through nulling algorithm formation composite iterative algorithm, can improve the precision of frequency departure and phase angle measurement simultaneously.

Composite iterative algorithm implementation step is as follows:

First step: choose N+L Window sampling signal, top n point forms sequence x 1(n), rear N number of formation sequence x 2(n).

Second step: three peak interpolation methods of applying are analyzed these two sequences respectively, obtain the phase place of first-harmonic in these two sections of sequences with 0 → k.

Third step: with with as phase place iterative initial value, substitute in formula (21) and calculate fundamental frequency deviation and by calculate the frequency departure of l subharmonic.

4th step: operation identical is as follows carried out to two sequences: obtain harmonic frequency peaks spectral line amplitude 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 | X w ( i l + &mu; ) | | X w ( i l ) | = &tau; h &tau; s Calculate Amplitude Ration;

5th step: calculate according to formula (36) and to obtain each harmonic phase angle with

6th step: set error precision as ε=10 -4if, then stop iteration; Otherwise k+1 → k, turns third step, until meet accuracy requirement.

7th step: the phase place exporting last iteration and frequency departure

8th step: will the Frequency Estimation result of each harmonic is obtained in substitution formula (21-2)

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

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

1.1) signal x (t) to be analyzed is added hanning window, N+L point of sampling, gets L=N/4 usually, and top n point forms sequence 1, x w1(n); Rear N number of formation sequence 2, x w2(n); Respectively FFT computing is done to two sequences;

1.2) peak value spectral line i under each harmonic frequency is determined in two sequences lwith two large position of spectral line i l+ 1, i l-1, and the amplitude X obtaining 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) y is made 1=X w1(i l-1), y 2=X w1(i l), y 3=X w1(i l+ 1), and by calculate β l;

1.4) by β lsubstitute into following formula and calculate harmonic frequency deviation δ l;

&delta; l = 0.0011 &beta; l 7 - 0.0079 &beta; l 6 + 0.0181 &beta; l 5 _ 0.0048 &beta; l 4 - 0.0785 &beta; l 3 + 0.0015 &beta; l 2 + 0.6665 &beta; l - 0.5

1.5) by δ lsubstitute into the amplitude A that following formula calculates each harmonic l;

A 1 = N ( y l + 2 y 2 + y 3 ) ( 1.4726 + 0.5891 &delta; l + 0.7221 &delta; l 2 + 0.288 &delta; l 3 + 0.2039 &delta; l 4 + 0.0895 &delta; l 5 + 0.049 &delta; l 6 )

2) recycle composite iterative algorithm determination harmonic phase and frequency parameter, step is as follows:

2.1) by formula calculate the fundamental phase of acquisition two sequences and as the initial value of complex iteration; Specification error limit ε, the usual limits of error is set as ε=10 -4, 0 → k;

2.2) by calculate fundamental frequency deviation, and by calculate other each harmonic frequency departures;

2.3) if 0≤δ l<0.5, then get μ=1; Otherwise get μ=-1, Amplitude Ration is calculated to sequence 1 amplitude Ration is calculated to sequence 2

2.4) Amplitude Ration obtained in conjunction with upper step again by following formula calculates the iterative value of two sequence each harmonic phase angles

2.5) if meet then stop iteration; Otherwise, go to step 2.2);

2.6) phase angle of the last iteration of each harmonic is exported and frequency departure and by formula f l=(i l+ δ l) △ f, l=1,2 ..., K, obtains the Frequency Estimation result of each harmonic

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

Getting harmonic signal model is:

x ( t ) = 100 sin ( 2 &pi; f 0 t + &pi; 3 ) + 2 sin ( 3 &times; &pi; f 0 t + &pi; 6 ) + 10 sin ( 5 &times; 2 &pi; f 0 t + &pi; 2 ) + sin ( 7 &times; 2 &pi; f 0 t + &pi; 2 )

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

Table 1 amplitude relative error

Adopt unimodal interpolation, double peak interpolation, three peak interpolations, phase difference correction method respectively, spectrum reveals opposition method and complex iteration method carries out frequency analysis to this power network signal, if the error in composite iterative 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 harmonic amplitude accuracy of detection, and adopt the computational accuracy of composite iterative algorithm to harmonic phase to have significantly to improve.

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

What finally illustrate is, above embodiment is only in order to illustrate technical scheme of the present invention and unrestricted, although with reference to preferred embodiment to invention has been detailed description, those of ordinary skill in the art is to be understood that, can modify to technical scheme of the present invention or equivalent replacement, and not departing from aim and the scope of technical solution of the present invention, it all should be encompassed in the middle of right of the present invention.

Claims (3)

1. based on a High-precision harmonic parameter estimation method for composite iterative algorithm, it is characterized in that: first utilize three peak interpolation algorithms directly to calculate harmonic amplitude, calculate the phase place initial value of harmonic phase as complex iteration simultaneously; Again phase place initial value is delivered in phase difference correction algorithm, utilizes phase difference correction algorithm to calculate and obtain frequency departure amount, then combine with spectrum feed-through nulling algorithm, obtain the phase value of first time iterative computation; So repeatedly in phase difference correction algorithm and spectrum feed-through nulling algorithm, carry out loop iteration, obtain the phase place iteration result meeting the limits of error and require until final.
2. High-precision harmonic parameter estimation method according to claim 1, is characterized in that: utilize three peak interpolation algorithms to calculate harmonic amplitude step as follows:
1.1) signal x (t) to be analyzed is added hanning window, N+L point (getting L=N/4) of sampling, top n point forms sequence 1, x w1(n); Rear N number of formation sequence 2, x w2(n); Respectively FFT computing is done to two sequences;
1.2) peak value spectral line i under each harmonic frequency is determined in two sequences lwith two large position of spectral line i l+ 1, i l-1, and the amplitude X obtaining 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) y is made 1=X w1(i l-1), y 2=X w1(i l), y 3=X w1(i l+ 1), and by calculate β l;
1.4) by β lsubstitute into following formula and calculate harmonic frequency deviation δ l;
&delta; l = 0.0011 &beta; l 7 - 0.0079 &beta; l 6 + 0.0181 &beta; l 5 _ 0.0048 &beta; l 4 - 0.0785 &beta; l 3 + 0.0015 &beta; l 2 + 0.6665 &beta; l - 0.5
1.5) by δ lsubstitute into the amplitude A that following formula calculates each harmonic l;
A l = N ( y 1 + 2 y 2 + y 3 ) ( 1.4726 + 0.5891 &delta; l + 0.7221 &delta; l 2 + 0.288 &delta; l 3 + 0.2039 &delta; l 4 + 0.0895 &delta; l 5 + 0.049 &delta; l 6 )
Utilize composite iterative algorithm determination harmonic phase and frequency parameter, step is as follows:
2.1) by formula calculate the fundamental phase of acquisition two sequences and as the initial value of complex iteration; Specification error limit ε, 0 → k;
2.2) by calculate fundamental frequency deviation, and by calculate other each harmonic frequency departures;
2.3) if 0≤δ l< 0.5, then get μ=1; Otherwise get μ=-1, Amplitude Ration is calculated to sequence 1 amplitude Ration is calculated to sequence 2
2.4) Amplitude Ration obtained in conjunction with upper step again by following formula calculates the iterative value of two sequence each harmonic phase angles
2.5) if meet then stop iteration; Otherwise, k+1 → k, goes to step 2.2);
2.6) phase angle of the last iteration of each harmonic is exported and frequency departure and by formula f l=(i l+ δ l) Δ f, l=1,2 ..., K, obtains the Frequency Estimation result f of each harmonic l k; Wherein Δ f=1/NT s, N is sampling number, T sfor the sampling period.
3. High-precision harmonic parameter estimation method according to claim 1 and 2, is characterized in that: the described limits of error is set as ε=10 -4.
CN201310354109.0A 2013-08-09 2013-08-09 A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm CN103399203B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201310354109.0A CN103399203B (en) 2013-08-09 2013-08-09 A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201310354109.0A CN103399203B (en) 2013-08-09 2013-08-09 A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm

Publications (2)

Publication Number Publication Date
CN103399203A CN103399203A (en) 2013-11-20
CN103399203B true CN103399203B (en) 2015-08-26

Family

ID=49562860

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201310354109.0A CN103399203B (en) 2013-08-09 2013-08-09 A kind of High-precision harmonic parameter estimation method based on composite iterative algorithm

Country Status (1)

Country Link
CN (1) CN103399203B (en)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103630743B (en) * 2013-12-16 2015-12-30 电子科技大学 The method of a kind of heterodyne system spectrum analyzer frequency correction
CN104849548B (en) * 2015-06-01 2018-05-25 海南大学 A kind of electric system instantaneous frequency monitoring method and system
CN104881394B (en) * 2015-06-03 2017-08-18 河南理工大学 Single-mode system harmonic balance subtraction unit
CN105758840B (en) * 2016-03-01 2019-01-29 华中科技大学 A method of molecular orbit tomography is realized using higher hamonic wave amplitude
CN107305223B (en) * 2016-04-19 2019-12-10 天津大学 Improved phase difference frequency estimation method
CN106841774B (en) * 2017-01-24 2019-10-25 海南大学 A kind of power system frequency acquisition methods and system based on double-layer lap generation
CN106918741B (en) * 2017-03-02 2019-04-23 浙江大学 Adaptively sampled phase difference correction method applied to frequency wide swings power grid
CN106970264B (en) * 2017-03-02 2020-02-21 浙江大学 Improved phase difference correction method considering power grid frequency change rate
CN109495195A (en) * 2018-10-15 2019-03-19 中国人民解放军战略支援部队信息工程大学 The combined estimation method and device of radio communication PCMA signal amplitude

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0921835A (en) * 1995-07-05 1997-01-21 Mitsubishi Heavy Ind Ltd Voltage phase detection device
CN103018557A (en) * 2012-11-30 2013-04-03 合肥工业大学 Normalization master-slave type harmonic wave and inter-harmonic wave real-time analysis method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH0921835A (en) * 1995-07-05 1997-01-21 Mitsubishi Heavy Ind Ltd Voltage phase detection device
CN103018557A (en) * 2012-11-30 2013-04-03 合肥工业大学 Normalization master-slave type harmonic wave and inter-harmonic wave real-time analysis method

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
FFT高精度谐波检测方法研究;段小华等;《江苏电器》;20071231(第6期);第54-58页 *
基于三谱线插值FFT的电力谐波分析算法;牛胜锁等;《中国电机工程学报》;20120605;第32卷(第16期);第130-136页 *
基于加汉宁窗的FFT高精度谐波检测改进算法;王刘旺等;《电力系统保护与控制》;20121216;第40卷(第24期);第28-33页 *
基于迭代傅里叶变换的交流电压基波幅值与相位检测;张华军等;《电工电气》;20130131(第1期);第14-16页,34页 *

Also Published As

Publication number Publication date
CN103399203A (en) 2013-11-20

Similar Documents

Publication Publication Date Title
Asiminoaei et al. A digital controlled PV-inverter with grid impedance estimation for ENS detection
CN102435851B (en) Method for measuring zero-sequence parameters of double-circuit transmission lines
Suonan et al. Distance protection for HVDC transmission lines considering frequency-dependent parameters
CN101902195B (en) Method for automatically calibration of excitation system modelling and PSS optimization
CN101806832B (en) Measuring method for frequencies of low-frequency signals
CN101603985B (en) Method for measuring sine signal with high accuracy
CN104391178B (en) A kind of time shift phase difference stable state harmonic signal bearing calibration based on Nuttall windows
CN102510263B (en) Method for identifying practical parameters of synchronous generator on basis of load rejection test and numerical difference
CN103257271B (en) A kind of micro-capacitance sensor harmonic wave based on STM32F107VCT6 and m-Acetyl chlorophosphonazo pick-up unit and detection method
CN101701982B (en) Method for detecting harmonic waves of electric system based on window and interpolated FFT
CN101807795B (en) Method for forming electric energy metering simulation system and device thereof
US20100072978A1 (en) Synchrophasor measuring device and inter-bus-line phase angle difference measurement unit using the same
CN101635457B (en) Electric network parameter estimation method based on parameter sensitivity of state estimation residual error
CN202339381U (en) Harmonic electric energy metering system based on Nuttall self-convolution window weighed FFT (Fast Fourier Transform)
CN102288807B (en) Method for measuring electric network voltage flicker
CN103576002B (en) A kind of computing method of capacitive insulator arrangement dielectric loss angle
CN102163844B (en) Method for detecting state of power system based on phasor measurement unit (PMU)
CN106443246B (en) The on-line identification method of small interference stability parameter based on PMU metric data
Mai et al. An adaptive dynamic phasor estimator considering DC offset for PMU applications
Zeng et al. Harmonic phasor analysis based on improved FFT algorithm
CN102435815B (en) Operating method of resistive current on-line monitoring system of metal oxide arrester (MOA)
CN104897961B (en) Three spectral line interpolation FFT harmonic analysis methods and system based on cross multiplication window function
CN101701984B (en) Fundamental wave and harmonic wave detecting method based on three-coefficient Nuttall windowed interpolation FFT
CN101261292A (en) Base wave and harmonic detection method based on fiver item Rife-Vincent(1)window double spectral line interpolation FFT
CN104020370B (en) The power transformer interior fault diagnostic method monitored based on virtual parameter change

Legal Events

Date Code Title Description
PB01 Publication
C06 Publication
SE01 Entry into force of request for substantive examination
C10 Entry into substantive examination
GR01 Patent grant
C14 Grant of patent or utility model
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20150826

Termination date: 20180809

CF01 Termination of patent right due to non-payment of annual fee