CN108037361B - High-precision harmonic parameter estimation method based on sliding window DFT - Google Patents
High-precision harmonic parameter estimation method based on sliding window DFT Download PDFInfo
- Publication number
- CN108037361B CN108037361B CN201711265070.XA CN201711265070A CN108037361B CN 108037361 B CN108037361 B CN 108037361B CN 201711265070 A CN201711265070 A CN 201711265070A CN 108037361 B CN108037361 B CN 108037361B
- Authority
- CN
- China
- Prior art keywords
- frequency
- signal
- estimation
- harmonic
- component
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 23
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 46
- 238000004458 analytical method Methods 0.000 claims abstract description 10
- 230000003595 spectral effect Effects 0.000 claims description 16
- 230000009466 transformation Effects 0.000 claims description 11
- 238000005070 sampling Methods 0.000 claims description 9
- 239000000654 additive Substances 0.000 claims description 3
- 230000000996 additive effect Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 abstract description 7
- 238000004364 calculation method Methods 0.000 description 7
- 230000008901 benefit Effects 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 239000011159 matrix material Substances 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000004888 barrier function Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 239000012456 homogeneous solution Substances 0.000 description 1
- 238000009776 industrial production Methods 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 239000004745 nonwoven fabric Substances 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000013021 overheating Methods 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/16—Spectrum analysis; Fourier analysis
Landscapes
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses a high-precision harmonic parameter estimation method based on sliding window DFT, which is suitable for estimating harmonic components of an electric power system, wherein after the electric power system contains harmonic components and power grid signals are subjected to SWDFT, the conversion result has the same linear relation with the original time domain signals, namely the converted sequence still retains the mathematical structure of multi-frequency signals, and the parameter estimation can be carried out by adopting the conventional harmonic analysis model. In addition, the SWDFT also enhances the signal-to-noise ratio of each component, and inhibits the interference of noise on the estimation performance. Therefore, the method firstly carries out SWDFT conversion on the harmonic signals in the time domain, and then adopts the Prony algorithm on the converted result to obtain the estimated values of the frequency, the amplitude and the phase of each harmonic component of the signals. In order to enhance the estimation performance and robustness of the algorithm, the algorithm introduces a complex least square method criterion for the original Prony algorithm for expansion, and the precision of parameter estimation is effectively improved.
Description
Technical Field
The invention relates to a harmonic parameter estimation method for stationary periodic signals, in particular to a harmonic parameter high-precision estimation method based on sliding window Fourier transform (SWDFT) and Prony algorithm, which is used for estimating harmonic frequency, amplitude and phase parameters and belongs to the technical field of signal harmonic analysis of power systems.
Background
With the development of power electronic technology and the expansion of industrial production, electronic components and nonlinear loads in power systems are increasing. The nonlinear load also injects a large amount of higher harmonics into the power grid while bringing huge economic benefits, harmonic pollution becomes more and more serious, and various damages are caused to a power system and users while the distortion degree of an electric signal is intensified, such as motor overheating, abnormal relay protection of the power system, increased loss of a power transmission line and the like. At present, measures for injecting reverse components of corresponding harmonic components into a system and methods for installing a filter and the like are mostly adopted for harmonic control, so that harmonic analysis is carried out on a power grid signal, and great research significance and practical value are provided for electric energy metering and electric energy quality analysis and control.
Discrete Fourier Transform (DFT) is the most common method for harmonic analysis in a power system, and the algorithm has the characteristics of simple operation and high calculation efficiency, and has a better estimation result on stable and noiseless signals. But the inherent defects of frequency aliasing, frequency spectrum leakage, barrier effect and the like exist, and the accuracy of harmonic analysis is influenced. The improved measures such as windowing and interpolation can reduce estimation errors caused by frequency spectrum leakage and barrier effect to a certain extent, but when the frequency interval of each harmonic component is close, the method cannot work effectively. In order to improve the estimation accuracy, the iterative DFT algorithm is developed, but the introduction of iteration generates additional calculation amount, which causes huge calculation burden. Besides DFT, there are some time-domain-based harmonic analysis methods, such as kalman filtering method, however, the implementation of these methods needs to be based on the knowledge of some characteristics of the signal data to be processed, and the state matrix needs to be accurately defined in advance, otherwise, the estimation effect in advance will not be achieved.
Disclosure of Invention
The purpose of the invention is as follows: aiming at the prior art, a high-precision harmonic parameter estimation method based on sliding window DFT is provided, which can greatly improve the estimation precision of parameters such as power grid signal frequency, amplitude and phase.
The technical scheme is as follows: a high-precision harmonic parameter estimation method based on sliding window DFT comprises the following steps:
step 1: collecting multi-frequency signals in an electric power system, and obtaining discrete signals x (n) at time n after sampling
Wherein a ism，ωm，φmRespectively showing the amplitude, digital frequency and phase of the mth harmonic, wherein M is the number of harmonics contained in the signal, s (n) is a multi-frequency harmonic signal part, q (n) is an additive white Gaussian noise part, and n is 0,1, …, L-1 and L is the number of discrete points in an observed time domain data range;
the multi-frequency harmonic signal is partially rewritten by the transformation relation between cosine and index to obtain an index signal model
And performing subsequent analysis and estimation by using the index model; wherein A isi、ωiAnd phiiThe amplitude, the digital frequency and the phase of the ith frequency component in the exponential signal model are shown, and P is the number of the multi-frequency exponential components;
step 2: performing N-point sliding window Fourier transform on a time domain sampling signal x (N) in an exponential form to obtain a sequence X (m);
and step 3: estimating the frequency of each component in the model by using a Prony algorithm expanded by a complex least square criterion for the transformed sequence X (m);
and 4, step 4: estimating the amplitude and the phase corresponding to each component of the exponential signal model by using an extended Prony algorithm on the basis of the frequency estimation value;
and 5: and recovering the harmonic component parameters of the initial power system signal according to the exponential model estimation value of the signal.
Further, the step 2 of obtaining the sequence x (m) by using a sliding window fourier transform comprises the following steps:
2.1) representing the frequency of each component as ωi＝2π/N(ki+δi) Wherein k isi∈{-N/2,…,0,…,N/2}，|δiThe | is less than or equal to 0.5 and respectively represents an integer part and a decimal part of each unknown frequency normalized by 2 pi/N, N is the window length of Fourier transform of a sliding window, and a spectral line value k corresponding to each secondary component is obtained through frequency rough estimationi；
2.2) calculating k separately for the signals to be analyzediThe result of the Fourier transform of the sliding window at the secondary spectral line is obtained as a sequence
2.3) adding the transformation results at the P frequency components to obtain the required sequence X (m).
Further, the step of performing frequency estimation by using the Prony algorithm after the complex least square criterion extension in the step 3 is as follows:
3.1) estimation models adopted according to the Prony algorithm
n-0, 1, …, L-1, sliding window fourier transform sequence still satisfies s (m) -a1S(m-1)+a2S(m-2)+…+aPS (m-P), thereby obtaining a characteristic equation a used for estimating parameters1z-1+a2z-2+…+aPz-PThe characteristic root of the equation is x (m) ═ a, where x (m) ═ a is the linear relationship between 1 and 0 and the frequency parameter to be estimated, taking noise interference into account1X(m-1)+a2X(m-2)+…+aPX (m-P); wherein, ci＝AiejφiIs a parameter containing time domain signal amplitude and phase information; z is a radical ofiIs the characteristic root of the equation; a isiI-1, 2, …, P representing the coefficients of the difference equation; s (m) is the spectral line k corresponding to each frequency component of the time-domain pure multi-frequency signal s (n)iProcessing a superimposed sequence of sliding window Fourier transform results; x (m) is the result of superposition of the SWDFT transform of the noisy multifrequency signal x (n) at the spectral lines corresponding to the P frequencies when the noise interference is considered;
3.2) defining vector Xm＝[X(m),X(m+1),…,X(m+L-N)]TThen Xm＝a1Xm-1+a2Xm-2+…+aPXm-PThe corresponding estimation error is
Searching for the overall mean square error according to the complex least squares criterion
Minimum value when estimating to obtain a characteristic polynomialCoefficient of (2)
i＝1,2,…,P；
3.3) substituting the estimated values of the equation coefficients into a characteristic polynomial for solving the signal, from the characteristic root z of the equationiObtaining frequency estimation values of the signal components
Further, in step 4, the amplitude and phase angle corresponding to each component are estimated, and the steps are as follows:
4.1) defining parameters containing amplitude and phase parameter information of each frequency component as
Substituting the estimation model of the Prony algorithm into the characteristic root and the sampled time domain signal to obtain an estimation value under the complex least square criterion
Wherein,
andamplitude information and phase parameters of each component in S (m) respectively, and corresponding to the signal component before transformation;
4.2) based on the estimated valueAnd calculating to obtain the amplitude and the phase of each frequency component in the exponential signal model by combining a mathematical expression that the sliding window Fourier exerts influence on the amplitude and the phase of each harmonic component.
Further, in step 5, for the estimated exponential model parameter, a component with a positive frequency estimation value is reserved, the corresponding amplitude estimation value is multiplied by 2, and the phase estimation value is kept unchanged, so that the original harmonic component parameter of the power system signal is obtained.
Has the advantages that: in the invention, a multi-frequency index model of the signal is derived from the parameter estimation problem of the power system signal containing harmonic pollution through the mathematical conversion relation between cosine and index, and the specific harmonic parameter estimation process of the scheme is carried out on the basis of the multi-frequency index model. Further, the signal is regarded as a superposition of several exponential components, i.e. the number of signal components is twice that of the actual signal, wherein half of the corresponding frequencies are negative values of the respective harmonic frequencies. The Sliding Window DFT (SWDFT) keeps the linear relation of the signal index model, inhibits the interference of noise, enhances the signal-to-noise ratio (SNR) of each component, and adopts a Prony model to carry out parameter estimation on the transformed sequence, so that the improvement of the SNR can bring about the improvement of the estimation precision. In addition, in order to enhance the robustness and estimation performance of the algorithm, the Prony algorithm adopted in the invention is established on the basis of the Complex Least Squares (CLS) criterion. Therefore, the harmonic parameter estimation method based on the SWDFT and the Prony algorithm can obtain higher estimation precision, and the robust anti-noise performance and the calculation complexity of the algorithm have obvious advantages.
Compared with the prior art, the invention has the following advantages: 1. the mathematical relationship between the SWDFT conversion result and the original time domain sequence and the signal-to-noise ratio gain brought by the mathematical relationship are fully utilized, and the anti-noise performance and the estimation precision of the parameter estimation method are improved. 2. Compared with the traditional Prony algorithm, the extended Prony algorithm adopting the CLS criterion reduces the equation solving error and enhances the robustness of the harmonic parameter estimation method.
Drawings
FIG. 1 is a diagram of mean square error of power grid harmonic signal frequency estimation as a function of signal-to-noise ratio;
fig. 2 is a graph of variation of amplitude and phase estimation errors of fundamental frequency components of a power grid signal with a signal-to-noise ratio, wherein fig. 2(a) is an amplitude estimation error curve, and fig. 2(b) is a phase estimation error curve;
fig. 3 is a graph of variation of amplitude and phase estimation error of 3-order harmonic component of a power grid signal with signal-to-noise ratio, where fig. 3(a) is an amplitude estimation error curve and fig. 3(b) is a phase estimation error curve;
fig. 4 is a graph of amplitude and phase estimation error of the power grid signal 5 th harmonic component as a function of signal-to-noise ratio, where fig. 4(a) is an amplitude estimation error curve and fig. 4(b) is a phase estimation error curve.
Detailed Description
The invention is further explained below with reference to the drawings.
A high-precision harmonic parameter estimation method based on sliding window Fourier transform (SWDFT) is disclosed, wherein a multifrequency signal of an electric power system under a noise background is as follows:
wherein s (t) represents a portion of the clean harmonic signal containing fundamental frequency components at time t, and q (t) represents a mean value of 0 and a variance of
M is the number of harmonics, fm、amAnd phimRespectively representing the frequency, amplitude and phase corresponding to the mth harmonic.
With a period TsDiscrete sampling is performed on the signal to obtain an L-point discrete sequence x (n):
in the formula of omegam＝2πfmTsS (n) is the pure multifrequency harmonic signal in the model at time n, and q (n) is the additive white gaussian noise part in the model. Since each sinusoidal component can be written as the addition of two complex exponential signals with frequencies opposite to each other, the signal model of the formula (2) can be rewritten as
Wherein A isi、ωiAnd phiiIs the amplitude, digital frequency and phase of the ith frequency component in an exponential signal model, and P is a multi-frequency finger in the model shown in formula (3)Number of component parts. Comparing with equation (2), it can be seen that the signal component in the exponential model becomes twice the original, i.e. P is 2M, the same negative frequency component is generated, the corresponding amplitude becomes 1/2, and the phase of the signal remains unchanged. Therefore, by estimating the model shown in formula (3), then retaining the partial components with positive frequency estimation values, and expanding the corresponding amplitudes to be twice of the estimation values, the harmonic estimation parameters of the initial power system signal can be obtained.
Note that the frequency of any component in the signal shown in equation (3) can be expressed as ωi＝2π/N(ki+δi) Wherein k isi∈{-N/2,…,0,…,N/2}，|δiAnd | is less than or equal to 0.5 is an integer part and a decimal part normalized by 2 pi/N of each unknown frequency, and N is an optional parameter and corresponds to the length of a sliding window of SWDFT transformation in the subsequent step. The number P of the multi-frequency exponential components and the normalized integer value k of the multi-frequency signal under the noise backgroundiAnd i is 1, …, and P can be accurately obtained by an FFT spectral peak search, Singular Value Decomposition (SVD) fixed-order equal-frequency rough estimation scheme.
For a pure multifrequency harmonic signal portion s (N) in the model, the transformation result at k-th spectral lines of the N-point SWDFT is as follows:
wherein m represents the starting time of the sliding window, and the starting time of the sliding window is delayed by m sampling points from the signal starting time; s (n + m) is a multi-frequency signal component corresponding to the time of n + m, and in the formula
From equation (4), the result of the SWDFT transform still satisfies the exponential model of equation (3), and the SWDFT transform changes only the amplitude and phase of each component. According to the Prony algorithm, the mathematical model adopted is a group of P exponential functions with arbitrary amplitude, phase, frequency and attenuation factor, and the function form of discrete time is
We use
As an approximation of s (n), in whichCharacteristic root of equation
αiAnd TsRespectively attenuation factor and sampling period.
The key point of the Prony algorithm is that the derivation finds that the fitting of the formula (6) is a homogeneous solution of a constant coefficient linear difference equation, namely a linear prediction difference equation satisfying the recursion
s(n)＝a1s(n-1)+a2s(n-2)+…+aPs(n-P) (7)
aiI-1, 2, …, P represents the coefficients of a difference equation or linear prediction coefficients, the corresponding characteristic polynomial being
Thus, from equation (4), the transformed sequence Sk(m) Linear prediction Structure which still conforms to the Prony Algorithm, i.e. Sk(m)＝a1Sk(m-1)+…+aPSk(m-P), we can estimate the parameters of each component using the Prony algorithm. Specifically, the Prony algorithm finds the coefficients of a characteristic polynomial from a matrix relation through an observed time sequence, calculates a characteristic root according to the coefficients, obtains the frequency parameter of the signal, and then obtains the amplitude and phase estimation values of the signal on the basis of the solution of the characteristic polynomial.
In combination with the formula (5), each time after SWDFT conversionThe amplitude of each frequency component decays as its frequency increases with distance from the corresponding angular frequency of the k spectral lines by 2 π k/N. For the ith frequency component, the integer value after frequency normalization is kiWhere k is equal to kiWith minimum time-amplitude attenuation, i.e. pure multifrequency signal s (n) at kiSWDFT results at the secondary line
The amplitude attenuation of the ith frequency component in the signal can be minimized.
Therefore, considering noise interference, for an actual noisy multi-frequency grid signal x (n), a spectral line k ═ k corresponding to the ith frequency componentiWhere the result of the N-point SWDFT transform is
Wherein m represents the number of delay samples of the sliding window relative to the initial signal time, and x (n + m) is the power grid signal corresponding to the n + m time. According to equations (4) and (5), the signal-to-noise ratio of the ith frequency component will also change compared to before the transform: in the time domain signal before SWDFT, the signal-to-noise ratio of the component is
The power of the transformed noise part q (n) becomes
The signal-to-noise ratio after conversion is thus
As shown in the formula (9), k is a value of a signal containing noise interferenceiSWDFT results at the secondary line
Enhanced angular frequency 2 pi kiThe signal-to-noise ratio of frequency components near/N and the corresponding parameter estimation performance are improved. If the frequency parameters of each component are to be accurately estimated, we need to do soThe signal is calculated and estimated P times by SWDFT, which brings a huge calculation burden. In order to improve the calculation efficiency as much as possible while ensuring the estimation accuracy, the addition situation of the P-path spectral line transformation results is considered next, and obviously, the sequence of pure harmonic parts
The linear relationship is still satisfied:
S(m)＝a1S(m-1)+a2S(m-2)+…+aPS(m-P) (10)
s (m) is spectral line k corresponding to each frequency component of time-domain multi-frequency signal s (n)iAnd (4) processing the superposed sequence of the Fourier transform results of the sliding window.
In order to analyze the signal-to-noise ratio change of the sequence before and after SWDFT in the noise interference environment, the S (m) is expanded in detail:
wherein k isjA normalized integer value representing the frequency of the jth signal component in the signal.
Therefore, for the transformation sequence under consideration of noise interference
In the case of a non-woven fabric,
for the signal to be analyzed at kiAs a result of the sliding window fourier transform at the secondary spectral line, the signal-to-noise ratio of the ith frequency component becomes larger than the original time-domain signal
The relation SNR can be found by mathematical derivationi≤SNRi,DFT≤SNRi,SWDFTThe inequality shows that the performance of parameter estimation using the sequence x (m) will be slightly lower than using the sequence x (m)
But significantly better than using the original time-domain sampled signal. Therefore, the algorithm in the invention adopts the sequence X (m) of the superposition of the transformation results at the spectral lines corresponding to the signal components for parameter estimation.
In addition, in order to enhance the robustness of the algorithm and reduce the error when solving the matrix equation in the Prony algorithm, a Complex Least Square (CLS) structure is adopted to estimate corresponding parameters, and the specific process is as follows:
(10) formula is scalable in noisy situations
Definition vector Xm＝[X(m),X(m+1),…,X(m+L-N)]TThe above formula can be written as
Xm＝a1Xm-1+a2Xm-2+…+aPXm-P(15)
The corresponding estimation error is
Then its self-conjugate matrix is
The objective of solving the characteristic equation shown in equation (8) for the CLS structure is to minimize the estimated overall mean square error
Order to
Then
The coefficients of the characteristic polynomial are calculated by solving the formula (17), and then the coefficients are substituted into the characteristic equation to solve P characteristic roots, and the corresponding frequency estimation value is
∠ (DEG) represents angle calculation, note that SWDFT conversion sequence X (m) is changed from L to L-N +1 in length compared with time domain signal x (N) by introducing N-point long sliding window conversion, so that parameter containing amplitude and phase information is obtained by Prony algorithm
Satisfy the requirement of
Similarly, the CLS criterion is still employed in estimating the magnitude and phase angle. Let X be [ X (0), X (1), …, X (L-N)]T，And b ═ b1,b2,…,bP]TThen the estimation equation using the CLS structure is
ZiIs a characteristic root vector containing frequency information, and b is a parameter containing amplitude and phase information
The expanded vector.
Calculating an estimated value from equation (21)Combining the formula (11) to obtain the amplitude and phase corresponding to the signal index model
In order to recover the harmonic component parameters of the initial power system signal, only the signal component with the frequency of a positive value in the estimation result needs to be selected, and the corresponding harmonic amplitude is corrected into an amplitude estimation value
Twice as much, the phase remains unchanged.
The invention is further illustrated below with reference to simulation results of an embodiment.
The harmonic power grid signal model is taken as follows:
x(t)＝20cos(2πf0t)+4cos(3×2πf0t+π/2)+3cos(5×2πf0t+π/3)+q(t)
wherein the fundamental frequency is f0Q (t) denotes the interference noise of the signal, which contains 3 rd and 5 th harmonics in addition to the fundamental frequency component. Taking the sampling frequency fs800Hz, the number of sampling points L is 250, the window length of SWDFT transform is selected to be N round (fs/f) 16, and round (cndot) represents rounding.
In a first group of simulations, the harmonic frequency of a power grid signal is estimated by using the algorithm provided by the invention, in order to clearly provide the advantages of the algorithm, the algorithm is compared with a frequency estimation result estimated by a Prony algorithm (under a CLS (common line-of-service) principle) by directly using a time domain signal, and a Mean Square Error (MSE) is used as a performance measurement value, wherein the larger the mean square error is, the worse the estimation performance of the algorithm is. Fig. 1 shows the frequency estimation error versus signal to noise ratio for the two algorithms described above. As shown in the figure, the solid line and the dotted line marked by the same color and the graph represent the frequency estimation mean square error of the algorithm disclosed by the invention and the time domain Prony algorithm respectively, and simulation shows that the algorithm provided by the invention can accurately estimate the frequency parameter of the harmonic component, wherein the estimation result of the fundamental frequency is particularly accurate, and the estimation error of the scheme is obviously lower than that of the time domain Prony algorithm.
On the basis of the frequency parameters, the amplitude and phase of each frequency component are estimated by using the Prony algorithm, and the amplitude and phase estimation errors of fundamental frequency, 3-order harmonic and 5-order harmonic components are respectively shown in FIGS. 2-4. The mean square error of the estimation of the same signal component by two algorithms in each simulation image is compared, and the estimation error of the proposed scheme in the aspects of harmonic component amplitude and phase estimation is smaller than that of a time domain Prony algorithm. In conclusion, the harmonic estimation algorithm provided by the invention has excellent performance in parameter estimation, and can accurately provide estimation results of frequency, amplitude and phase angle of each harmonic component.
The technical scheme of the invention can be applied to harmonic analysis, electric energy metering and electric energy quality monitoring of the electric power system.
Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present invention and not for limiting, and although the present invention has been described in detail with reference to the simulation results of the preferred embodiments, it should be understood by those skilled in the art that modifications or equivalent substitutions can be made on the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, which should be covered by the claims of the present invention.
Claims (3)
1. A high-precision harmonic parameter estimation method based on sliding window DFT is characterized by comprising the following steps:
step 1: collecting multi-frequency signals in an electric power system, and obtaining discrete signals x (n) at time n after samplingWherein a ism，ωm，φmRespectively representing the amplitude, digital frequency and phase of the mth harmonic, M is the number of harmonic contained in the signal, s (n) is a multi-frequency harmonic signal part, q (n) is an additive white Gaussian noise part, and the observed time domain data rangeSelecting n as 0,1, … and L-1, wherein L is the number of discrete points;
the multi-frequency harmonic signal is partially rewritten by the transformation relation between cosine and index to obtain an index signal model
And performing subsequent analysis and estimation by using the index model; wherein A isi、ωiAnd phiiThe amplitude, the digital frequency and the phase of the ith frequency component in the exponential signal model are shown, and P is the number of the multi-frequency exponential components;
step 2: performing N-point sliding window Fourier transform on a time domain sampling signal x (N) in an exponential form to obtain a sequence X (m);
and step 3: estimating the frequency of each component in the model by using a Prony algorithm expanded by a complex least square criterion for the transformed sequence X (m);
and 4, step 4: estimating the amplitude and the phase corresponding to each component of the exponential signal model by using an extended Prony algorithm on the basis of the frequency estimation value;
and 5: recovering harmonic component parameters of the initial power system signal according to the exponential model estimation value of the signal;
the step of performing frequency estimation by using the Prony algorithm after the complex least square criterion expansion in the step 3 comprises the following steps:
3.1) estimation models adopted according to the Prony algorithm
The sliding window fourier transform sequence still satisfies s (m) ═ a1S(m-1)+a2S(m-2)+…+aPS (m-P), thereby obtaining a characteristic equation a used for estimating parameters1z-1+a2z-2+…+aPz-PThe characteristic root of the equation is x (m) ═ a, where x (m) ═ a is the linear relationship between 1 and 0 and the frequency parameter to be estimated, taking noise interference into account1X(m-1)+a2X(m-2)+…+aPX (m-P); wherein, ci＝AiejφiIs a parameter containing time domain signal amplitude and phase information; z is a radical ofiIs a characteristic of an equationGetting root; a isiI-1, 2, …, P representing the coefficients of the difference equation; s (m) is the spectral line k corresponding to each frequency component of the time-domain pure multi-frequency signal s (n)iProcessing a superimposed sequence of sliding window Fourier transform results; x (m) is the result of superposition of the SWDFT transform of the noisy multifrequency signal x (n) at the spectral lines corresponding to the P frequencies when the noise interference is considered;
3.2) defining vector Xm＝[X(m),X(m+1),…,X(m+L-N)]TThen Xm＝a1Xm-1+a2Xm-2+…+aPXm-PThe corresponding estimation error is
Searching for the overall mean square error according to the complex least squares criterion
Minimum value at which the coefficients of the characteristic polynomial are estimated
3.3) substituting the estimated values of the equation coefficients into a characteristic polynomial for solving the signal, from the characteristic root z of the equationiObtaining frequency estimation values of the signal components
In step 4, the amplitude and phase angle corresponding to each component are estimated, and the steps are as follows:
4.1) defining parameters containing amplitude and phase parameter information of each frequency component as
Substituting the estimation model of the Prony algorithm into the characteristic root and the sampled time domain signal to obtain an estimation value under the complex least square criterionWherein A isi"and phii"the amplitude information and phase parameter of each component in S (m), respectively, and corresponding to the signal component before transformation;
4.2) based on the estimated value
And calculating to obtain the amplitude and the phase of each frequency component in the exponential signal model by combining a mathematical expression that the sliding window Fourier exerts influence on the amplitude and the phase of each harmonic component.
2. The sliding window DFT-based high precision harmonic parameter estimation method as claimed in claim 1, wherein: the step 2 of obtaining the sequence x (m) by using sliding window fourier transform comprises the following steps:
2.1) representing the frequency of each component as ωi＝2π/N(ki+δi) Wherein k isi∈{-N/2,…,0,…,N/2}，|δiThe | is less than or equal to 0.5 and respectively represents an integer part and a decimal part of each unknown frequency normalized by 2 pi/N, N is the window length of Fourier transform of a sliding window, and a spectral line value k corresponding to each secondary component is obtained through frequency rough estimationi；
2.2) calculating k separately for the signals to be analyzediThe result of the Fourier transform of the sliding window at the secondary spectral line is obtained as a sequence
2.3) adding the transformation results at the P frequency components to obtain the required sequence X (m).
3. The sliding window DFT-based high precision harmonic parameter estimation method as claimed in claim 1, wherein: and 5, reserving components with positive frequency estimation values for the estimated exponential model parameters, multiplying the corresponding amplitude estimation values by 2, and keeping the phase estimation values unchanged to obtain the original harmonic component parameters of the power system signals.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711265070.XA CN108037361B (en) | 2017-12-05 | 2017-12-05 | High-precision harmonic parameter estimation method based on sliding window DFT |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711265070.XA CN108037361B (en) | 2017-12-05 | 2017-12-05 | High-precision harmonic parameter estimation method based on sliding window DFT |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108037361A CN108037361A (en) | 2018-05-15 |
| CN108037361B true CN108037361B (en) | 2020-02-07 |
Family
ID=62095278
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711265070.XA Active CN108037361B (en) | 2017-12-05 | 2017-12-05 | High-precision harmonic parameter estimation method based on sliding window DFT |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108037361B (en) |
Families Citing this family (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109358223A (en) * | 2018-09-29 | 2019-02-19 | 海特尔机电工程技术（马鞍山）有限公司 | A kind of sliding window DFT harmonic current detecting method and device, storage medium |
| CN109490625B (en) * | 2018-11-14 | 2020-11-20 | 国网江苏省电力有限公司徐州供电分公司 | Harmonic signal analysis method based on sliding window and semi-definite programming |
| CN109683016A (en) * | 2019-01-18 | 2019-04-26 | 江苏理工学院 | A kind of harmonic amplitude analysis method |
| CN110008434B (en) * | 2019-03-20 | 2020-11-17 | 华中科技大学 | High-precision simple harmonic signal parameter estimation method |
| CN110806534B (en) * | 2019-12-04 | 2021-08-17 | 绵阳市维博电子有限责任公司 | Railway 25Hz phase-sensitive track signal detection method and system based on channel multiplexing |
| CN113032716B (en) * | 2019-12-24 | 2024-06-18 | 南京理工大学 | Harmonic and inter-harmonic analysis method based on windowed interpolation and Prony algorithm |
| CN112269054A (en) * | 2020-09-16 | 2021-01-26 | 国网安徽省电力有限公司六安供电公司 | Power adaptive algorithm based on improved Prony |
| CN112362968B (en) * | 2020-11-18 | 2021-07-02 | 华中科技大学 | Single-phase harmonic real-time extraction method based on pre-modulation CDSC and SDFT |
| CN112557751B (en) * | 2020-12-03 | 2023-07-18 | 东南大学 | Harmonic parameter estimation method based on DFT iteration method |
| CN112883318A (en) * | 2020-12-16 | 2021-06-01 | 中国空气动力研究与发展中心设备设计及测试技术研究所 | Multi-frequency attenuation signal parameter estimation algorithm of subtraction strategy |
| CN113449255B (en) * | 2021-06-15 | 2022-11-11 | 电子科技大学 | Improved method and device for estimating phase angle of environmental component under sparse constraint and storage medium |
| CN113552415A (en) * | 2021-07-28 | 2021-10-26 | 南京灿能电力自动化股份有限公司 | Ultrahigh harmonic measurement device and measurement method |
| CN114624513B (en) * | 2022-01-27 | 2024-06-25 | 清华大学 | Method and device for detecting phase of periodic signal with anti-harmonic interference |
| CN114859120B (en) * | 2022-06-17 | 2023-09-15 | 烟台东方威思顿电气有限公司 | Harmonic Analysis Method |
| CN115950529B (en) * | 2023-03-10 | 2023-06-09 | 天津大学 | Micro-angle resonance signal estimation method and device based on spectrum enhancement and electronic equipment |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101566649A (en) * | 2009-05-27 | 2009-10-28 | 重庆大学 | Harmonic detection method in a power system |
| CN103630742A (en) * | 2013-12-16 | 2014-03-12 | 国家电网公司 | Dynamic signal parameter acquisition method |
| CN107085140A (en) * | 2017-04-25 | 2017-08-22 | 东南大学 | Nonequilibrium system frequency estimating methods based on improved SmartDFT algorithms |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8116710B2 (en) * | 2009-06-04 | 2012-02-14 | Telefonaktiebolaget L M Ericsson (Publ) | Continuous sequential scatterer estimation |
-
2017
- 2017-12-05 CN CN201711265070.XA patent/CN108037361B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101566649A (en) * | 2009-05-27 | 2009-10-28 | 重庆大学 | Harmonic detection method in a power system |
| CN103630742A (en) * | 2013-12-16 | 2014-03-12 | 国家电网公司 | Dynamic signal parameter acquisition method |
| CN107085140A (en) * | 2017-04-25 | 2017-08-22 | 东南大学 | Nonequilibrium system frequency estimating methods based on improved SmartDFT algorithms |
Non-Patent Citations (2)
| Title |
|---|
| 基于总体最小二乘改进的SDFT三相交流电频率估计算法;刘勇 等;《东南大学学报(自然科学版)》;20171120;第47卷(第6期);1129-1134 * |
| 采用Prony算法的谐波和间谐波分析;常晓颖，吴茜琼;《低压电器》;20110330(第2011年06期);38-41 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108037361A (en) | 2018-05-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108037361B (en) | High-precision harmonic parameter estimation method based on sliding window DFT | |
| CN110333389B (en) | Sinusoidal signal frequency estimation method based on interpolation DFT | |
| CN107102255B (en) | Single ADC acquisition channel dynamic characteristic test method | |
| CN107085140B (en) | Nonequilibrium system frequency estimating methods based on improved SmartDFT algorithm | |
| CN104142425B (en) | Phase matching method for sinusoidal signal frequency estimation | |
| CN109856455B (en) | Real-time repeated conversion type attenuation signal parameter estimation method | |
| CN111222088B (en) | Improved method for estimating weighted power harmonic amplitude of flat-top self-convolution window | |
| Wang et al. | Estimation of damping factor and signal frequency for damped sinusoidal signal by three points interpolated DFT | |
| CN114781196A (en) | Harmonic detection method based on sparse acquisition model | |
| CN112881796A (en) | Multi-frequency real signal frequency estimation algorithm for spectrum leakage correction | |
| CN114842867A (en) | DFT-based audio sinusoidal signal frequency estimation method and system | |
| CN114280366B (en) | Sinusoidal signal frequency estimation method based on improved frequency interpolation algorithm | |
| JP2014153354A (en) | Method for estimating frequencies and phases in three phase power system | |
| CN112557751B (en) | Harmonic parameter estimation method based on DFT iteration method | |
| CN112883318A (en) | Multi-frequency attenuation signal parameter estimation algorithm of subtraction strategy | |
| Djurovic | Estimation of the sinusoidal signal frequency based on the marginal median DFT | |
| CN112816779A (en) | Harmonic real signal parameter estimation method for analytic signal generation | |
| CN108801296B (en) | Sensor frequency response function calculation method based on error model iterative compensation | |
| CN113468474B (en) | Power grid frequency estimation method based on root Mini-Norm | |
| Ruan et al. | Improved Prony method for high-frequency-resolution harmonic and interharmonic analysis | |
| CN115032453A (en) | Multi-frequency dynamic phasor measurement method | |
| CN115166354B (en) | Noisy real sinusoidal signal frequency estimation method based on zero-padding difference DFT | |
| Xiao et al. | A neural network model for power system inter-harmonics estimation | |
| CN113392813B (en) | Method and system for accurately identifying main frequency of vibration signal | |
| CN108173259B (en) | Sine frequency estimation method based on unit constraint minimum mean square error |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |