CN102435844B - Sinusoidal signal phasor calculating method being independent of frequency - Google Patents

Abstract

The invention provides a sinusoidal signal phasor calculating method being independent of frequency. The method comprises the following steps: (1) integrating signal frequency and calculation precision to determine sampling frequency fs, and using a fixed sampling interval Ts=1/fs to obtain sample data of a signal; (2) calculating signal fundamental wave frequency fx (or period), employing a zero-crossing detection method or a Fourier frequency measurement method; (3) according to signal frequency (or period), calculating a largest sampling point number N of the signal needed by a fundamental frequency period, wherein, N=floor (fa/fx), fs is a sampling rate, fx is signal frequency and floor is a function which assumes an integral value of a decimal downwards; (4) selecting signal N point and N+1 point sampling data respectively, and calculating signal N point DFT and N+1 point DFT; (5) by utilizing two DFT calculation results, according to the following two expressions, calculating a real part and an imaginary part at a signal real frequency point fx position, and calculating amplitude of each subharmonic and an instantaneous phase position of the signal. The method is suitable or an occasion with high requirements of precision and real-time property and random change of signal frequency.

Description

A kind of sinusoidal signal phasor calculating method of frequency-independent

Technical field

The invention belongs to digital signal processing and power automation field, be used for calculating the phasor of sinusoidal signal.The present invention is applicable to computational accuracy and requirement of real-time higher, but signal frequency has the occasion of random variation.

Background technology

These three characteristic quantities of the amplitude of signal, phase angle and frequency are important parameters of reflection Operation of Electric Systems characteristic.In time, exactly measuring system frequency and phasor can prognoses system whether stable operation, thereby trigger the safe operation that the system that guarantees is moved in relay protection.And by the time-domain sampling of signal, with DFT, signal is carried out the phasor parameter of the each harmonic that analysis of spectrum just can picked up signal.When electrical network is in power frequency, DFT algorithm based on the fixed sampling interval technique technology has good harmonic characteristic, measurement result is very accurate, but, when mains frequency departs from 50hz, because bringing frequency domain to leak to cause the traditional measurement algorithm, non-synchronous sampling is difficult to meet simultaneously that calculated amount is little, tracking velocity is fast and the high requirement of computational accuracy.

For improving the conventional DFT arithmetic accuracy, multiple improving one's methods proposed at present, revise traditional DFT algorithm, reduce the impact of spectrum leakage.These algorithms are removed when middle step is derived some round-off errors may occur, outside calculated amount is larger, also, because use the data of a plurality of signal periods to cause the algorithm real-time to descend, has limited its application.The present invention proposes a kind of track algorithm based on frequency, effectively simplify the computing of frequency domain interpolation in the past, utilize twice DFT result to carry out linear interpolation, the quick phasor parameter of the each harmonic of picked up signal, can meet simultaneously that calculated amount is little, tracking velocity is fast and the high requirement of computational accuracy.

The more method of research mainly comprises following a few class at present:

1. with windowed function and interpolating method, reduce leakage errors ^{[5] [6] [7] [8]}

In the non-synchronous sampling situation, windowing is blocked to time-domain signal, will cause spectrum leakage.In order to improve computational accuracy, can to sampled data, be weighted by selecting different window functions, thereby signal energy is dropped in main lobe as far as possible, reduce the spectrum leakage, then the sequence after windowing is carried out the FFT computing and carry out double spectral line interpolation obtaining the time and frequency parameter that precision is higher.

2. the synchronization of non-synchronous sampling point ^[4]

Based on the synchronous sampling by software of dual rate or in the situation that the known signal fundamental frequency by original sampled signal is carried out to Lagrange's interpolation, the syncul sequence that obtains being similar to, then apply the DFT algorithm and calculate the every time-frequency parameter of this signal

3. dynamically adjust sampling rate and realize synchronized sampling

Adopt PHASE-LOCKED LOOP PLL TECHNIQUE to realize the hardware synchronization sampling or adjust in real time the sampling interval to signal according to current frequency, then apply the DFT algorithm and calculate the every time-frequency parameter of this signal, can reduce widely the truncation error and the circular error that by non-synchronous sampling, are caused.

Above-mentioned several method can reduce to a certain extent spectrum leakage but weak point is arranged:

1. windowing interpolation method need to be used the data of a plurality of signal periods for guaranteeing precision, and not only operand is large but also real-time and dynamic response performance algorithm are bad

2. these class methods of the synchronization of non-synchronous sampling point need to be carried out the synchronization interpolation to sample sequence in time domain, random variation due to frequency, the operand that becomes interpolator while causing is larger when frequency deviation signal is excessive, the positioning of interpolation point can occur, and at this moment interpolation formula can produce very large error.And algorithm usings signal zero crossing as the starting point of ideal, when a plurality of zero crossing of the excessive appearance of signal harmonic, the interpolation mistake will appear.

3. due to system frequency random variation within the specific limits, adopt equal angular strategy to need dynamically to adjust sampling interval, this can cause chain reaction and cause system unstable concerning the synchronous sampling system that involves a plurality of equipment, traditional digital filter also can't be suitable for.And in actual conditions, the speed of A/D sampling can't meet synchronized sampling so exactly, especially when signal is interfered the occurrence frequency fluctuation.

List of references

1. Liu builds woods, masses troops a kind of the 2011st the 4th phase of volume of high precision simple signal frequency estimation algorithm radio engineering

2. the self-adaptation such as the outstanding Han Ying's tongued bell of Min Yong Ding Ren is adjusted phasor on-line measurement algorithm research Power System and its Automation journal in the April, 2010 of sampling rate

3. Wu Jun Wu Chong sky is based on improvement Fu Shi Frequency Measurement Algorithm Jiangsu electrical engineering in the July, 2008 of Kaiser window

4. more than Zeng Zehao, there is spirit to be permitted the dimension victory based on the spectrum leakage analysis of interpolation synchronized algorithm and analogue system emulation technology in October, 2005

5. yellow pure, windowed interpolation improvement algorithm Proceedings of the CSEE in the August, 2005 of river subgroup frequency analysis

6. Li Xinfu Liu believer Cui Yulong electrical equipment algorithm of harmonics analysis Hebei University of Technology journal in October, 2003 of detecting

7. Pang Hao, Li Dongxia, ancient sacrificial utensil sky, wait application FFT to carry out the improvement algorithm [J] of harmonic analysis in power system. Proceedings of the CSEE 2003,23

8. Liang Zhi state, Sun Yu. a kind of digital measuring method [J] of signal period. Chinese journal of scientific instrument, 2003,24 (4, supplementary issue): 195-198

9. high cloud roc Teng Zhao victory minister in ancient times Bai Yuan is based on harmonic analysis method electronic letters, vol in the Dec, 2000 of Kaiser window double spectral line interpolation FFT

10. Yang Hong ploughs favour brocade Hou Peng Harmonious Waves in Power Systems and harmonic detection method summary Automation of Electric Systems in October, 1998

Summary of the invention

The objective of the invention is to propose a kind of real-time phasor calculating method that is not subjected to the electric system sinusoidal signal of frequency influence.The present invention proposes a kind of track algorithm based on frequency, effectively simplify the computing of frequency domain interpolation in the past, utilize twice DFT result to carry out linear interpolation, the quick phasor parameter of the each harmonic of picked up signal, can meet simultaneously that calculated amount is little, tracking velocity is fast and the high requirement of computational accuracy.

The present invention realizes by such scheme: a kind of phasor calculating method of electric system sinusoidal signal of frequency-independent,

1). integrated signal frequency and computational accuracy are determined sample frequency f _s, and with the fixed sample interval T _s=1/f _sObtain the sample data of signal;

2). calculate signal fundamental frequency f _x(or cycle), can adopt zero passage detection method or Fu Shi Measuring Frequency Method ^[3]Etc. method ^{[1] [8]}

3). according to signal frequency (or cycle), calculate the required maximum sampling number N of a fundamental frequency cycles of signal;

N=floor (f _s/ f _x) f _s-sample rate f _xThe function that-signal frequency floor-rounds decimal downwards

4). choose respectively signal N point and N+1 point sampling data, calculate signal N point DFT and N+1 point DFT;

5). utilize twice DFT result of calculation to calculate the true frequency f of signal by following two kinds of expression formulas _xReal part and the imaginary part at place;

$y_x=\frac{y_2\·(f_x-f_1)+y_1\·(f_2-f_x)}{f_2-f_1}$

Or

$y_x=y_2\·\frac{N}{N_x}\·(N+1-N_x)+y_1\·\frac{N+1}{N_x}(N_x-N)$

N _x＝f _s/f _x N＝floor(f _s/f _x)

In formula: N _x-per cycle sampling number, the integer sampling number of the per cycle maximum of N-

F ₂Gained fundamental frequency (f when the corresponding N of employing point DFT calculates ₂=1/NT _s)

F ₁Gained fundamental frequency (f when the corresponding N+1 of employing point DFT calculates ₁=1/[(N+1) T]))

F _s-sample rate f _xThe function that-signal frequency floor-rounds decimal downwards

Y ₁Real part or the imaginary part of first-harmonic or harmonic wave in-respective signal N point DFT result

Y ₂Real part or the imaginary part result of calculation of first-harmonic or harmonic wave in-respective signal N+1 point DFT result

Y _xReal part or the imaginary part of the final first-harmonic of-interpolation gained or harmonic wave

6). calculate amplitude and the instantaneous phase of signal each harmonic.

The invention has the beneficial effects as follows: the collection of the real-time phasor of electric system sinusoidal signal is parameter important and basic in power grid operation, and it is very important needing the real-time phasor of accurate acquisition fast,

1. the present invention calculates the sampled data only need one-period, and the phasor value that (signal period) can accurate picked up signal within the shortest time is better than similar algorithm dynamic response characteristic;

2. algorithm is simple, can directly use DFT or by after the sampled point zero padding, utilizing fft algorithm;

3. memory consumption is few, and the sampled data sequence only needs to store the data of one-period left and right;

4. by increasing sampling rate (reducing sampling interval), can increase computational accuracy but not affect dynamic perfromance;

5. result of calculation is not subjected to the impact of signal frequency

6. the employing equal interval sampling, do not need dynamically to change sampling interval with the adaptation signal frequency change, can use traditional digital filtering technique to guarantee simultaneously the stable of sampling system.

The accompanying drawing explanation

Fig. 1 is signal phasor interpolation schematic diagram of the present invention

Fig. 2 is the real part/imaginary part interpolation schematic diagram of signal phasor of the present invention

Table 1 is the error analysis table of comparisons of the present invention, and in table, sampling rate unit is: sampled point/second, error is relative error.

Specific implementation

1. the ultimate principle of algorithm

According to digital signal processing theory, the N of signal point DFT conversion is f in frequency domain resolution _s=1/ (NT _s), the one-period that sequence length is signal if get, its fundamental frequency that corresponds to signal is f _Xh=1/ (NT _s).In like manner, N+1 point DFT conversion is f in frequency domain resolution _s=1/[(N+1) T]), the respective signal frequency is f _Xl=1/[(N+1) T]).For fundamental frequency, be f _xSignal, can find suitable sequence length N, makes following formula set up:

F _Xl<f _x<f _XhIn formula: f _Xh=1/NT _s, f _Xl=1/[(N+1) T _s]

When the frequency of signal at [f _Xl, f _Xh] while changing in interval, the phasor value of signal (and corresponding real part/imaginary part) forms one with frequency f _xFor the continuous curve of variable, frequency f _xThe phasor at place can be passed through at frequency separation [f _Xl, f _Xh] linear interpolation obtains, and closed interval [f _Xl, f _Xh] large young pathbreaker determine the error of interpolation.

2. signal frequency is calculated

The algorithm that the present invention proposes need to first record out signal frequency, and the simplest method is carry out the picked up signal cycle and then extrapolate signal frequency by the position that checks two zero crossings in burst, and its ultimate principle is as follows:

(1). search signal sequence x (n), if find that there is x (n)≤0 and x (n+1)>=0, at sampling instant t _nAnd t _N+1Between have a zero crossing;

(2). on time shaft, by Lagrange's interpolation, obtain the moment t of current zero crossing _cur(t _n≤ t _cur≤ t _N+1);

(3). by calculating current zero point and previous zero point t _preBetween time interval, can obtain signal period T _x=t _cur-t _pre

(4). record current null position and go to step 1 and carry out new search at zero point

The impact that is noise reduction and harmonic wave when carrying out searching for zero point can be carried out bandpass filtering to sequence.The calculating of signal frequency can certainly utilize document ^{[1] [3] [8]}The method of middle proposition, as the Fu Shi Measuring Frequency Method.

3. maximum cycle sampling number N

If the signal fundamental frequency cycles is T _x(respective frequencies f _x=1/T _x), sampling interval is T _s(corresponding sample rate f _s=1/T _s), the sampled data that obtains of per cycle count into:

N _x＝T _x/T _s＝f _s/f _x (1)

N when meeting the synchronized sampling condition _xFor integer, do not leak.But due to the signal frequency disturbance, will cause sampling not meet weekly a phase integer point (is N _xBe not integer) requirement, if still with integral point, calculate the near frequency domain of DFT, leak.

Algorithm need to be determined each fundamental frequency cycles maximum integer sampling number N, calculates in the following way:

N＝floor(T _x/T _s)＝floor(f _s/f _x) (2)

In following formula, function f loor (x) means to round downwards, and obviously following formula is set up:

N≤N _x≤N+1

(3)

4.DFT calculate

The present invention uses N sampled data of rectangular window intercept signal sampled value sequence to calculate simply, according to the definition of discrete Fourier transform (DFT) (DFT), can directly by following expression formula, calculate signal N point DFT

$\begin{array}{ccc}X\left(k\right)=\frac{2}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}x\left(n\right)\&CenterDot;w\left(k\right)& w\left(k\right)=e^{-j2\mathrm{k\&pi;}/N}& 0\&le;k\&le;N-1\end{array}---\left(4\right)$

When the few operand of overtone order of needs analysis is little, can directly utilize following formula to calculate.Similarly method can computational length N+1 point DFT.When the overtone order of analyzing when needs is many, can by the N sequence length after, mend 0, make to mend length after 0 and reach certain power of 2 and re-use FFT calculating.

5. linear interpolation

The present invention adopts the method for linear interpolation to obtain frequency f _xThe signal phasor at place.If only to the amplitude of signal and phase angle is interested can be by accompanying drawing 1, directly interpolation obtains frequency f _xAmplitude and the angle at place.If need simultaneously signal frequency point f _xPlace's real part and imaginary part can be by accompanying drawings 2, to real part and the imaginary part difference interpolation of phasor.

Frequency f _xThe y of the phasor at place _xValue can obtain by Lagrange interpolation polynomial, that is:

$y_x=\frac{y_2\&CenterDot;(f_x-f_1)+y_1\&CenterDot;(f_2-f_x)}{f_2-f_1}---\left(5\right)$

In formula: f ₂The corresponding f of fundamental frequency when the corresponding N of employing point DFT calculates _Xh, f ₁The corresponding f of fundamental frequency when the corresponding N+1 of employing point DFT calculates _Xl.

The interpolation formula that is converted into sampling number is

$y_x=\frac{y_2\&CenterDot;(\frac{1}{N_x{T}_s}-\frac{1}{(N+1)T_s})+y_1\&CenterDot;(\frac{1}{{\mathrm{NT}}_s}-\frac{1}{N_x{T}_s})}{\frac{1}{{\mathrm{NT}}_s}-\frac{1}{(N+1)T_s}}$ T _s--sampling interval (6)

The following formula abbreviation is:

$y_x=y_2\&CenterDot;\frac{N}{N_x}\&CenterDot;(N+1-N_x)+y_1\&CenterDot;\frac{N+1}{N_x}(N_x-N)---\left(7\right)$

Formula (7) shows, when meeting the synchronized sampling condition, N to be arranged _x=N+1 (correspondence: f _x=f ₁, y _x=y ₁) or N _x=N (correspondence: f _x=f ₂, y _x=y ₂), while namely meeting the synchronized sampling condition, algorithm is equal to conventional DFT.

According to formula (4), the k subharmonic of signal is corresponding to the k root spectral line of DFT result, so can obtain corresponding to kf in the lump by the method for Frequency domain interpolation _xThe signal each harmonic parameter at place.

6. amplitude and phase calculation

According to formula (4), by DFT result of calculation, be complex vector, i.e. the real part of phasor and imaginary part, and along with the corresponding k subharmonic of the value Different Results of k.For obtaining signal at frequency f _xThe phasor at place, can carry out respectively linear interpolation by formula (7) to real part and the imaginary part of DFT result, then from the result of interpolation, utilizing the methods such as plural delivery to obtain amplitude and phase angle.

7. applicating example

This sentences the civil power AC signal is example, can be different to the identical just sampling rate of other frequency signal operation stepss.If sample rate f _s=3200 times/second (sampling interval is T _s=312.5 μ S are 64 per cycles of sampling point during signal frequency 50hz), practical operation step of the present invention is as follows:

(1). with sampling interval T _sSignal is carried out to equal interval sampling

(2). utilize 2 given methods to calculate the current frequency of signal, suppose that herein the gained frequency is f _x=49.5hz;

(3). maximum integer sampled point number of per cycle of picked up signal and required sequence length f _s/ f _x=3200/49.5=64.6464, round to obtain N=floor (f downwards _s/ f _x)=64

(4). from current sampling point forward, get respectively and 65 data at 64 and build data windows, and calculate and DFT at 64 at 65;

(5). by formula (7), calculate real part and the imaginary part of signal fundametal compoment according to the result of step 4, and calculate thus amplitude and phase angle;

(6). similarly method is calculated needed harmonic wave real part and imaginary part, and derives amplitude and the phase angle of harmonic wave.

8. error analysis and countermeasure

From Fig. 1 and Fig. 2, finding out, work as f ₁And f ₂Between the less error of spacing less, and f ₁And f ₂Relevant with sequence length.Increase per cycle sampling number (correspondence reduces sampling interval), can reduce f ₁And f ₂Between spacing, cause improving precision, but increase sampling number, can increase operand and storage overhead, during application, should do balance according to precision and system arithmetic capability.

In the occasion that actual ADC sampling rate can't increase, can before calculating, use the Sinc method of interpolation to carry out the time domain interpolation, to between every two points of crude sampling sequence, carrying out interpolation, increase sampling rate, its result is equivalent to increase per cycle sampling number, thereby reach the raising precision, the method does not have the synchronized problem of the described non-synchronous sampling point of preamble.

Utilize Matlab to carry out error scanning to this algorithm, during analysis, amplitude normalization is 1 (now obtaining the amplitude error relative error), take frequency and initial phase angle to be variable, sampling interval T _sFor parameter builds and to contain first-harmonic and 2,3,5 harmonic sine sequences by following formula, recycle the amplitude of above-mentioned algorithm calculating signal and phase place and theoretical value and compare error.

x(n)＝sin(2πf _xnT _s+φ ₁)+sin(4πf _xnT _s+φ ₂)+sin(6πf _xnT _s+φ ₃)+sin(10πf _xnT _s+φ ₅)

In formula: T _s-sampling interval, f _x-signal frequency, φ ₁, φ ₂, φ ₃, φ ₅-harmonic wave initial phase angle

Table 1 provides the error information under electric system several sampling rates commonly used; even showing under the condition of sampling rate 3200 times/second, result still can meet the accuracy requirement of electric system to measure and control device; due to algorithm, only use the data in single cycle of signal, so the dynamic response real-time also can finely meet the requirement of protective device.

The 45hz-55hz sinusoidal signal algorithm simulating error table of comparisons under the different sampling rates of table 1

Sampling rate	3200	4000	6400	8000	10000	12800
							Phase error	0.001149	0.000728	0.000291	0.0001858	0.0001189	0.0000715
Range error	0.0011	0.0006959	0.0003	0.0001769	0.000113151	0

Claims (2)

1. the sinusoidal signal phasor calculating method of a frequency-independent is characterized in that step is as follows:

1). determine sample frequency f _s, with interval T _s=1/f _sObtain the sample data of signal;

2). calculate sinusoidal signal frequency f _xOr cycle T _x

3). according to sinusoidal signal frequency or computation of Period, go out the required maximum sampling number N of half period of signal;

4). get respectively sinusoidal signal N point and N+1 point sampling data, calculate N point half-wave DFT and N+1 point half-wave DFT;

5). utilize twice half-wave DFT result of calculation by the true frequency f of interpolation calculation signal _xReal part and the imaginary part at place;

6). according to amplitude and the instantaneous phase of real part and imaginary part calculating signal;

The maximum integer sampled point number N of per semiperiod of sinusoidal signal calculates by following method:

N=floor (f _s/ 2f _x)=floor (T _x/ 2T _s) wherein: f _s-sampling rate, f _x-signal frequency, T _x-signal actual cycle, T _s-sampling interval, the function that floor-rounds decimal downwards.

2. the sinusoidal signal phasor half-wave computing method of frequency-independent according to claim 1 is characterized in that: the half-wave DFT according to sinusoidal signal N point and N+1 are ordered, utilize expression to calculate the true frequency f of signal _xReal part and the imaginary part at place:

$X_r^{T_x}=X_r^{(N+1)T_s}\&CenterDot;V_{T_x}+X_r^{N{T}_s}(1-V_{T_x})$

$X_i^{T_x}=X_i^{(N+1)T_s}\&CenterDot;V_{T_x}+X_i^{N{T}_s}(1-V_{\mathrm{Tx}})$

$V_{T_x}=\frac{T_x}{2T_s}-N$

N＝floor(f _s/f _x)＝floor(T _x/2T _s)

In formula: T _x-signal actual cycle, T _s-sampling interval,

-corresponding the cycle is T _xSignal half-wave algorithm linear interpolation weights,

The integer sampling number of the per semiperiod maximum of N-,

Gained real part when the-corresponding N of employing point half-wave DFT calculates,

Gained imaginary part when the-corresponding N of employing point half-wave DFT calculates,

Gained real part when the-corresponding N of employing point half-wave DFT calculates,

Gained imaginary part when the-corresponding N of employing point half-wave DFT calculates,

-corresponding the cycle is T _xSignal half-wave DFT gained real part actual value while calculating,

-corresponding the cycle is T _xSignal half-wave DFT gained imaginary part actual value while calculating.

Images