CN102680785A - Synchronous phasor measurement method based on self-adoption variable window - Google Patents

Abstract

The invention provides a synchronous phasor measurement method based on a self-adoption variable window, which is used for the field of a synchronous phase measurement of a power system. When being variable, a signal frequency is measured in real time, and the signal frequency in the next sampling time is estimated. The method used by the invention is used for keeping a sampling frequency unchanged rather than adjusting the sampling frequency fs. A data length N belonging to current calculation can be automatically adjusted according to an estimation frequency, so as to keep N*Ts=T (Ts=1/fs is a fixed sampling period, fs is the sampling frequency, T=1/f is an actual signal period value, and f is a signal frequency). Thus, the error compensation precision is improved; the frequency estimation precision is greatly improved. The synchronous phasor measurement method provided by the invention is applicable to synchronous phasor measurement calculation and frequency calculation correction in a ship power system.

Description

Become the synchronous phasor measurement method of window based on self-adaptation

Technical field

The present invention relates to the synchronous phasor measurement field, specifically is a kind of synchronous phasor measurement method based on self-adaptation change window that is used for shipboard power system.

Background technology

Existing most of phasor measurement algorithm, its theoretical foundation all is DFT DFT, when the skew of electrical network medium frequency was little, no matter to the measurement of amplitude and phase place, DFT all can provide satisfied result.But when frequency shift (FS) was big, because the spectrum leakage that non-synchronous sampling causes, the precision of phasor measurement can descend rapidly.

Summary of the invention

The present invention provides a kind of synchronous phasor measurement method that becomes window based on self-adaptation; Some shortcomings to present phasor measurement algorithm; Improve on this basis; Improve the error compensation precision, improved the Frequency Estimation precision greatly, be applicable to that synchronous phasor measurement calculates and the frequency computation part correction in the shipboard power system.

A kind of synchronous phasor measurement method based on self-adaptation change window comprises the steps:

(1) variable-definition

φ: be instantaneous phasor phase angle

X, x: instantaneous phasor amplitude

N: each power frequency period sampling number

f _s: SF

(2) the actual frequency f of next sampling instant signal estimates

$f=\frac{\mathrm{\Δ\Φ}}{2\mathrm{\π\Δt}}$

$f=\frac{\mathrm{d\Φ}}{2\mathrm{\πdt}}$

$f=\frac{3\mathrm{\Φ}\left(t\right)+-4\mathrm{\Φ}(t-\mathrm{\Δt})+\mathrm{\Φ}(t-2\mathrm{\Δt})}{4\mathrm{\π\Δt}}$

Φ (t) is a t phase place constantly, and Φ (t-Δ t), Φ (t-2 Δ t) are its preceding two moment phase places;

(3) t _kAnd the phase place Φ in preceding two moment ^(k), Φ ^(k-1), Φ ^(k-2), t _kFrequency f constantly ^(k):

$f^{\left(k\right)}=\frac{{3\mathrm{\Φ}}^{\left(k\right)}-4{\mathrm{\Φ}}^{(k-1)}+{\mathrm{\Φ}}^{(k-2)}}{4\mathrm{\π}{T}_s}$

Current frequency change rate f ^(k)'

$f^{\left(k\right)\&prime;}=\frac{{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">d</figure-callout>}}^2\mathrm{\&Phi;}}{2\mathrm{\&pi;}{\mathrm{<figure-callout\; id="2"\; label="dt"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">dt</figure-callout>}}^2}=\frac{{2\mathrm{\&Phi;}}^{\left(k\right)}-5{\mathrm{\&Phi;}}^{(k-1)}+{4\mathrm{\&Phi;}}^{(k-2)}-{\mathrm{\&Phi;}}^{(k-3)}}{2\mathrm{\&pi;}{\mathrm{<figure-callout\; id="2"\; label="T\; s"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">T</figure-callout>}}_s^{}}$

Frequency predication value

$f_p^{(k+1)}=f^{\left(k\right)}+f^{\left(k\right)\&prime;}g{T}_s=\frac{{7\mathrm{\&Phi;}}^{\left(k\right)}-14{\mathrm{\&Phi;}}^{(k-1)}+{9\mathrm{\&Phi;}}^{(k-2)}-2{\mathrm{\&Phi;}}^{(k-3)}}{4\mathrm{\&pi;}{T}_s}$

By the frequency predication value

Adjustment t _K+1Data window length N constantly _K+1:

Wherein [g] _EvenThe even number computing is got in expression.

(4) calculating of signal Synchronization phasor

If t _kSelected sample sequence constantly is { x ^(k)(n) }, its ring shift right sequence does Corresponding phasor does The ring shift left sequence does

Corresponding phasor does

A) work as N _K+1＞N _kThe time

$X^{(k+1)}=W_{N_{k+1}}^m{X}_{\mathrm{LS}}^{(k+1)};$

$X_{\mathrm{LS}}^{(k+1)}=\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

$=\frac{\sqrt{2}}{N_{k+1}}(\underset{n=0}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+\underset{n=N_k}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n)$

$\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}\underset{n=N_k}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

B) work as N _K+1=N _kThe time, use circulation DFT method

C) work as N _K+1＜N _kThe time

$X^{(k+1)}=W_{N_{k+1}}^{-m}{X}_{\mathrm{RS}}^{(k+1)}$

$X_{\mathrm{RS}}^{(k+1)}=\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

$=\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+(\underset{n=m}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

$+\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n)-\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

$\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

${\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}=\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

${\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}=-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

$\stackrel{r=n-N_{k+1}}{=}-\underset{r=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}(r+N_{k+1})W_{N_{k+1}}^{r+N_{k+1}}$

$=-\underset{r=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(r\right)W_{N_{k+1}}^r$

$=-\underset{n=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

Because N _K+1With N _kBe even number, N _K+1-N _k>=2, work as N _K+1-N _k=2 o'clock,

$X_{\mathrm{RS}}^{(k+1)}\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-2x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[x_{\mathrm{RS}}^{(k+1)}\left(0\right)-2x^{\left(k\right)}\left(0\right)$

$+(x_{\mathrm{RS}}^{(k+1)}\left(1\right)-2x^{\left(k\right)}\left(1\right))W_{N_{k+1}}^1]$

Work as N _K+1-N _k＞2 o'clock,

$X_{\mathrm{RS}}^{(k+1)}\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-2x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

$+\underset{n=m}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

Work as N like this _K+1≈ N _kThe time, use above-mentioned formula cause X ^(k)Recursion X ^(k+1)

(5) the total distortion degree THD of definition signal { x (n) } is:

${\mathrm{THD}}_x=\sqrt{\frac{P_x-{A_X}^2}{{A_{\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">X</figure-callout>}}}^2}}\&times;100\%=\sqrt{\frac{\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x\left(n\right)\right|}^2-{A_X}^2}{{{\mathrm{<figure-callout\; id="2"\; label="A\; X"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">A</figure-callout>}}_X}^2}}\&times;100\%$

A wherein _XBe the corresponding fundamental phasors amplitude of signal x (n);

Set up error-factor of influence numerical value partial differential (Error-Factors Partial Differential, EFPD) table.Each element definition in the EFPD table is the partial derivative of each measuring error to each factor of influence, and promptly each element is in the EFPD table:

$(\&PartialD;{\mathrm{\&delta;}}_A/\&PartialD;\left(\mathrm{\&Delta;f}\right),\&PartialD;{\mathrm{\&delta;}}_A/\&PartialD;\left(\mathrm{THD}\right),{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;\mathrm{\&Phi;},{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\left(\mathrm{\&Delta;f}\right),{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\left(\mathrm{THD}\right),{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\mathrm{\&Phi;})$

δ wherein _ABe amplitude error, δ _ΦBe phase error;

(6) t _kThe initial instantaneous phasor of constantly using RDFT to calculate signal is X ^(k), just obtained the estimated value A of amplitude ^(k)Estimated value Φ with phase place ^(k), estimate frequency f ^(k), estimate again

${\mathrm{THD}}_x^{\left(k\right)}=\sqrt{\frac{\frac{1}{N}\stackrel{N-1}{\underset{n=0}{\mathrm{\&Sigma;}}}{\left|x^{\left(k\right)}\left(n\right)\right|}^2-{\left(A_X^{\left(k\right)}\right)}^2}{{\left(A_X^{\left(k\right)}\right)}^2}}\&times;100\%$

$=\sqrt{\frac{\frac{1}{N}(\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x^{(k-1)}\left(n\right)\right|}^2+{\left|x^{\left(k\right)}(N-1)\right|}^2-{\left|x^{(k-1)}\left(0\right)\right|}^2)-{\left(A_X^{\left(k\right)}\right)}^2}{{\left(A_X^{\left(k\right)}\right)}^2}}\&times;100\%$

Compensated instantaneous phasor magnitude

and the phase

$A_{\mathrm{cmp}}^{\left(k\right)}=A^{\left(k\right)}+{\mathrm{\&delta;}}_A^{\left(k\right)}+{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;{\left(\mathrm{\&Delta;f}\right)}^{\left(k\right)}g({\mathrm{\&Delta;f}}^{\left(k\right)}-{\mathrm{\&Delta;f}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;{\left(\mathrm{THD}\right)}^{k}g({\mathrm{THD}}^{\left(k\right)}-{\mathrm{THD}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_A/{\&PartialD;\mathrm{\&Phi;}}^{\left(k\right)}g({\mathrm{\&Phi;}}^{\left(k\right)}-{\mathrm{\&Phi;}}_T^{\left(k\right)})$

${\mathrm{\&Phi;}}_{\mathrm{cmp}}^{\left(k\right)}={\mathrm{\&Phi;}}^{\left(k\right)}+{\mathrm{\&delta;}}_{\mathrm{\&Phi;}}^{\left(k\right)}+{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;{\left(\mathrm{\&Delta;f}\right)}^{\left(k\right)}g({\mathrm{\&Delta;f}}^{\left(k\right)}-{\mathrm{\&Delta;f}}_T^{\left(k\right)})$

$+\&PartialD;{\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;{\left(\mathrm{THD}\right)}^{k}g({\mathrm{THD}}^{\left(k\right)}-{\mathrm{THD}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/{\&PartialD;\mathrm{\&Phi;}}^{\left(k\right)}g({\mathrm{\&Phi;}}^{\left(k\right)}-{\mathrm{\&Phi;}}_T^{\left(k\right)})$

(7) by t _kConstantly

With E [ΔΦ _Survey(k)] replacement ΔΦ, the Frequency Estimation formula:

The present invention is used for the synchronous phase measuring in power system field, and when signal frequency changed, the frequency of the next sampling instant of the frequency of measuring-signal, and estimated signal in real time, the method that the present invention adopts were not to regulate SF f _s, but keep SF constant, regulate automatically according to estimated frequency and include the data length N in the current calculating in, to keep NT _s=T (T _s=1/f _sBe fixing sampling period, f _sBe SF, T=1/f is the actual cycle value of signal, and f is a signal frequency), improved the error compensation precision, improved the Frequency Estimation precision greatly, be applicable to that synchronous phasor measurement calculates and the frequency computation part correction in the shipboard power system.

Description of drawings

Fig. 1 is the estimated frequency error of the present invention figure that is scattered;

The maximum error of measuring of the improved self-adapting data window of Fig. 2 the present invention length method amplitude figure that is scattered;

The maximum error of measuring of the improved self-adapting data window of Fig. 3 the present invention length method phase place figure that is scattered.

Embodiment

To combine the accompanying drawing among the present invention below, the technical scheme among the present invention will be carried out clear, intactly description.

The adaptively sampled synchronized phasor method that the present invention adopted is used for the shipboard power system synchronous phasor measurement, and this measurement result is carried out the difference compensation, has greatly improved the precision of synchronous phasor measurement, and said step phase metering method comprises the steps:

(1) variable-definition

φ: be instantaneous phasor phase angle

X, x: instantaneous phasor amplitude

N: each power frequency period sampling number

f _s: SF

(2) the actual frequency f of next sampling instant signal estimates

Make T _s=Δ t, then:

$f=\frac{\mathrm{\&Delta;\&Phi;}}{2\mathrm{\&pi;\&Delta;t}}$

When Δ t → 0, can obtain

$f=\frac{\mathrm{d\&Phi;}}{2\mathrm{\&pi;dt}}$

Utilize the consequent difference formula of newton's second order, can get

$f=\frac{3\mathrm{\&Phi;}\left(t\right)+-4\mathrm{\&Phi;}(t-\mathrm{\&Delta;t})+\mathrm{\&Phi;}(t-2\mathrm{\&Delta;t})}{4\mathrm{\&pi;\&Delta;t}}$

(3) self-adaptation adjustment data window length makes up

Obtain t according to the DFT algorithm _kPhasor X constantly ^(k), and obtain t _kAnd the phase place Φ in preceding two moment ^(k), Φ ^(k-1), Φ ^(k-2).Obtain t by formula (2.3.4) then _kFrequency f constantly ^(k):

$f^{\left(k\right)}=\frac{{3\mathrm{\&Phi;}}^{\left(k\right)}-4{\mathrm{\&Phi;}}^{(k-1)}+{\mathrm{\&Phi;}}^{(k-2)}}{4\mathrm{\&pi;}{T}_s}$

Differential formulas by asking second derivative in the consequent differential formulas of Lagrange obtains current frequency change rate f ^(k)'

$f^{\left(k\right)\&prime;}=\frac{{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">d</figure-callout>}}^2\mathrm{\&Phi;}}{2\mathrm{\&pi;}{\mathrm{<figure-callout\; id="2"\; label="dt"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">dt</figure-callout>}}^2}=\frac{{2\mathrm{\&Phi;}}^{\left(k\right)}-5{\mathrm{\&Phi;}}^{(k-1)}+{4\mathrm{\&Phi;}}^{(k-2)}-{\mathrm{\&Phi;}}^{(k-3)}}{2\mathrm{\&pi;}{T}_s^2}$

So

$f_p^{(k+1)}=f^{\left(k\right)}+f^{\left(k\right)\&prime;}g{T}_s=\frac{{7\mathrm{\&Phi;}}^{\left(k\right)}-14{\mathrm{\&Phi;}}^{(k-1)}+{9\mathrm{\&Phi;}}^{(k-2)}-2{\mathrm{\&Phi;}}^{(k-3)}}{4\mathrm{\&pi;}{T}_s}$

By the frequency predication value Adjustment t _K+1Data window length N constantly _K+1:

Wherein [g] _EvenThe even number computing is got in expression.

(4) calculating of signal Synchronization phasor

If t _kSelected sample sequence constantly is { x ^(k)(n) }, its ring shift right sequence does

Corresponding phasor does

The ring shift left sequence does Corresponding phasor does

A) work as N _K+1＞N _kThe time, { x ^(k+1)(n) } m=(N that need move to left _K+1-N _kA)/2-1 sampling period obtains

Character by DFT has:

$X^{(k+1)}=W_{N_{k+1}}^m{X}_{\mathrm{LS}}^{(k+1)}$

Because at this moment

Preceding N _kIndividual point and { x ^(k)(n) } corresponding point is identical, so just have

$X_{\mathrm{LS}}^{(k+1)}=\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

$=\frac{\sqrt{2}}{N_{k+1}}(\underset{n=0}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+\underset{n=N_k}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n)$

$\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}\underset{n=N_k}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

B) work as N _K+1=N _kThe time, use circulation DFT method

C) work as N _K+1＜N _kThe time, { x ^(k+1)(n) } m=(N that need move to right _k-N _K+1A)/2+1 sampling period obtains

Character by DFT has:

$X^{(k+1)}=W_{N_{k+1}}^{-m}{X}_{\mathrm{RS}}^{(k+1)}$

Because at this moment

M+1 o'clock to N _K+1Individual point and { x ^(k)(n) } corresponding point is identical, so just have

$X_{\mathrm{RS}}^{(k+1)}=\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

$=\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+(\underset{n=m}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

$+\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n)-\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

$\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

Wherein,

${\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}=\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

${\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}=-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

$\stackrel{r=n-N_{k+1}}{=}-\underset{r=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}(r+N_{k+1})W_{N_{k+1}}^{r+N_{k+1}}$

$=-\underset{r=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(r\right)W_{N_{k+1}}^r$

$=-\underset{n=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

Because N _K+1With N _kBe even number.So N _K+1-N _k>=2.So, work as N _K+1-N _k=2 o'clock,

$X_{\mathrm{RS}}^{(k+1)}\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-2x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[x_{\mathrm{RS}}^{(k+1)}\left(0\right)-2x^{\left(k\right)}\left(0\right)$

$+(x_{\mathrm{RS}}^{(k+1)}\left(1\right)-2x^{\left(k\right)}\left(1\right))W_{N_{k+1}}^1]$

Work as N _K+1-N _k＞2 o'clock,

$X_{\mathrm{RS}}^{(k+1)}\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-2x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

$+\underset{n=m}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

Work as N like this _K+1≈ N _kThe time, use above-mentioned formula cause X ^(k)Recursion X ^(k+1)

(5) set up off-line error compensation tables and error-factor of influence numerical value partial differential table.

The improved Error Compensation Algorithm based on interpolation of this joint is the multiple influence of having taken all factors into consideration frequency shift (FS), distorted signals and amplitude fluctuations, through setting up error compensation tables, initial measurement result is tabled look-up earlier and interpolation realization error compensation.Before PMU formally measures; Produce one group of sufficiently long normalized sinusoidal sequence of amplitude that contains other composition through certain measure (waveform generator is through AD conversion or software simulation); Its frequency shift (FS) Δ f (in ± 2.5Hz) and THD (in 0～5% scope) are known in advance; PMU adopts certain algorithm (like DFT, average, self-adapting data window length etc.) that these group data are handled then; Provide the statistics of corresponding phasor measurements (amplitude, phase place) after intact all the input data of algorithm process, the phasor of this result and known array is compared, obtain to amplitude error under the phase sequence of fixed step size and phase error.Give the phase sequence of fixed step size be meant span (π, π] an equal difference phase sequence in interval.As phase place Φ get respectively 5 π/6 ,-4 π/6 ..., 4 π/6,5 π/6, π }.

The total distortion degree THD of definition signal { x (n) } is:

${\mathrm{THD}}_x=\sqrt{\frac{P_x-{A_X}^2}{{A_{\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">X</figure-callout>}}}^2}}\&times;100\%=\sqrt{\frac{\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x\left(n\right)\right|}^2-{A_X}^2}{{A_{\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="DEST\_PATH\_HDA00001748125200011.png,DEST\_PATH\_HDA00001748125200012.png"\; state="\{\{state\}\}">X</figure-callout>}}}^2}}\&times;100\%$

A wherein _XBe the corresponding fundamental phasors amplitude of signal x (n).

If (Δ f THD) is one group of parameter of given sinusoidal sequence, existing hypothesis give the Δ f sequence of fixed step size be 2.5 ,-2.4 ..., 2.3,2.4,2.5}, be for the THD sequence of fixed step size 0,0.5% ..., 4.0%, 4.5%, 5.0%}.The parameter matrix table that then defines input signal does

$\left[\begin{array}{ccccc}(-\mathrm{2.5,0}\%)& (-\mathrm{2.5,0.5}\%)& L& L& (-\mathrm{2.5,5.0}\%)\\ (-\mathrm{2.4,0}\%)& (-\mathrm{2.4,0.5}\%)& L& L& (-\mathrm{2.4,5.0}\%)\\ M& O& O& O& M\\ M& O& O& O& M\\ (\mathrm{2.5,0}\%)& (\mathrm{2.5,0.5}\%)& L& L& (\mathrm{2.5,5.0}\%)\end{array}\right]$

For each element in the parameter matrix table (being each group parameter), producing corresponding list entries respectively, and obtain to amplitude error and phase error under the Φ sequence of fixed step size.

So just, obtain the error compensation tables of a three-dimensional, the factor of influence variable of each dimension representative is that (Φ), each element in the table is a two-dimensional array for Δ f, THD, respectively corresponding amplitude error δ _AWith phase error δ _Φ

After error compensation tables is set up, set up error-factor of influence numerical value partial differential (Error-Factors Partial Differential, EFPD) table according to it.Each element definition in the EFPD table is the partial derivative of each measuring error to each factor of influence, and promptly each element is in the EFPD table:

$(\&PartialD;{\mathrm{\&delta;}}_A/\&PartialD;\left(\mathrm{\&Delta;f}\right),\&PartialD;{\mathrm{\&delta;}}_A/\&PartialD;\left(\mathrm{THD}\right),{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;\mathrm{\&Phi;},{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\left(\mathrm{\&Delta;f}\right),{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\left(\mathrm{THD}\right),{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\mathrm{\&Phi;})$

For the each point that is designated as down lower boundary in the error compensation tables, be ± each point of 5Hz like Δ f.Obtain each corresponding element in the EFPD table with 3 forward difference formula of newton.For example ask the formula of

to be:

$\frac{{\&PartialD;\mathrm{\&delta;}}_A}{\&PartialD;\left(\mathrm{\&Delta;f}\right)}\&ap;\frac{-3{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)+4{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f}+h_{\mathrm{\&Delta;f}},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)-{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f}+2h_{\mathrm{\&Delta;f}},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)}{2h_{\mathrm{\&Delta;f}}}$

THD _CAnd Φ _CThese two subscripts keep constant in the expression error compensation tables, h _{Δ f}Δ f subscript increases by 1 in the expression error compensation tables.EFPD formula between other variable is consistent.

For the each point that is designated as non-border in the error compensation tables down; Available newton's central-difference formula is asked each corresponding element in the EFPD table, for example asks the formula of

to be:

$\frac{{\&PartialD;\mathrm{\&delta;}}_A}{\&PartialD;\left(\mathrm{\&Delta;f}\right)}\&ap;\frac{{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f}+h_{\mathrm{\&Delta;f}},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)-{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f}-h_{\mathrm{\&Delta;f}},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)}{{2h}_{\mathrm{\&Delta;f}}}$

For the each point that is designated as the coboundary in the error compensation tables down, the consequent difference formula of available newton is asked each corresponding element in the EFPD table:

$\frac{{\&PartialD;\mathrm{\&delta;}}_A}{\&PartialD;\left(\mathrm{\&Delta;f}\right)}\&ap;\frac{3{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)-4{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f}-h_{\mathrm{\&Delta;f}},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)+{\mathrm{\&delta;}}_A(\mathrm{\&Delta;f}-2h_{\mathrm{\&Delta;f}},{\mathrm{THD}}_C,{\mathrm{\&Phi;}}_C)}{2h_{\mathrm{\&Delta;f}}}$

(6) online phasor calculation and interpolation error compensation

t _kThe initial instantaneous phasor of constantly using RDFT to calculate signal is X ^(k), just obtained the estimated value A of amplitude ^(k)Estimated value Φ with phase place ^(k)Estimate frequency f ^(k), estimate with following formula again

${\mathrm{THD}}_x^{\left(k\right)}=\sqrt{\frac{\frac{1}{N}\stackrel{N-1}{\underset{n=0}{\mathrm{\&Sigma;}}}{\left|x^{\left(k\right)}\left(n\right)\right|}^2-{\left(A_X^{\left(k\right)}\right)}^2}{{\left(A_X^{\left(k\right)}\right)}^2}}\&times;100\%$

$=\sqrt{\frac{\frac{1}{N}(\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x^{(k-1)}\left(n\right)\right|}^2+{\left|x^{\left(k\right)}(N-1)\right|}^2-{\left|x^{(k-1)}\left(0\right)\right|}^2)-{\left(A_X^{\left(k\right)}\right)}^2}{{\left(A_X^{\left(k\right)}\right)}^2}}\&times;100\%$

So obtain factor of influence variable (Δ f ^(k), THD ^(k), Φ ^(k)) estimated value, estimated value is looked into error compensation tables and EFPD table thus again, obtains nearest subscript array and the element value of two table middle distances, is made as

And

Instantaneous phasor amplitude after being compensated according to following linear interpolation formula then

And phase place

$A_{\mathrm{cmp}}^{\left(k\right)}=A^{\left(k\right)}+{\mathrm{\&delta;}}_A^{\left(k\right)}+{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;{\left(\mathrm{\&Delta;f}\right)}^{\left(k\right)}g({\mathrm{\&Delta;f}}^{\left(k\right)}-{\mathrm{\&Delta;f}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;{\left(\mathrm{THD}\right)}^{k}g({\mathrm{THD}}^{\left(k\right)}-{\mathrm{THD}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_A/{\&PartialD;\mathrm{\&Phi;}}^{\left(k\right)}g({\mathrm{\&Phi;}}^{\left(k\right)}-{\mathrm{\&Phi;}}_T^{\left(k\right)})$

${\mathrm{\&Phi;}}_{\mathrm{cmp}}^{\left(k\right)}={\mathrm{\&Phi;}}^{\left(k\right)}+{\mathrm{\&delta;}}_{\mathrm{\&Phi;}}^{\left(k\right)}+{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;{\left(\mathrm{\&Delta;f}\right)}^{\left(k\right)}g({\mathrm{\&Delta;f}}^{\left(k\right)}-{\mathrm{\&Delta;f}}_T^{\left(k\right)})$

$+\&PartialD;{\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;{\left(\mathrm{THD}\right)}^{k}g({\mathrm{THD}}^{\left(k\right)}-{\mathrm{THD}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/{\&PartialD;\mathrm{\&Phi;}}^{\left(k\right)}g({\mathrm{\&Phi;}}^{\left(k\right)}-{\mathrm{\&Phi;}}_T^{\left(k\right)})$

(7) based on the frequency estimation algorithm of adding up

By t _kMoment E [ΔΦ _Survey(k)] formula has

This formula only need be done subtraction and multiplication operation just can be accomplished E [ΔΦ _Survey(k)] calculating.

Consider with E [ΔΦ _Survey(k)] replacement ΔΦ.Obtain new Frequency Estimation formula:

ΔΦ _SurveyThe frequency of vibration is two times of survey frequency.Therefore to all ΔΦs in the measuring period _SurveyAsk average with to all ΔΦs in half measuring period _SurveyAsk average effect basic identical.And ask ΔΦ in half measuring period _SurveyAverage, the real-time of algorithm is better.Therefore formula can change into

Through simulation result, estimated frequency error is as shown in Figure 1 with the curve that Δ f changes.It is thus clear that when Δ f was in ± 3Hz, the absolute value of estimated frequency error can be greater than 0.3Hz.When Δ f was in ± 5Hz, the absolute value of estimated frequency error can be greater than 0.5Hz.

After using new Frequency Estimation formula, the maximum error of measuring of self-adapting data window length method amplitude is as shown in Figure 2, and the maximum error of measuring of phase place is as shown in Figure 3.Visible by figure, when Δ f was in ± 5Hz, the amplitude relative error can be greater than 0.6%, and the phase place absolute error is in 2.5 °.Special, when Δ f was in ± 3Hz, the absolute error of phase place was controlled in 1 °, satisfies the needs of the requirement and the practical application of variety of protocol so fully.

Claims (1)

1. the synchronous phasor measurement method based on self-adaptation change window is characterized in that: comprise the steps:

(1) variable-definition

φ: be instantaneous phasor phase angle;

X, x: instantaneous phasor amplitude;

N: each power frequency period sampling number;

f _s: SF;

(2) the actual frequency f of next sampling instant signal estimates

$f=\frac{\mathrm{\&Delta;\&Phi;}}{2\mathrm{\&pi;\&Delta;t}}$

$f=\frac{\mathrm{d\&Phi;}}{2\mathrm{\&pi;dt}}$

$f=\frac{3\mathrm{\&Phi;}\left(t\right)+-4\mathrm{\&Phi;}(t-\mathrm{\&Delta;t})+\mathrm{\&Phi;}(t-2\mathrm{\&Delta;t})}{4\mathrm{\&pi;\&Delta;t}}$

Φ (t) is a t phase place constantly, and Φ (t-Δ t), Φ (t-2 Δ t) are its preceding two moment phase places;

(3) t _kAnd the phase place Φ in preceding two moment ^(k), Φ ^(k-1), Φ ^(k-2), t _kFrequency f constantly ^(k):

$f^{\left(k\right)}=\frac{{3\mathrm{\&Phi;}}^{\left(k\right)}-4{\mathrm{\&Phi;}}^{(k-1)}+{\mathrm{\&Phi;}}^{(k-2)}}{4\mathrm{\&pi;}{T}_s}$

Current frequency change rate f ^(k)'

$f^{\left(k\right)\&prime;}=\frac{d^2\mathrm{\&Phi;}}{2\mathrm{\&pi;}{\mathrm{dt}}^2}=\frac{{2\mathrm{\&Phi;}}^{\left(k\right)}-5{\mathrm{\&Phi;}}^{(k-1)}+{4\mathrm{\&Phi;}}^{(k-2)}-{\mathrm{\&Phi;}}^{(k-3)}}{2\mathrm{\&pi;}{T}_s^2}$

Frequency predication value

$f_p^{(k+1)}=f^{\left(k\right)}+f^{\left(k\right)\&prime;}g{T}_s=\frac{{7\mathrm{\&Phi;}}^{\left(k\right)}-14{\mathrm{\&Phi;}}^{(k-1)}+{9\mathrm{\&Phi;}}^{(k-2)}-2{\mathrm{\&Phi;}}^{(k-3)}}{4\mathrm{\&pi;}{T}_s}$

By the frequency predication value

Adjustment t _K+1Data window length N constantly _K+1:

Wherein [g] _EvenThe even number computing is got in expression.

(4) calculating of signal Synchronization phasor:

If t _kSelected sample sequence constantly is { x ^(k)(n) }, its ring shift right sequence does

Corresponding phasor does

The ring shift left sequence does

Corresponding phasor does

Concrete,

A) work as N _K+1＞N _kThe time

$X^{(k+1)}=W_{N_{k+1}}^m{X}_{\mathrm{LS}}^{(k+1)};$

$X_{\mathrm{LS}}^{(k+1)}=\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

$=\frac{\sqrt{2}}{N_{k+1}}(\underset{n=0}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+\underset{n=N_k}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n)$

$\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}\underset{n=N_k}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{LS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

B) work as N _K+1=N _kThe time, use circulation DFT method;

C) work as N _K+1＜N _kThe time

$X^{(k+1)}=W_{N_{k+1}}^{-m}{X}_{\mathrm{RS}}^{(k+1)};$

$X_{\mathrm{RS}}^{(k+1)}=\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n$

$=\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+(\underset{n=m}{\overset{N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{RS}}^{(k+1)}\left(n\right)W_{N_{k+1}}^n+\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

$+\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n)-\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

$\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}];$

${\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}=\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n;$

${\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}=-\underset{n=N_{k+1}}{\overset{N_k-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n$

$\stackrel{r=n-N_{k+1}}{=}-\underset{r=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}(r+N_{k+1})W_{N_{k+1}}^{r+N_{k+1}}$

$=-\underset{r=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(r\right)W_{N_{k+1}}^r$

$=-\underset{n=0}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n;$

Because N _K+1With N _kBe even number, N _K+1-N _k>=2, work as N _K+1-N _k=2 o'clock,

$X_{\mathrm{RS}}^{(k+1)}\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-2x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[x_{\mathrm{RS}}^{(k+1)}\left(0\right)-2x^{\left(k\right)}\left(0\right)$

$+(x_{\mathrm{RS}}^{(k+1)}\left(1\right)-2x^{\left(k\right)}\left(1\right))W_{N_{k+1}}^1]$

Work as N _K+1-N _k＞2 o'clock,

$X_{\mathrm{RS}}^{(k+1)}\&ap;X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[{\mathrm{\&Delta;}}_1{X}_{\mathrm{RS}}^{(k+1)}+{\mathrm{\&Delta;}}_2{X}_{\mathrm{RS}}^{(k+1)}]$

$=X^{\left(k\right)}+\frac{\sqrt{2}}{N_{k+1}}[\underset{n=0}{\overset{m-1}{\mathrm{\&Sigma;}}}(x_{\mathrm{RS}}^{(k+1)}\left(n\right)-2x^{\left(k\right)}\left(n\right))W_{N_{k+1}}^n$

$+\underset{n=m}{\overset{N_k-N_{k+1}-1}{\mathrm{\&Sigma;}}}{x}^{\left(k\right)}\left(n\right)W_{N_{k+1}}^n]$

Work as N like this _K+1≈ N _kThe time, use above-mentioned formula cause X ^(k)Recursion X ^(k+1)

(5) the total distortion degree THD of definition signal { x (n) } is:

${\mathrm{THD}}_x=\sqrt{\frac{P_x-{A_X}^2}{{A_X}^2}}\&times;100\%=\sqrt{\frac{\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x\left(n\right)\right|}^2-{A_X}^2}{{A_X}^2}}\&times;100\%$

A wherein _XBe the corresponding fundamental phasors amplitude of signal x (n);

Set up error-factor of influence numerical value partial differential (Error-Factors Partial Differential, EFPD) table.Each element definition in the EFPD table is the partial derivative of each measuring error to each factor of influence, and promptly each element is in the EFPD table:

$(\&PartialD;{\mathrm{\&delta;}}_A/\&PartialD;\left(\mathrm{\&Delta;f}\right),\&PartialD;{\mathrm{\&delta;}}_A/\&PartialD;\left(\mathrm{THD}\right),{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;\mathrm{\&Phi;},{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\left(\mathrm{\&Delta;f}\right),{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\left(\mathrm{THD}\right),{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;\mathrm{\&Phi;})$

δ wherein _ABe amplitude error, δ _ΦBe phase error;

(6) t _kThe initial instantaneous phasor of constantly using RDFT to calculate signal is X ^(k), just obtained the estimated value A of amplitude ^(k)Estimated value Φ with phase place ^(k), estimate frequency f ^(k), estimate again

${\mathrm{THD}}_x^{\left(k\right)}=\sqrt{\frac{\frac{1}{N}\stackrel{N-1}{\underset{n=0}{\mathrm{\&Sigma;}}}{\left|x^{\left(k\right)}\left(n\right)\right|}^2-{\left(A_X^{\left(k\right)}\right)}^2}{{\left(A_X^{\left(k\right)}\right)}^2}}\&times;100\%$

$=\sqrt{\frac{\frac{1}{N}(\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left|x^{(k-1)}\left(n\right)\right|}^2+{\left|x^{\left(k\right)}(N-1)\right|}^2-{\left|x^{(k-1)}\left(0\right)\right|}^2)-{\left(A_X^{\left(k\right)}\right)}^2}{{\left(A_X^{\left(k\right)}\right)}^2}}\&times;100\%$

Compensated instantaneous phasor magnitude

and the phase

$A_{\mathrm{cmp}}^{\left(k\right)}=A^{\left(k\right)}+{\mathrm{\&delta;}}_A^{\left(k\right)}+{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;{\left(\mathrm{\&Delta;f}\right)}^{\left(k\right)}g({\mathrm{\&Delta;f}}^{\left(k\right)}-{\mathrm{\&Delta;f}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_A/\&PartialD;{\left(\mathrm{THD}\right)}^{k}g({\mathrm{THD}}^{\left(k\right)}-{\mathrm{THD}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_A/{\&PartialD;\mathrm{\&Phi;}}^{\left(k\right)}g({\mathrm{\&Phi;}}^{\left(k\right)}-{\mathrm{\&Phi;}}_T^{\left(k\right)})$

${\mathrm{\&Phi;}}_{\mathrm{cmp}}^{\left(k\right)}={\mathrm{\&Phi;}}^{\left(k\right)}+{\mathrm{\&delta;}}_{\mathrm{\&Phi;}}^{\left(k\right)}+{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;{\left(\mathrm{\&Delta;f}\right)}^{\left(k\right)}g({\mathrm{\&Delta;f}}^{\left(k\right)}-{\mathrm{\&Delta;f}}_T^{\left(k\right)})$

$+\&PartialD;{\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/\&PartialD;{\left(\mathrm{THD}\right)}^{k}g({\mathrm{THD}}^{\left(k\right)}-{\mathrm{THD}}_T^{\left(k\right)})$

$+{\&PartialD;\mathrm{\&delta;}}_{\mathrm{\&Phi;}}/{\&PartialD;\mathrm{\&Phi;}}^{\left(k\right)}g({\mathrm{\&Phi;}}^{\left(k\right)}-{\mathrm{\&Phi;}}_T^{\left(k\right)})$

(7) by t _kConstantly

With E [ΔΦ _Survey(k)] replacement ΔΦ, the Frequency Estimation formula:

Images