CN1207679C - High-accuracy synchronous phasor measuring method - Google Patents

Abstract

The present invention relates to a synchronous phasor measuring method with high accuracy, which belongs to the technical field of automatic measurement in an electric power system. The present invention is characterized in that the influence of frequency deviation is taken into account so as to enable the measurement of a frequency and a phasor to have the high accuracy on the basis of a fixed interval sampling cycle DFT when a system frequency deviates at a rated frequency. A field test indicates that when a sampling frequency is 4800Hz and a phasor correction calculating frequency is 200Hz, the method has the steady-state precision indexes: phasor amplitude error is less than 0.5%, the phase angle error of the phasor is less than one DEG and frequency calculation error is less than 0.01Hz; the method has the dynamic precision indexes: phasor amplitude error is less than 1%, the phase angle error of the phasor is less than one DEG, and frequency calculation error is less than 0.01Hz. The synchronous phasor measuring method with high accuracy completely satisfies the requirements of technical specification for a real-time dynamic monitoring (controlling) system in an electric power system.

Description

A kind of synchronous phasor measurement method of compensation cycle discrete Fourier transform (DFT) error

Technical field

A kind of synchronous phasor measurement method of compensation cycle discrete Fourier transform (DFT) error belongs to electric system automatic measurement technology field.

Background technology

Traditional synchronized phasor algorithm is when signal frequency and rated frequency have deviation, and there are bigger fluctuation (error) in the phasor amplitude that obtains, phase angle and frequency.For improving phasor measurement unit (PMU, Phasor Measurement Unit) precision, this patent propose a kind of new synchronized phasor algorithm on the basis of fixed sampling interval technique cyclic dispersion Fourier transform, to overcome the shortcoming of traditional algorithm, obtain high-precision frequency and synchronized phasor.In recent years, along with the theoretical research of phasor measurement unit and carrying out in a deep going way of application, the raising of phasor measurement precision has more importance and urgency.

Existing measuring method is based on the cyclic dispersion Fourier transform, and in principle, this method is when the system frequency excursion rated frequency, and the error of calculation increases with side-play amount, thereby the measurement range that is suitable for is little, and precision is not ideal enough.When mains frequency skew 50Hz was far away, computational accuracy did not reach the requirement of " electric system real-time dynamic monitoring (control) systems technology standard ".

Summary of the invention

The object of the present invention is to provide a kind of synchronous phasor measurement method of compensation cycle discrete Fourier transform (DFT) error, this method has been taken into account the influence of frequency departure to the cyclic dispersion fourier transform algorithm, thereby has higher measuring accuracy when the system frequency excursion rated frequency; This method also is applicable to phasor measurement unit sample frequency and the different situation of phasor correction calculation frequency.

The present invention is a kind of synchronous phasor measurement method of compensation cycle discrete Fourier transform (DFT) error, it is characterized in that the method includes the steps of:

At first carry out initialization of variable:

N: the ripple AC sampling is counted weekly, the integer greater than 12;

Δ T: in the AC sampling time interval, initial value is the 20/N millisecond;

The initial zero clearing of each sampling buffer;

With Δ T serves as to gather three-phase voltage (or electric current) value at interval, when sampling number is less than N, then continues to wait for new sampled point; When sampling number is more than or equal to N, whenever obtain a new sampled point, it and preceding N-1 sampled point form a new data window;

To the N point data of each new data window, order is carried out following steps;

1) adopt the cyclic dispersion fourier transform algorithm, calculate the phasor of A, B, each phase of C three-phase, as centre frequency, result of calculation is not considered the influence of frequency departure to the cyclic dispersion Fourier transform with rated frequency, and computing formula is:

$\hat{x}\left(r\right)=\hat{x}(r-1)+j\frac{\sqrt{2}}{N}[\tilde{x}(r+N-1)-\tilde{x}(r-1)]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="C0314595200131.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}(r-1)};$

Just entered the last samples value of data window, $\tilde{x}(r+N-1)=\tilde{x}\left(k\right){|}_{k=r+N-1},$

Just withdrawed from the sampled value of data window, $\tilde{x}(r-1)=\tilde{x}\left(k\right){|}_{k=r-1};$

Wherein Be defined as follows:

X: signal effective value;

f ₀: rated frequency;

Δ f: the deviation of signal bias rated frequency;

: the initial phase angle of signal;

The synchronized phasor that last cyclic dispersion Fourier transform calculates is used for the cyclic dispersion Fourier transformation computation;

The synchronized phasor that this cyclic dispersion Fourier transform calculates;

2) utilize A, B, the synthetic positive and negative and zero sequence phasor of C three-phase phasor, computing formula is:

${\tilde{x}}_1\left(r\right)=\frac{1}{3}[{\tilde{x}}_A\left(r\right)+a{\tilde{x}}_B\left(r\right)+a^2{\tilde{x}}_C\left(r\right)];$

${\tilde{x}}_2\left(r\right)=\frac{1}{3}[{\tilde{x}}_A\left(r\right)+a^2{\tilde{x}}_B\left(r\right)+a{\tilde{x}}_C\left(r\right)];$

${\tilde{x}}_0\left(r\right)=\frac{1}{3}[{\tilde{x}}_A\left(r\right)+{\tilde{x}}_B\left(r\right)+{\tilde{x}}_C\left(r\right)];$

The cyclic dispersion Fourier transform calculates A, B, the C three-phase synchronized phasor of gained;

By three-phase phasor synthetic positive and negative, zero sequence phasor;

a＝e ^j2π/3；

3) utilize the error of calculation of bringing owing to frequency departure in 3 each preface phasor of point calibration algorithm correction, modification method is:

A. frequency departure calculates, and uses 3 equally spaced uncorrected synchronized phasors, asks for the deviation of current frequency and rated frequency by 3 point calibration algorithms:

$g(r+m)=\frac{\hat{x}(r+m)}{\hat{x}\left(r\right)};$

$g(r+2m)=\frac{\hat{x}(r+2m)}{\hat{x}(r+m)};$

$f(r+2m)=\frac{g(r+2m)+e^{-j\frac{4\mathrm{\&pi;}}{N}m}/g(r+m)}{2};$

${\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)=f(r+2m)+\sqrt{{\left[f(r+2m)\right]}^2-e^{-j\frac{4\mathrm{\&pi;}}{N}m}};$

Φ(r)＝atan{Re[α ^m(Φ(r))]，Im[α ^m(Φ(r))]}/m；

$\mathrm{\&Delta;f}\left(r\right)=\frac{f_0\mathrm{N\&Phi;}\left(r\right)}{2\mathrm{\&pi;}};$

M: the ratio of sampling rate and phasor correction calculation frequency, promptly every sampling m point is done the phasor corrected Calculation one time;

B. calculate the correction of phasor according to frequency departure, and phasor is proofreaied and correct, bearing calibration is:

Order

$c_2\left(r\right)=\frac{\hat{x}(r+m)-\hat{x}\left(r\right){\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)}{{\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)-{\mathrm{\&alpha;}}^{-m}\left(\mathrm{\&Phi;}\left(r\right)\right)e^{-j\frac{4\mathrm{\&pi;}}{N}m}};$

$c_1\left(r\right)=\hat{x}\left(r\right)+c_2\left(r\right);$

C wherein ₂(r) and c ₁(r) intermediate variable for defining, c ₂(r) the expression phase angle compensation factor, c ₁(r) the middle phasor of revising of the phase angle compensation factor has been taken into account in expression;

$\stackrel{\&OverBar;}{x}\left(r\right)=c_1\left(r\right)\frac{N\sin (\mathrm{\&Phi;}\left(r\right)/2)}{\sin (\mathrm{\&Phi;}\left(r\right)N/2)};$

X (r): revised high-precise synchronization phasor;

4) utilize frequency departure to calculate the phasor frequency, adopt method of difference calculated rate rate of change, computing method are:

f(r)＝f ₀+Δf(r)；

$\frac{\mathrm{df}}{\mathrm{dt}}\left(r\right)={\mathrm{Nf}}_0[f\left(r\right)-f(r-1)];$

Frequency change rate;

5) phasor and frequency measurement are carried out smothing filtering, phasor is carried out the level and smooth of two cycles, and frequency is carried out the level and smooth of one-period, and filtering algorithm is:

${\left|\stackrel{\&OverBar;}{x}\left(r\right)\right|}_{\mathrm{flt}}={\left|\stackrel{\&OverBar;}{x}(r-1)\right|}_{\mathrm{flt}}+\frac{\left|\stackrel{\&OverBar;}{x}\left(r\right)\right|-\left|\stackrel{\&OverBar;}{x}(r-2N)\right|}{2N};$

${\mathrm{\&theta;}\left(r\right)}_{\mathrm{flt}}=\mathrm{\&theta;}{(r-1)}_{\mathrm{flt}}+\frac{\mathrm{\&theta;}\left(r\right)-\mathrm{\&theta;}(r-2N)}{2N};$

$f{\left(r\right)}_{\mathrm{flt}}=f{(r-1)}_{\mathrm{flt}}+\frac{f\left(r\right)-f(r-N)}{N};$

| x (r) | _Flt: the amplitude of the synchronized phasor behind the smothing filtering;

θ (r) _Flt: the phase angle of the synchronized phasor behind the smothing filtering;

F (r) _Flt: the frequency behind the smothing filtering.

Operate in software in industrial computer and the intelligent acquisition card and realized the synchronous phasor measurement method of this compensation cycle discrete Fourier transform (DFT) error.The site test test shows, when sampling rate is 4800Hz, when phasor correction calculation frequency was 200Hz, the stable state accuracy index of algorithm was: the phasor amplitude error is less than 0.5%, and the phasor phase angle error is less than 1 °, and the frequency computation part error is less than 0.01Hz; The dynamic accuracy index of algorithm is: the phasor amplitude error is less than 1%, and the phasor phase angle error is less than 1 °, and the frequency computation part error is less than 0.01Hz.As seen, this synchronous phasor measurement method has higher measuring accuracy, satisfies the requirement of " electric system real-time dynamic monitoring (control) systems technology standard " fully.

Description of drawings

Fig. 1 is for realizing a cover measurement mechanism synoptic diagram of the present invention.

Fig. 2 is for realizing algorithm block diagram of the present invention.

Fig. 3 calculates and phasor correcting process figure for frequency departure.

Embodiment

The synchronous phasor measurement method of the compensation cycle discrete Fourier transform (DFT) error that the present invention puts forward can adopt the multiple hardwares scheme to realize, this example is introduced the measuring system based on industrial computer and DSP capture card that we have realized, comprises relevant hardware configuration and software flow.

The hardware configuration of measuring system as shown in Figure 1, comprise industrial computer, A, B, C threephase potential transformer summation current transformer, and from the A/D capture card on the industrial computer of being integrated in of the secondary side of secondary voltage mutual inductor and secondary current mutual inductor output image data.

At machine end or substation bus bar place, the voltage transformer (VT) summation current transformer is measured bus three-phase voltage and outlet three-phase current respectively, obtain ± the interior ac voltage signal of 120V scope, deliver to secondary voltage mutual inductor and secondary current mutual inductor, convert to ± the interior voltage signal of 5V scope by it; The A/D capture card carries out low-pass filtering, AC sampling and digital-to-analog conversion with secondary voltage mutual inductor and secondary current mutual inductor secondary side voltage signal, and the digital signal that obtains is transferred to industrial computer CPU by pci bus, handle by corresponding software programs, finish the measurement of synchronized phasor and the calculating of frequency, and then output to user's monitoring interface.

Measuring method proposed by the invention mainly is embodied on the software that operates in A/D capture card and the industrial computer, the algorithm structure of software as shown in Figure 2, operate on the industrial computer the flow process that frequency departure calculates and phasor is proofreaied and correct as shown in Figure 3.

Measuring method may further comprise the steps:

1. initialization

A. the sampling number N that points out the user to set every power frequency period (50Hz system), this routine N=96 represents every power frequency period (20 milliseconds) sampling 96 points, thus sampling interval Δ T=20/96 millisecond.

B. point out the relevant configured parameter of user's input measurement device, comprising: (this example is the voltage transformer (VT) no-load voltage ratio

Current transformer ratio (this example is 600A/120V), secondary voltage mutual inductor ratio (this example is 120V/5V), secondary current mutual inductor ratio (this example is 120V/5V), filtering scale-up factor (this example is 0.95) and A/D conversion coefficient (this example is 5V/3FFFh) etc.

C. parameter initialization in the algorithm: the ratio m of sample frequency and phasor emending frequency, this routine m=24 represents that per 24 samplings carry out the phasor correction calculation one time; Each sampling buffer zero clearing:

u _a(t+kΔT-(N-1)ΔT)＝0，u _b(t+kΔT-(N-1)ΔT)＝0，u _c(t+kΔT-(N-1)ΔT)＝0，

i _a(t+kΔT-(N-1)ΔT)＝0，i _b(t+kΔT-(N-1)ΔT)＝0，i _c(t+kΔT-(N-1)ΔT)＝0，k＝0，...，N-1，

They are used for circulating and deposit the synchronous generator output phase voltage u that collects _a, u _b, u _c, current i _a, i _b, i _cSignal forms the sampled point sequence, and wherein t is current sampling instant, and k is the numbering of sample sequence, and that k=N-1 represents is nearest, be t sampled point constantly.

2. under the driving of AC sampling data, carry out following link successively:

A. obtain new sampling number certificate from A/D sampling card, at first they are marked change, that is: consider and the factors such as no-load voltage ratio, filtering scale-up factor and A/D conversion coefficient of PT/CT and secondary PT/CT the digital quantity that obtains is converted to the actual physics value of measuring bus or circuit place correspondence.

After obtaining the actual value of voltage, electric current, the value of voltage/current sample sequence is before pushed away one successively, nearest sampled point leaves on the data point of k=N-1 correspondence.

B. when sampling number is less than N, then continue to wait for new sampled point, otherwise, whenever obtaining a new sampled point, it and preceding N-1 sampled point form a new data window; To the N point data of each new data window, order is carried out following steps and is measured calculating.

1) utilize the cyclic dispersion Fourier transform to ask for and do not consider the synchronized phasor that frequency departure influences, the calculating of voltage and current phasor is identical, therefore adopts uniform way

Express sampled value.

$\hat{x}\left(r\right)=\hat{x}(r-1)+j\frac{\sqrt{2}}{N}[\tilde{x}(r+N-1)-\tilde{x}(r-1)]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="C0314595200131.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}(r-1)}$

2) an order phasor is calculated in every m sampling, utilizes A, B, the synthetic positive sequence phasor of C three-phase phasor, and this routine calculated rate deviation and phasor are proofreaied and correct and all adopted the positive sequence phasor.Synthetic positive sequence phasor deposits initial phasor sequence in, and sequence length is 3, is used for 3 point calibration algorithm computation frequency departures and phasor and proofreaies and correct:

${\hat{x}}_1\left(r\right)=\frac{1}{3}[{\hat{x}}_A\left(r\right)+a{\hat{x}}_B\left(r\right)+a^2{\hat{x}}_C\left(r\right)]$

3) utilize 3 uniformly-spaced positive sequence phasor calculation frequency deviation f:

3.1) use continuous 3 initial positive sequence phasors in the initial phasor sequence, calculate intermediate variable g (r+m), g (r+2m), f (r+2m) by following formula.

$g(r+m)=\frac{\hat{x}(r+m)}{\hat{x}\left(r\right)};$

$g(r+2m)=\frac{\hat{x}(r+2m)}{\hat{x}(r+m)};$

$f(r+2m)=\frac{g(r+2m)+e^{-j\frac{4\mathrm{\&pi;}}{N}m}/g(r+m)}{2};$

3.2) calculating intermediate variable α ^m(Φ (r)) and Φ (r).

${\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)=f(r+2m)+\sqrt{{\left[f(r+2m)\right]}^2-e^{-j\frac{4\mathrm{\&pi;}}{N}m}};$

Φ(r)＝atan{Re[α ^m(Φ(r))]，Im[α ^m(Φ(r))]}/m；

3.3) calculate Δ f according to Φ (r):

$\mathrm{\&Delta;f}\left(r\right)=\frac{f_0\mathrm{N\&Phi;}\left(R\right)}{2\mathrm{\&pi;}};$

4) proofread and correct synchronized phasor according to frequency departure:

4.1) calculating intermediate variable c ₁(r), c ₂(r).

$c_2\left(r\right)=\frac{\hat{x}(r+m)-\hat{x}\left(r\right){\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)}{{\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)-{\mathrm{\&alpha;}}^{-m}\left(\mathrm{\&Phi;}\left(r\right)\right)e^{-j\frac{4\mathrm{\&pi;}}{N}m}};$

$c_1\left(r\right)=\hat{x}\left(r\right)+c_2\left(r\right);$

4.2) the correction synchronized phasor:

$\stackrel{\&OverBar;}{x}\left(r\right)=c_1\left(r\right)\frac{N\sin (\mathrm{\&Phi;}\left(r\right)/2)}{\sin (\mathrm{\&Phi;}\left(r\right)N/2)};$

5) calculate real-time frequency and frequency change rate according to frequency departure:

5.1) calculating real-time frequency f (r).

f(r)＝f ₀+Δf(r)；

5.2) utilize method of difference calculated rate rate of change

$\frac{\mathrm{df}}{\mathrm{dt}}\left(r\right)={\mathrm{Nf}}_0[f\left(r\right)-f(r-1)];$

6) synchronized phasor and frequency are carried out smothing filtering:

${\left|\stackrel{\&OverBar;}{x}\left(r\right)\right|}_{\mathrm{flt}}={\left|\stackrel{\&OverBar;}{x}(r-1)\right|}_{\mathrm{flt}}+\frac{\left|\stackrel{\&OverBar;}{x}\left(r\right)\right|-\left|\stackrel{\&OverBar;}{x}(r-2N)\right|}{2N};$

${\mathrm{\&theta;}\left(r\right)}_{\mathrm{flt}}=\mathrm{\&theta;}{(r-1)}_{\mathrm{flt}}+\frac{\mathrm{\&theta;}\left(r\right)-\mathrm{\&theta;}(r-2N)}{2N};$

${f\left(r\right)}_{\mathrm{flt}}=f{(r-1)}_{\mathrm{flt}}+\frac{f\left(r\right)-f(r-N)}{N};$

Above-mentioned measuring method all can obtain very high precision at power system mesomeric state and transient state process, and can be on the monitor of industrial computer dynamic refresh synchronized phasor and frequency measurement.

Claims (1)

1. the synchronous phasor measurement method of a compensation cycle discrete Fourier transform (DFT) error, it is characterized in that: on fixed sampling interval technique cyclic dispersion Fourier transform basis, when system frequency excursion during as the rated frequency of Fourier transform centre frequency, taken into account the influence of frequency departure, make the measurement of frequency and phasor have higher precision, it contains following steps successively: the 1st step: initialization of variable

Set: N counts the integer greater than 12 for ripple AC sampling weekly;

Δ T is the AC sampling time interval, and initial value is the 20/N millisecond;

M is the ratio of sampling rate and phasor correction calculation frequency, and promptly every sampling m point is done the phasor corrected Calculation one time;

The initial zero clearing of each sampling buffer; The 2nd step: make up data window

With Δ T serves as to gather three-phase voltage or current value with mutual inductor at interval, when sampling number is less than N, waits for new sampled point; When sampling number is more than or equal to N, whenever obtain a new sampled point, it and preceding N-1 sampled point form a new data window; The 3rd step: use the cyclic dispersion fourier transform algorithm, with the phasor of industrial COMPUTER CALCULATION A, B, each phase of C three-phase;

$\hat{x}\left(r\right)=\hat{x}(r-1)+j\frac{\sqrt{2}}{N}[\tilde{x}(r+N-1)-\tilde{x}(r-1)]e^{-j\frac{2\mathrm{\&pi;}}{N}(r-1)};$

Wherein: For just entering the last samples value of data window, $\tilde{x}(r+N-1)=\tilde{x}\left(k\right){|}_{k=r+N-1},$

For just withdrawing from the sampled value of data window, $\tilde{x}(r-1)=\tilde{x}\left(k\right){|}_{k=r-1};$

Wherein: x is the signal effective value;

f ₀Be rated frequency;

Δ f is the deviation of signal bias rated frequency;

is the initial phase angle of signal;

Synchronized phasor for last cyclic dispersion Fourier transform calculates is used for the cyclic dispersion Fourier transformation computation:

The synchronized phasor that calculates for this cyclic dispersion Fourier transform; The 4th step: utilize A, B, the synthetic positive and negative and zero sequence phasor of C three-phase phasor, positive and negative and zero sequence phasor is followed successively by:

${\hat{x}}_1\left(r\right)=\frac{1}{3}[{\hat{x}}_A\left(r\right)+a{\hat{x}}_B\left(r\right)+a^2{\hat{x}}_C\left(r\right)];$

$\widehat{x_2}\left(r\right)=\frac{1}{3}[{\hat{x}}_A\left(r\right)+a^2{\hat{x}}_B\left(r\right)+a{\hat{x}}_C\left(r\right)];$

${\hat{x}}_0\left(r\right)=\frac{1}{3}[{\hat{x}}_A\left(r\right)+{\hat{x}}_B\left(r\right)+{\hat{x}}_C\left(r\right)];$

Wherein,

Be respectively by above-mentioned steps (3) and calculate resulting A, B, C three-phase synchronized phasor;

a＝e ^j2π/3；

The 5th step: utilize the error of calculation of bringing owing to frequency departure in 3 each preface phasor of point calibration algorithm correction;

The 5.1st step: the calculated rate deviation, with 3 equally spaced uncorrected synchronized phasors

Ask for the deviation of current frequency and rated frequency by 3 point calibration algorithms:

Order: $(r+m)=\frac{\hat{x}(r+m)}{\hat{x}\left(r\right)};$

$g(r+2m)=\frac{\hat{x}(r+2m)}{\hat{x}(r+m)};$

$f(r+2m)=\frac{g(r+2m)+e^{-j\frac{4\mathrm{\&pi;}}{N}m}/g(r+m)}{2};$

${\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)=f(r+2m)+\sqrt{{\left[f(r+2m)\right]}^2-e^{-j\frac{4\mathrm{\&pi;}}{N}m}};$

Φ(r)＝αtan{Re[α ^m(Φ(r))]，Im[α ^m(Φ(r))]}/m；

Frequency deviation f (r) is:

$\mathrm{\&Delta;f}\left(r\right)=\frac{f_0\mathrm{N\&Phi;}\left(r\right)}{2\mathrm{\&pi;}};$

The 5.2nd step: calculate the correction of phasor according to frequency departure, and phasor is proofreaied and correct:

Order

$c_2\left(r\right)=\frac{\hat{x}(r+m)-\hat{x}\left(r\right){\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)}{{\mathrm{\&alpha;}}^m\left(\mathrm{\&Phi;}\left(r\right)\right)-{\mathrm{\&alpha;}}^{-m}\left(\mathrm{\&Phi;}\left(r\right)\right)e^{-j\frac{4\mathrm{\&pi;}}{N}m}};$

$c_1\left(r\right)=\hat{x}\left(r\right)+c_2\left(r\right);$

C wherein ₂(r) and c ₁(r) intermediate variable for defining, c ₂(r) the expression phase angle compensation factor, c ₁(r) the middle phasor of revising of the phase angle compensation factor has been taken into account in expression;

Revised high-precise synchronization phasor

$\stackrel{\&OverBar;}{x}\left(r\right)=c_1\left(r\right)\frac{N\sin (\mathrm{\&Phi;}\left(r\right)/2)}{\sin (\mathrm{\&Phi;}\left(r\right)N/2)};$

The 6th step: utilize frequency departure to calculate phasor frequency f (r), adopt method of difference calculated rate rate of change

f(r)＝f ₀+Δf(r)；

$\frac{\mathrm{df}}{\mathrm{dt}}\left(r\right)=N{f}_0[f\left(r\right)-f(r-1)];$

The 7th step: the measured value to phasor and frequency carries out smothing filtering, and phasor is carried out the level and smooth of two cycles, and frequency is carried out the level and smooth of one-period;

The amplitude of the synchronized phasor behind the smothing filtering is | x (r) | _Flt:

${\left|\stackrel{\&OverBar;}{x}\left(r\right)\right|}_{\mathrm{flt}}={\left|\stackrel{\&OverBar;}{x}(r-1)\right|}_{\mathrm{flt}}+\frac{\left|\stackrel{\&OverBar;}{x}\left(r\right)\right|-\left|\stackrel{\&OverBar;}{x}(r-2N)\right|}{2N};$

The phase angle of the synchronized phasor behind the smothing filtering is θ (r) _Flt:

${\mathrm{\&theta;}\left(r\right)}_{\mathrm{flt}}={\mathrm{\&theta;}(r-1)}_{\mathrm{flt}}+\frac{\mathrm{\&theta;}\left(r\right)-\mathrm{\&theta;}(r-2N)}{2N};$

Frequency behind the smothing filtering is f (r) _Flt:

${f\left(r\right)}_{\mathrm{flt}}=f{(r-1)}_{\mathrm{flt}}+\frac{f\left(r\right)-f(r-N)}{N}$ $;$

Images