CN103454494B - A kind of high-precision harmonic analysis method - Google Patents

Abstract

The invention discloses a kind of high-precision harmonic analysis method improving realization on standard synchronous DFT basis.The method comprises: the position that the accurate synchronous DFT of application carries out the sampling of frequency analysis time-frequency domain changes according to the drift of signal frequency, and namely described frequency domain sampling position is μ 2 π/N, wherein: μ is the drift of signal frequency, without μ during drift, is 1.The present invention includes a thought: variable fence, i.e. the position of frequency analysis time-frequency domain sampling changes according to the drift of signal frequency; Technology of the present invention, the field contributing to the application frequency analyses such as electric energy quality monitoring, electronic product production testing, electric appliances monitoring obtains the information such as the amplitude of each harmonic, initial phase angle and frequency more accurately.

Description

A kind of high-precision harmonic analysis method

The application is divisional application, the application number of original application: 201110245638.8, invention and created name: a kind of harmonic analysis method, the applying date: 2011-8-24.

Technical field

The present invention relates to a kind of high-precision harmonic analysis method.

Background technology

Frequency analysis technology is widely used in various fields such as electric energy quality monitoring, electronic product production testing, electric appliances monitoring, is the important technical of carrying out power system monitor, quality inspection, monitoring of tools.The most widely used technology of current frequency analysis is discrete Fourier transformation (DFT) and Fast Fourier Transform (FFT) (FFT).The frequency analysis technology that quasi-synchronous sampling technique and DFT technology combine can improve the precision of frequency analysis, and its formula is:

$\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\cos \left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\sin \left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,$

In formula: k is the number of times (as first-harmonic k=1,3 subharmonic k=3) needing the harmonic wave obtained; Sin and cos is respectively sine and cosine functions; And a _kand b _kbe respectively real part and the imaginary part of k subharmonic; N is iterations; W is determined by integration method, when adopting muiltiple-trapezoid integration method, and W=nN; γ _iit is a weighting coefficient; for all weighting coefficient sums; I-th sampled value that f (i) is analysis waveform; N is sampling number in the cycle.

In engineer applied, frequency analysis is always carried out the sampling of finite point and is difficult to accomplish the synchronized sampling of stricti jurise.Like this, when applying accurate synchronous DFT and carrying out frequency analysis, the long scope caused due to truncation effect will be there is and to leak and the short scope that causes due to fence effect is leaked, make analysis result precision not high, even not credible.

Summary of the invention

The technical problem to be solved in the present invention is to provide a kind of high-precision harmonic analysis method, effectively to improve the analytical error of accurate synchronous DFT frequency analysis technology, obtain high-precision frequency analysis result, thus improve the quality of field instrument and equipment and the validity of condition adjudgement such as electric energy quality monitoring, electronic product production testing, electric appliances monitoring based on frequency analysis theory.

For solving the problems of the technologies described above, high-precision harmonic analysis method provided by the invention comprises: the position that the accurate synchronous DFT of application carries out the sampling of frequency analysis time-frequency domain changes according to the drift of signal frequency, namely described frequency domain sampling position is μ 2 π/N, wherein: μ is the drift of signal frequency, during nothing drift, μ is 1.

Harmonic analysis method of the present invention, based on the thought of variable fence, is realized by 5 analytical procedures.

The thought of variable fence: the main cause of accurate synchronous DFT analytical error is the drift of the signal frequency position that causes spectrum peak to occur and ideal position generation deviation, if the analysis result still obtained to carry out sampling in a frequency domain according to 2 π/N is extremely incorrect.Variable fence refers to: the position of frequency domain sampling be not 2 fixing π/N, but change according to the drift of signal frequency, namely frequency sampling position is μ 2 π/N (μ is the drift of signal frequency).Frequency domain sampling fence changes the position that accurately can estimate each harmonic peak value and occur along with the drift of signal frequency, and then obtains high-precision amplitude and phase angle information.

Frequency analysis step of the present invention is as follows:

(1) equal interval sampling W+2 sampling number is according to { f (i), i=0,1, (W is determined by selected integration method W+1}, the present invention does not specify a certain integration method, and conventional integration method has muiltiple-trapezoid integration method W=nN, complexification rectangular integration method W=n (N-1), iterative Simpson integration method W=n (N-1)/2 etc., can select suitable integration method according to the actual conditions of the present invention's application.General more satisfactory with muiltiple-trapezoid integration method effect.); (2) the accurate synchronous DFT formula of application from sampled point i=0 $\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\cos \left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\sin \left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,$ Analyze W+1 data acquisition fundamental information with

(3) accurate synchronous DFT formula is applied from sampled point i=1 $\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f(i+1)\cos \left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f(i+1)\sin \left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,$ Analyze W+1 data acquisition fundamental information with

(4) application of formula $\mathrm{\μ}=N\frac{{\mathrm{tg}}^{-1}\[\frac{F_{a0}^n\left(1\right)}{F_{b0}^n\left(1\right)}\]-{\mathrm{tg}}^{-1}\[\frac{F_{a0}^n\left(0\right)}{F_{b0}^n\left(0\right)}\]}{2\mathrm{\π}}$ Calculate the frequency drift μ of signal;

(5) application of formula $\left\{\begin{array}{c}{a}_k=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\cos \left(k\frac{\mathrm{\μ}2\mathrm{\π}}{N}i\right)\\ b_k=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\sin \left(k\frac{\mathrm{\μ}2\mathrm{\π}}{N}i\right)\end{array}\right.$ Calculate amplitude and the phase angle of each harmonic.

Accurate synchronous DFT frequency analysis can suppress long scope to be leaked effectively, the main cause of its spectrum leakage is that the short scope that signal frequency drift causes is leaked, and signal frequency drift causes principal character that short scope leaks is spectrum peak-to-peak value occurs position along with signal frequency drift synchronous change, so variable fence frequency domain sample effectively can catch the position of spectrum peak-to-peak value appearance according to signal drift, thus obtain high-precision harmonic information.

Equal interval sampling is cycle T and the frequency f (if power frequency component frequency f is 50Hz, the cycle is 20mS) that basis carries out the ideal signal of frequency analysis, and N point of sampling in one-period, namely sample frequency is f _s=Nf, and N>=64.

Sampling W+2 described sampling number certificate does corresponding selection according to selected integration method, according to muiltiple-trapezoid integration method, then W=nN; According to complexification rectangular integration method, then W=n (N-1); According to iterative Simpson integration method, then W=n (N-1)/2.Then according to sample frequency f _s=Nf, acquisition sampled point data sequence f (i), i=0,1 ..., W+1}, n>=3, finally carry out frequency analysis to this data sequence.

An iteration coefficient γ _idetermined by integration method, ideal period sampled point N and iterations n, concrete derivation see document [Dai Xianzhong. quasi-synchro sampling application in some problem [J]. electrical measurement and instrument, 1988, (2): 2-7.].

for all weighting coefficient sums.

A _kand b _kfor imaginary part and the real part of k subharmonic, according to a _kand b _kjust can obtain harmonic amplitude and initial phase angle.

The drift μ of signal frequency obtains according to the fixed relationship of sampling number N in neighbouring sample point first-harmonic phase angle difference and ideal period, and the drift μ of signal frequency also can be used for the frequency f revising first-harmonic ₁with the frequency f of higher hamonic wave _k.

Adopt above-mentioned high precision frequency analysis technology, also namely based on the frequency analysis technology of variable fence thought, there is following technical advantage:

(1) high-precision frequency analysis result.No matter the analysis result that frequency analysis technology of the present invention obtains is that amplitude or phase angle error improve more than 4 orders of magnitude.

(2) frequency analysis technology of the present invention fundamentally solves the low problem of accurate synchronous DFT analysis precision, and without the need to carrying out complicated inverting and correction, algorithm is simple.

(3) relative to the synchronous DFT of standard, frequency analysis technology of the present invention only needs increase sampled point just to solve the large problem of accurate synchronous DFT analytical error, is easy to realize.

(4) applying the present invention and improve existing instrument and equipment, is technically feasible, and does not need any hardware spending of increase that analysis result just can be made can to improve more than 4 orders of magnitude.

(5) variable fence thought is also applicable to the frequency analysis process of carrying out successive ignition and non-once iteration too, now only needs a Breaking Recurrently to become successive ignition realization just passable.One time iteration is the same with successive ignition in essence, and just when calculating, successive ignition carries out decoupled method, and an iteration is that the process of successive ignition is merged into iteration coefficient γ _iin once calculated, so the present invention is equally applicable to successive ignition process.

Embodiment

A kind of high precision frequency analysis technology of the present invention, comprises the following steps:

First, an equal interval sampling W+2 sampled point, with obtain analyzed signal discrete series f (k), k=0,1 ..., W+1}.W is determined jointly by sampling number N in integration method, iterations n and ideal period.Equal interval sampling refers to determines sample frequency f according to the frequency f (if power frequency component frequency is 50Hz, the cycle is 20mS) of the ideal signal carrying out frequency analysis _s=Nf, at sample frequency f _seffect under to sample equably in one-period N point.Usually, periodic sampling point N=64 or more just can obtain good frequency analysis result, and iterations n=3-5 just can obtain comparatively ideal frequency analysis result.Integration method has muiltiple-trapezoid integration method W=nN, complexification rectangular integration method W=n (N-1), Simpson integration method W=n (N-1)/2 etc. multiple, can select according to actual conditions.

Secondly, the accurate synchronous DFT formula of application from sampled point k=0

$\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\cos \left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\sin \left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,$ Analyze W+1 data acquisition fundamental information with wherein, an iteration coefficient γ _idetermined by integration method, ideal period sampled point N and iterations n, and for all weighting coefficient sums.

Again, accurate synchronous DFT formula is applied from sampled point k=1 $\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f(i+1)\cos \left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f(i+1)\sin \left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,$ Analyze W+1 data acquisition fundamental information with

Then, application of formula $\mathrm{\μ}=N\frac{{\mathrm{tg}}^{-1}\[\frac{F_{a0}^n\left(1\right)}{F_{b0}^n\left(1\right)}\]-{\mathrm{tg}}^{-1}\[\frac{F_{a0}^n\left(0\right)}{F_{b0}^n\left(0\right)}\]}{2\mathrm{\π}}$ Calculate the frequency drift μ of signal.After obtaining frequency drift μ, can according to sample frequency f _scalculate with sampling number N in ideal period and obtain the first-harmonic of analyzed signal and the frequency f of higher hamonic wave.

Finally, apply $\left\{\begin{array}{c}{a}_k=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\cos \left(k\frac{\mathrm{\μ}2\mathrm{\π}}{N}i\right)\\ b_k=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)\sin \left(k\frac{\mathrm{\μ}2\mathrm{\π}}{N}i\right)\end{array}\right.$ Calculate the real part a of k subharmonic _kwith imaginary part information b _k, and then according to formula: calculate amplitude P _k, and according to formula: calculate initial phase angle

Those skilled in the art will be appreciated that, above embodiment is only used to the present invention is described, and not as limitation of the invention, the present invention can also be changing into more mode, as long as in spirit of the present invention, all will drop in Claims scope of the present invention the change of the above embodiment, modification.

Claims (1)

1. a harmonic analysis method, is characterized in that comprising the following steps:

(1), an equal interval sampling W+2 sampling number certificate: f (i), i=0,1 ..., W+1};

(2), the accurate synchronous DFT formula of application from sampled point i=0:

$\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)cos\left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)sin\left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,$ Analyze W+1 data acquisition fundamental information with

(3), accurate synchronous DFT formula is applied from sampled point i=1:

$\left\{\begin{array}{c}{a}_k=2F_{ak}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f(i+1)cos\left(k\frac{2\mathrm{\π}}{N}i\right)\\ b_k=2F_{bk}^n=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f(i+1)sin\left(k\frac{2\mathrm{\π}}{N}i\right)\end{array}\right.,,$

Analyze W+1 data acquisition fundamental information

(4), application of formula: $\mathrm{\μ}=N\frac{{\mathrm{tg}}^{-1}\[\frac{F_{a0}^n\left(1\right)}{F_{b0}^n\left(1\right)}\]-{\mathrm{tg}}^{-1}\[\frac{F_{a0}^n\left(0\right)}{F_{b0}^n\left(0\right)}\]}{2\mathrm{\π}},$

Calculate the frequency drift μ of signal;

(5) application of formula: $\left\{\begin{array}{c}{a}_k=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)cos\left(k\frac{\mathrm{\μ}2\mathrm{\π}}{N}i\right)\\ b_k=\frac{2}{Q}\underset{i=0}{\overset{W}{\Σ}}{\mathrm{\γ}}_{i}f\left(i\right)sin\left(k\frac{\mathrm{\μ}2\mathrm{\π}}{N}i\right)\end{array}\right.$

, calculate amplitude and the phase angle of each harmonic;

Described equal interval sampling is cycle T and the frequency f that basis carries out the ideal signal of frequency analysis, and N point of sampling in one-period, namely sample frequency is f _s=Nf, and N>=64; Frequency f is 50Hz;

Sampling W+2 described sampling number certificate does corresponding selection according to selected integration method, according to muiltiple-trapezoid integration method, then W=nN; According to complexification rectangular integration method, then W=n (N-1); According to iterative Simpson integration method, then W=n (N-1)/2;

Then according to sample frequency f _s=Nf, acquisition sampled point data sequence f (i), i=0,1 ..., W+1}, n>=3, n are iterationses, and i-th sampled value that f (i) is analysis waveform, finally carries out frequency analysis to this data sequence;

for all weighting coefficient sums; γ _iit is an iteration coefficient; a _kand b _kfor imaginary part and the real part of k subharmonic, according to a _kand b _kjust can obtain harmonic amplitude and initial phase angle;

The drift μ of signal frequency obtains according to the fixed relationship of sampling number N in neighbouring sample point first-harmonic phase angle difference and ideal period, and the drift μ of signal frequency also can be used for the frequency f revising first-harmonic ₁with the frequency f of higher hamonic wave _k.