A kind of sinusoidal signal phasor half-wave computing method of frequency-independent
One, 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.
Two, 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.At present, the principle of Most protection is based on the stable state fundamental component after fault, how fast and accurately Fundamental-frequency Current, voltage to be estimated the subject matter that Microcomputer Protection Calculation Method faces from the fault transient signal.Usually, the length of data window is depended in the time-delay of calculating, and the shorter algorithm response time of data window is faster.Although all integer harmonic components of all-wave fourier algorithm energy filtering, and good stability, its data window needs 1 cycle at least, so response speed is slower; And Half-cycle Fourier Algorithmfor is only used the sampled data of half period, and response is fast, the waveform recognition of widespread use and transformer differential quick-break, transformer excitation flow
[10]In the QA electric system Microcomputer Protection Calculation Method of needs.In addition on the one hand, when electrical network is in power frequency, by the time domain equal interval sampling to signal, traditional half-wave DFT algorithm fast picked up signal the phasor parameter, 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 satisfy simultaneously that calculated amount is little, tracking velocity is fast and the high requirement of computational accuracy.
The research emphasis of half-wave arithmetic accuracy is to eliminate the impact of even-order harmonic and aperiodic component at present, needs is improved the occasion of the frequency domain leakage errors that brings because of non-synchronous sampling, mainly adopts the all-wave fourier transform algorithm.Now existing multiple improving one's methods is used for revising traditional DFT algorithm, reduces the impact of spectrum leakage.These algorithms are removed when the step of centre 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, and the inapplicable half-wave algorithm of these algorithms.The present invention proposes a kind of track algorithm based on frequency, utilize twice half-wave DFT result to carry out linear interpolation, the quick phasor parameter of picked up signal can satisfy 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. reduce leakage errors with windowed function and interpolating method
[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 be weighted sampled data by selecting different window functions, reduce the spectrum leakage in main lobe thereby signal energy is dropped on as far as possible, then the sequence after windowing is carried out the FFT computing and carry out double spectral line interpolation obtaining the higher time and frequency parameter of precision.
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 Lagrange's interpolation, then the syncul sequence that obtains being similar to is used the DFT algorithm and is calculated 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 hardware synchronization sampling or the sampling interval adjusted signal real-time according to current frequency, then use the DFT algorithm and calculate the every time-frequency parameter of this signal, can reduce widely the truncation error and the circular error that are caused by non-synchronous sampling.
4. improvement half-wave Fourier algorithm
[11] [12] [13]
Traditional half-wave algorithm is improved, and the fast filtering algorithm by difference low pass or employing constrictive band-pass filter combine with the half-wave Fourier suppresses the impact of low frequency component and higher hamonic wave, improves the half-wave arithmetic accuracy.
Above-mentioned several method respectively has deficiency, and does not all make a concrete analysis of frequency change to the impact of half-wave algorithm.Method 1 is all a plurality of cycles for signal, can reduce to a certain extent spectrum leakage but weak point is arranged, not only operand is large but also real-time and dynamic response performance algorithm are bad, method 2 need to be carried out the synchronization interpolation to sample sequence in time domain, random variation due to frequency, the operand that becomes interpolator when 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.Method 3 adopts equal angular strategy to need dynamically to adjust sampling interval, though can solve the leakage problem that the half-wave algorithm runs into, 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 satisfy synchronized sampling so exactly.Method 4 mainly solves low frequency component to the impact of half-wave algorithm, leaks and pay close attention to frequency domain, and the scene when not considering that the frequency jitter scope surpasses the bandwidth of wave filter.
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. 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 more than Ceng Zehao
5. yellow pure, windowed interpolation improvement algorithm Proceedings of the CSEE in the August, 2005 of river subgroup frequency analysis
6. algorithm of harmonics analysis Hebei University of Technology's journal lO month in 2003 of detecting of Li Xinfu Liu believer Cui Yulong electrical equipment
7. Pang Hao, Li Dongxia, ancient sacrificial utensil sky wait and use the improvement algorithm [J] that FFT carries out 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 the harmonic analysis method electronic letters, vol 2OOO Dec of Kaiser window double spectral line interpolation FFT
10. Han Zheng celebrates the refined Lie group Zhan of Grolsch based on excitation flow recognition method Automation of Electric Systems in the July, 2005 of half-wave Fourier algorithm
11. Li Bin, Li Yongli, the 21 Annual Conference collection of thesis of the improved half-wave Fourier filtering of He Jiali algorithm Power System and its Automation specialty
12. military Xiao of fourth book literary composition Zhang Chengxue Gong Qing meets improvement Automation of Electric Systems in the March, 1999 of first Half-cycle Fourier Algorithmfor
13. Zou's wisdom, Li Xiaocong, Luo Xiaofen, the improvement half-wave Fourier algorithm protecting electrical power system of the effective filtering even-order harmonic of Fan Chuizheng with control on October 16th, 2009
Three, summary of the invention
The objective of the invention is to propose a kind of real-time phasor half-wave computing method that are not subjected to the sinusoidal signal of frequency influence.The present invention proposes a kind of track algorithm based on frequency, utilize twice half-wave DFT result to carry out linear interpolation, the quick phasor parameter of picked up signal can satisfy 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 half-wave computing method of 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 frequency f
x(or cycle T
x), can adopt zero passage detection method or Fu Shi Measuring Frequency Method
[3]Etc. method
[1] [8]
3). calculate the required maximum sampling number N of half period of signal according to signal frequency (or cycle);
N=floor (f
s/ 2f
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, press half-wave DFT algorithm and calculate N point DFT and N+1 point DFT;
5). utilize twice half-wave DFT result of calculation to calculate the true frequency f of signal by following expression formula
xReal part and the imaginary part at place;
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 when calculating
-corresponding the cycle is T
xSignal half-wave DFT gained imaginary part actual value when calculating
6). calculate amplitude and the instantaneous phase of signal.
The invention has the beneficial effects as follows:
1. calculate the sampled data that only needs half period, get final product the phasor value of accurate picked up signal, better than similar algorithm dynamic response characteristic;
2. result of calculation is not subjected to the impact of signal frequency;
3. memory consumption is few, and the sampled data sequence only needs to store the data of half period left and right;
4. 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;
5. by increasing sampling rate (reducing sampling interval), can increase computational accuracy but not affect dynamic perfromance;
Be used for calculating the phasor of sinusoidal signal, be applicable to computational accuracy and requirement of real-time higher, but the occasion of signal frequency random variation.The present invention is at first by equal interval sampling picked up signal sample data, then calculate the fundamental frequency of signal and the required maximum sampling number N of half period of definite signal according to sampled data, then press half-wave DFT algorithm and calculate signal N point DFT and N+1 point DFT, recycle twice DFT result of calculation and carry out linear interpolation and obtain real part and imaginary part on the true frequency of signal, calculate at last signal amplitude and instantaneous phase.
Four, description of drawings
Fig. 1 is signal non-synchronous sampling 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.
Five, specific implementation
A kind of phasor half-wave computing method of sinusoidal signal of frequency-independent, concrete methods of realizing is as follows.
1. the theory of algorithm is derived
1.1. traditional half-wave algorithm
To be the symmetry of utilizing sinusoidal signal develop from the basis of all-wave fourier transform algorithm tradition half-wave fourier transform algorithm, and the actual cycle of establishing signal is T (the corresponding semiperiod is T/2).
In formula:
--half-wave Fourier algorithm gained real part when the signal period is T
--half-wave Fourier algorithm gained imaginary part when the signal period is T
1.2. half-wave interpolation algorithm principle
If the actual cycle of signal is T, certainly exist the following formula establishment that Integer N makes
NT
s≤T/2≤(N+1)T
s (3)
By
formula 1 and
formula 2, when semiperiod of signal at [NT
s, (N+1) T
s] between when changing continuously, the real part of half-wave algorithm
And imaginary part
To consist of one section continuous curve, real part
Two end points of curve are respectively
Imaginary part
Two end points of curve are respectively
According to the Lagrange interpolation formula theorem, in error allowed band situation, the actual value of half-wave algorithm real part and imaginary part can be respectively at closed interval [NT
s, (N+1) T
s] obtain amplitude and phase angle based on the half-wave new algorithm by linear interpolation.
2. signal frequency is calculated
The method that the present invention proposes need to first record out signal frequency, and the simplest method is come 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 x (n)≤0 and x (n+1) 〉=0 are arranged, at sampling instant t
nAnd t
n+1Between have a zero crossing;
(2). by Lagrange's interpolation, obtain the moment t of current zero crossing on time shaft
cur(t
n≤ t
cur≤ t
n+1);
(3). by calculating current zero point and previous zero point t
preBetween time interval, can get signal period T
x=t
cur-t
pre
(4). record current null position and turn 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 is as the Fu Shi Measuring Frequency Method.Because the frequency of electrical network is relatively stable, frequency change is little within a signal period, can be with the measured value in the previous cycle foundation as current calculating.
3. maximum half cycle sampling number N calculates
Be T to fundamental frequency cycles
xSignal (respective frequencies f
x=1/T
x), if sampling interval is T
s(corresponding sample rate f
s=1/T
s), the required data window of half-wave algorithm is counted:
N
x=(T
x/2)/T
s=f
s/2f
x (8)
N when satisfying the synchronized sampling condition
xBe integer, do not leak.But due to the signal frequency disturbance, will cause sampling not satisfy a semiperiod integer point (is N
xBe not integer) requirement, will cause frequency domain to leak if still calculate half-wave DFT with integral point.Accompanying drawing 1 is non-synchronous sampling schematic diagram in the signal half period, in this case, in half period still to be spaced apart T
sSampling can't collect an integer sampling point.
According to formula 3, formula 8 result of calculations are rounded downwards can obtain the required definite semiperiod maximum integer sampling number N of algorithm:
N=floor(T
x/2T
s)=floor(f
s/2f
x) (9)
4. linear interpolation
The present invention uses N sampled data of rectangular window intercept signal sampled value sequence to calculate simply, definition according to discrete half-wave Fourier transform (DFT), with formula 1 and formula 2 discretizes, use discrete summation instead of integration, can obtain respectively real part and the imaginary part of N point half-wave DFT.
--real part (10) during N point data window
--imaginary part (11) during N point data window
Similarly, use the data of N+1 point of signal, estimate that real part and imaginary part can get following result
Real part (12) during-
N+1 point data window
Imaginary part (13) during-
N+1 point data window
According to Lagrange interpolation formula, in the situation that error allows, (fundamental frequency cycles is T to signal
x) real part and the imaginary part actual value of half-wave algorithm can obtain by linear interpolation respectively, namely
Order
(according to formula 3,
) (15)
Utilize the discretize result of formula 10 and formula 12, formula 15 substitution formulas 14 had:
In like manner can get the imaginary part interpolation
The real part that calculates according to formula 16 and formula 17 and imaginary part can further utilize the method such as plural delivery to obtain amplitude and phase angle, and closed interval [NT
s, (N+1) T
s] large young pathbreaker determine the error of interpolation, accompanying drawing 2 is the interpolation schematic diagram of real part imaginary part.
Must have when satisfying the synchronized sampling condition
Have this moment
When namely satisfying the synchronized sampling condition, new algorithm is equal to traditional half-wave DFT.
5. 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 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 semiperiod of picked up signal and required sequence length
f
s/ 2f
x=3200/99=32.323232 rounds to get N=floor (f downwards
s/ 2f
x)=32
(4). from current sampling point forward, get respectively and 33 data at 32 and build data windows, and calculate real part and the imaginary part of and half-wave DFT at 32 at 33;
(5). calculate signal half-wave algorithm real part and imaginary part actual value by formula (16) and (17) according to the result of step 4, and calculate thus amplitude and phase angle;
6. error analysis and countermeasure
Only use the sampled data of half period due to Half-cycle Fourier Algorithmfor, although response is fast, filter capacity relatively a little less than, if when signal contains even-order harmonic component or aperiodic component, algorithm will produce larger error, about the method for eliminating this impact referring to document
[11] [12] [13]
According to Lagrange interpolation formula about the analysis of error term and in conjunction with can find out in Fig. 1 and Fig. 2, sampling interval T
sLess error is less.Obviously increase per cycle sampling number (correspondence reduces sampling interval), can improve precision, can increase operand and storage overhead but increase sampling number, should do balance according to precision and system's arithmetic capability during application.
Can use the Sinc method of interpolation to carry out the time domain interpolation before calculating in the occasion that actual ADC sampling rate can't increase, increase sampling rate to carrying out interpolation between every two points of crude sampling sequence, its result is equivalent to increase per cycle sampling number, thereby reaches the raising precision.
Utilize Matlab to carry out error scanning to this algorithm, during analysis, amplitude normalization is 1 (obtaining the amplitude error relative error this moment), take frequency and initial phase angle as variable, and sampling interval T
sPress following formula for parameter and build sinusoidal sequence, 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+φ
1)
In formula: T
s-sampling interval, f
x-signal frequency, φ
1-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 satisfy electric system to the accuracy requirement of measure and control device; only use the data of signal half period due to algorithm, so the dynamic response real-time also can finely satisfy 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.000001 |
Although the present invention discloses as above with preferred embodiment, so it is not to limit the present invention.The persond having ordinary knowledge in the technical field of the present invention, without departing from the spirit and scope of the present invention, when being used for a variety of modifications and variations.Therefore, protection scope of the present invention is as the criterion when looking claims person of defining.