CN104897961A - Three spectral line interpolation FFT harmonic wave analysis method and system based on multiplication window function - Google Patents

Abstract

The invention discloses a three spectral line interpolation FFT (Fast Fourier Transform) harmonic wave analysis method and system based on a multiplication window function, relating to the harmonic wave analysis field. The method comprises the steps of: constructing a multiplication window function; pretreating signals; determining three spectral lines; calculating a correction formula of a three spectral line interpolation algorithm; calculating a fundamental wave parameter; determining a harmonic wave parameter; and performing error analysis. The invention provides a new window function construction mode, utilizes a routine window function to perform different multiplication combinations, firstly realizes construction of a multiplication window function, and applies the multiplication window function to a three spectral line interpolation FFT algorithm to perform harmonic wave analysis and calculate a harmonic wave parameter. According to a plurality of common cosine window function calculation instances, the multiplication window function of the invention has higher accuracy compared with a routine window function interpolation algorithm, and can realize high precision harmonic wave measurement.

Description

Based on three spectral line interpolation FFT harmonic analysis method and systems of cross multiplication window function

Technical field

The present invention relates to frequency analysis field, specifically relate to a kind of three spectral line interpolation FFT harmonic analysis method and systems based on cross multiplication window function.

Background technology

Along with the use of a large amount of non-linear power electronic devices, make harmonic pollution problems in electrical network day by day serious.Harmonic problem not only worsens the quality of power supply, also affects greatly the safety and stability of electrical network and economical operation.Therefore, high-acruracy survey is carried out to harmonic parameters in system, for minimizing Harmfulness Caused by Harmonics, safeguard electricity net safety stable, Effec-tive Function is very necessary.

The key issue of harmonic detecting is: under non-synchronous sampling, how to solve spectral leakage and fence effect.At present, common detection harmonic wave and the method for m-Acetyl chlorophosphonazo mainly contain: Fourier transform, machine learning method, wavelet transformation, power Spectral Estimation, Hilbert-Huang transform.

Analyze in harmonic wave in Fourier transform, method therefor is: use various special window function, block signal, then carries out frequency analysis in conjunction with spectral line interpolation FFT (Fast Fourier Transform, Fast Fourier Transform (FFT)).

Conventional window function has: the Chinese peaceful (Hanning) window function, Blacknam (Blackman) window function, this (Blackman-Harris) window function of the Blacknam Chinese, Nuttall (Nuttall) window function, Lay husband Vincent (Rife-Vincent) window function and various combination window function.

On windowing basis, the people such as D.Agrez and Pang Hao propose the correction algorithm of two spectral line separately, and the people such as Wu Jing, Niu Shengsuo and Huang Dongmei propose three spectral line correction algorithms.These improve the impact reducing spectrum leakage and fence effect, improve the accuracy of frequency analysis.

But in the actual use of engineering, conventional window function interpolation algorithm still cannot meet high-precision frequency analysis requirement.

Summary of the invention

The object of the invention is the deficiency in order to overcome above-mentioned background technology, a kind of three spectral line interpolation FFT harmonic analysis method and systems based on cross multiplication window function being provided, with custom window function ratio, there is higher precision.

The invention provides a kind of three spectral line interpolation FFT harmonic analysis methods based on cross multiplication window function, comprise the following steps:

S1, structure cross multiplication window function:

The general formula of multiplication window function is produced by multiple window function product, and the general formula of multiplication window function w (n) is:

Wherein, n is the ordinal number of sampled point, and n is natural number; M is the number of the basic window function participating in structure multiplication window, and m is positive integer; w _in () is i-th basic window function expression formula, i=1,2 ... m; c _ifor participating in the number of i-th basic window function, be called w _ithe sub-order of (n); c is total order of cross multiplication window function; According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), namely construct multi-form multiplication window function;

Work as w _iwhen () expression formula is identical n, when namely adopted basic window is identical, w (n) is multiplication window function, works as w _iwhen () expression formula is different n, when namely adopted basic window is different, w (n) is cross multiplication window function;

The expression formula of three conventional window functions is respectively w ₁(n), w ₂(n), w ₃(n), total order c=c ₁+ c ₂+ c ₃, to c maximal value, order: c=3;

According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), construct following 9 kinds of multiplication window functions:

Work as c ₁=3, c ₂=0, c ₃when=0, construct and take advantage of 300 window functions mutually;

Work as c ₁=2, c ₂=1, c ₃when=0, construct and take advantage of 210 window functions mutually;

Work as c ₁=2, c ₂=0, c ₃when=1, construct and take advantage of 201 window functions mutually;

Work as c ₁=1, c ₂=2, c ₃when=0, construct and take advantage of 120 window functions mutually;

Work as c ₁=1, c ₂=0, c ₃when=2, construct and take advantage of 102 window functions mutually;

Work as c ₁=0, c ₂=3, c ₃when=0, construct and take advantage of 030 window function mutually;

Work as c ₁=1, c ₂=1, c ₃when=1, construct and take advantage of 111 window functions mutually;

Work as c ₁=0, c ₂=2, c ₃when=1, construct and take advantage of 021 window function mutually;

Work as c ₁=0, c ₂=0, c ₃when=3, construct and take advantage of 003 window function mutually;

Wherein, take advantage of 300 window functions mutually, take advantage of 030 window function mutually, take advantage of 003 window function mutually, these three kinds of window functions belong to multiplication window function, all the other the 6 kinds cross multiplication window functions for constructing;

S2, Signal Pretreatment:

Mutual inductor gathers power network signal, and power network signal x (n) collected by mutual inductor, is transferred to host computer; Add cross multiplication window function w (n) to power network signal x (n) to block, obtain windowing signal x _w(n):

x _w(n)＝x(n)w(n) (2)

After FFT conversion is carried out to the windowing signal of formula (2), obtain windowing FFT frequency spectrum:

Wherein, the frequency spectrum that W () is window function, k is positive integer, and X (k) represents the frequency spectrum of kth subharmonic, and Ak is the amplitude of kth subharmonic, and j represents imaginary unit, and e is the truth of a matter of natural logarithm, for the initial phase of kth subharmonic, harmonic wave is first-harmonic for the first time, f _sfor sample frequency, f ₀for fundamental frequency, Δ f is discrete frequency intervals, and Δ f=f _s/ N;

Order: k ₀=f ₀/ Δ f, k ₀for the position of spectral line of real frequency spectrum,

Ignore the impact of negative frequency point place secondary lobe, formula (3) becomes:

S3, determine three spectral lines:

At the windowing FFT spectrum peak near zone that S2 obtains, k ₀three spectral lines that place's Frequency point is larger are respectively: kth ₁, k ₂, k ₃root, k ₁, k ₂, k ₃be positive integer, k ₁~ k ₃pass be: k2=k ₁+ 1, k ₃=k ₂+ 1, the amplitude that these three spectral lines are corresponding is respectively y ₁, y ₂, y ₃;

Note variable α=k-k ₂, then-0.5≤α≤0.5;

Another note variable $\mathrm{\β}=\frac{y_3-y_1}{y_2}---\left(5\right)$

S4, calculate the correction formula of three spectral line interpolation algorithms:

Obtain according to formula (4) and (5):

$\mathrm{\β}=\frac{\left|W(-\mathrm{\α}+1)\right|-\left|W(-\mathrm{\α}-1)\right|}{\left|W(-\mathrm{\α})\right|}---\left(6\right)$

Polynomial Approximation Method is adopted to calculate odd function β=g ^-1(α), expression formula is:

α≈p ₁₁×β+p ₁₃×β ³+…+p _1pβ ^p(7)

P ₁₁, p ₁₃; P _1pfor the odd item coefficient of approximation by polynomi-als, p is odd number;

According to formula (4), try to achieve the amplitude A of power network signal i-th subharmonic _i:

$A_i=2y_i\left|W(k-\frac{f_0}{\mathrm{\Δf}})\right|---\left(8\right)$

Wherein, i is positive integer, y _ifor the amplitude of the i-th subharmonic after windowing FFT;

Consider y ₂be from the nearest spectral line of true spectral line point, obtain:

$A_1=\frac{2(y_3+2y_2+y_1)}{\left|W(-\mathrm{\α}+1)\right|+2\left|W(-\mathrm{\α})\right|+\left|W(-\mathrm{\α}-1)\right|}---\left(9\right)$

During N > 1000, window function coefficient is real coefficient, and formula (9) is expressed as:

A ₁＝N ^-1(y ₃+2y ₂+y ₁)u(α)

Wherein, u (α) is correction formula, and is even function, and approximating polynomial is not containing odd item;

Three spectral line correction approximating polynomials are as follows:

u(α)＝(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (10)

In formula (10), p ₂₀, p ₂₂p _2dfor the even term coefficient of approximation by polynomi-als, d is the highest order of matching, and d is even number;

S5, calculating first-harmonic parameter:

Calculate fundamental frequency f ₀, fundamental voltage amplitude A ₁:

f ₀＝k·Δf＝(k ₂+α)Δf (11)

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (12)

According to formula (4), calculate the phase place of first-harmonic:

Copy asking for of first-harmonic parameter, carry out the analysis of each harmonic parameter according to formula (6), (7), (9), (11), (12), (13);

S6, determine harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, step S3 ~ S5 is repeated, until all harmonic parameters calculate complete;

S7, carry out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compare with custom window function interpolation arithmetic accuracy.

On the basis of technique scheme, described power network signal comprises current signal, voltage signal.

On the basis of technique scheme, the position of spectral line k of described real frequency spectrum ₀for decimal.

The present invention also provides a kind of three spectral line interpolation FFT frequency analysis systems based on cross multiplication window function, comprise cross multiplication window function tectonic element, Signal Pretreatment unit, spectral line determining unit, correction formula computing unit, first-harmonic parameter calculation unit, harmonic parameters determining unit, error analysis unit, wherein:

Described cross multiplication construction of function unit, for constructing cross multiplication window function:

The general formula of multiplication window function is produced by multiple window function product, and the general formula of multiplication window function w (n) is:

Wherein, n is the ordinal number of sampled point, and n is natural number; M is the number of the basic window function participating in structure multiplication window, and m is positive integer; w _in () is i-th basic window function expression formula, i=1,2 ... m; c _ifor participating in the number of i-th basic window function, be called w _ithe sub-order of (n); c is total order of cross multiplication window function; According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), namely construct multi-form multiplication window function;

Work as w _iwhen () expression formula is identical n, when namely adopted basic window is identical, w (n) is multiplication window function, works as w _iwhen () expression formula is different n, when namely adopted basic window is different, w (n) is cross multiplication window function;

The expression formula of three conventional window functions is respectively w ₁(n), w ₂(n), w ₃(n), total order c=c ₁+ c ₂+ c ₃, to c maximal value, order: c=3;

According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), construct following 9 kinds of multiplication window functions:

Work as c ₁=3, c ₂=0, c ₃when=0, construct and take advantage of 300 window functions mutually;

Work as c ₁=2, c ₂=1, c ₃when=0, construct and take advantage of 210 window functions mutually;

Work as c ₁=2, c ₂=0, c ₃when=1, construct and take advantage of 201 window functions mutually;

Work as c ₁=1, c ₂=2, c ₃when=0, construct and take advantage of 120 window functions mutually;

Work as c ₁=1, c ₂=0, c ₃when=2, construct and take advantage of 102 window functions mutually;

Work as c ₁=0, c ₂=3, c ₃when=0, construct and take advantage of 030 window function mutually;

Work as c ₁=1, c ₂=1, c ₃when=1, construct and take advantage of 111 window functions mutually;

Work as c ₁=0, c ₂=2, c ₃when=1, construct and take advantage of 021 window function mutually;

Work as c ₁=0, c ₂=0, c ₃when=3, construct and take advantage of 003 window function mutually;

Wherein, take advantage of 300 window functions mutually, take advantage of 030 window function mutually, take advantage of 003 window function mutually, these three kinds of window functions belong to multiplication window function, all the other the 6 kinds cross multiplication window functions for constructing;

Described Signal Pretreatment unit, for carrying out pre-service to signal:

Mutual inductor gathers power network signal, and power network signal x (n) collected by mutual inductor, is transferred to host computer; Add cross multiplication window function w (n) to power network signal x (n) to block, obtain windowing signal x _w(n):

x _w(n)＝x(n)w(n) (2)

After FFT conversion is carried out to the windowing signal of formula (2), obtain windowing FFT frequency spectrum:

Wherein, the frequency spectrum that W () is window function, k is positive integer, and X (k) represents the frequency spectrum of kth subharmonic, A _kfor the amplitude of kth subharmonic, j represents imaginary unit, and e is the truth of a matter of natural logarithm, for the initial phase of kth subharmonic, harmonic wave is first-harmonic for the first time, f _sfor sample frequency, f0 is fundamental frequency, and Δ f is discrete frequency intervals, and Δ f=f _s/ N;

Order: k ₀=f ₀/ Δ f, k ₀for the position of spectral line of real frequency spectrum,

Ignore the impact of negative frequency point place secondary lobe, formula (3) becomes:

Described spectral line determining unit, for determining three spectral lines:

At the windowing FFT spectrum peak near zone that S2 obtains, k ₀three spectral lines that place's Frequency point is larger are respectively: kth ₁, k ₂, k ₃root, k ₁, k ₂, k ₃be positive integer, k ₁~ k ₃pass be: k ₂=k ₁+ 1, k ₃=k ₂+ 1, the amplitude that these three spectral lines are corresponding is respectively y ₁, y ₂, y ₃;

Note variable α=k-k ₂, then-0.5≤α≤0.5;

Another note variable $\mathrm{\β}=\frac{y_3-y_1}{y_2}---\left(5\right);$

Described correction formula computing unit, for calculating the correction formula of three spectral line interpolation algorithms:

Obtain according to formula (4) and (5):

$\mathrm{\β}=\frac{\left|W(-\mathrm{\α}+1)\right|-\left|W(-\mathrm{\α}-1)\right|}{\left|W(-\mathrm{\α})\right|}---\left(6\right)$

Polynomial Approximation Method is adopted to calculate odd function β=g ^-1(α), expression formula is:

α≈p ₁₁×β+p ₁₃×β ³+…+p _1pβ ^p(7)

P ₁₁, p ₁₃; P _1pfor the odd item coefficient of approximation by polynomi-als, p is odd number;

According to formula (4), try to achieve the amplitude A of power network signal i-th subharmonic _i:

$A_i=2y_i\left|W(k-\frac{f_0}{\mathrm{\Δf}})\right|---\left(8\right)$

Wherein, i is positive integer, y _ifor the amplitude of the i-th subharmonic after windowing FFT;

Consider y ₂be from the nearest spectral line of true spectral line point, obtain:

$A_1=\frac{2(y_3+2y_2+y_1)}{\left|W(-\mathrm{\α}+1)\right|+2\left|W(-\mathrm{\α})\right|+\left|W(-\mathrm{\α}-1)\right|}---\left(9\right)$

During N > 1000, window function coefficient is real coefficient, and formula (9) is expressed as:

A ₁＝N ^-1(y ₃+2y ₂+y ₁)u(α)

Wherein, u (α) is correction formula, and is even function, and approximating polynomial is not containing odd item;

Three spectral line correction approximating polynomials are as follows:

u(α)＝(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (10)

In formula (10), p ₂₀, p ₂₂p _2dfor the even term coefficient of approximation by polynomi-als, d is the highest order of matching, and d is even number;

Described first-harmonic parameter calculation unit, for calculating first-harmonic parameter:

Calculate fundamental frequency f ₀, fundamental voltage amplitude A ₁:

f ₀＝k·Δf＝(k ₂+α)Δf (11)

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (12)

According to formula (4), calculate the phase place of first-harmonic:

Copy asking for of first-harmonic parameter, carry out the analysis of each harmonic parameter according to formula (6), (7), (9), (11), (12), (13);

Described harmonic parameters determining unit, for determining harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, described spectral line determining unit, correction formula computing unit, first-harmonic parameter calculation unit repeat to calculate, until all harmonic parameters calculate complete;

Described error analysis unit, for carrying out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compares with custom window function interpolation arithmetic accuracy.

On the basis of technique scheme, described power network signal comprises current signal, voltage signal.

On the basis of technique scheme, the position of spectral line k0 of described real frequency spectrum is decimal.

Compared with prior art, advantage of the present invention is as follows:

The present invention proposes a kind of make of new window function, utilize custom window function to carry out different multiplicative combination, realize structure cross multiplication window function first, cross multiplication window function is applied in three spectral line interpolation FFT algorithms, carry out frequency analysis, calculate harmonic parameters.Multiple conventional cosine window function calculated examples shows, the cross multiplication window function that the present invention constructs, relative to custom window function interpolation algorithm, has higher accuracy, achieves the high-acruracy survey of harmonic wave, has higher practical value.

Accompanying drawing explanation

Fig. 1 is the process flow diagram based on three spectral line interpolation FFT harmonic analysis methods of cross multiplication window function in the embodiment of the present invention.

Fig. 2 is the oscillogram of harmonic signal.

Fig. 3 is the amplitude versus frequency characte figure taking advantage of 210 window functions mutually.

Fig. 4 is the amplitude versus frequency characte figure taking advantage of 201 window functions mutually.

Fig. 5 is the amplitude versus frequency characte figure taking advantage of 120 window functions mutually.

Fig. 6 is the amplitude versus frequency characte figure taking advantage of 102 window functions mutually.

Fig. 7 is the amplitude versus frequency characte figure taking advantage of 111 window functions mutually.

Fig. 8 is the amplitude versus frequency characte figure taking advantage of 021 window function mutually.

Embodiment

Below in conjunction with drawings and the specific embodiments, the present invention is described in further detail.

Shown in Figure 1, the embodiment of the present invention provides a kind of three spectral line interpolation FFT harmonic analysis methods based on cross multiplication window function, comprises the following steps:

S1, structure cross multiplication window function:

The general formula of multiplication window function is produced by multiple window function product, and the general formula of multiplication window function w (n) is:

Wherein, n is the ordinal number of sampled point, and n is natural number; M is the number of the basic window function participating in structure multiplication window, and m is positive integer; w _in () is i-th basic window function expression formula, i=1,2 ... m; c _ifor participating in the number of i-th basic window function, be called w _ithe sub-order of (n); c is total order of cross multiplication window function.According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), namely construct multi-form multiplication window function.

Work as w _iwhen () expression formula is identical n, when namely adopted basic window is identical, w (n) is multiplication window function, works as w _iwhen () expression formula is different n, when namely adopted basic window is different, w (n) is cross multiplication window function.

Below for three conventional window functions, expression formula is respectively w ₁(n), w ₂(n), w ₃(n).Total order c=c ₁+ c ₂+ c ₃.Because c is larger, mean that the basic window function participating in multiplication window is more, the multiplication window function item number of formation is more, can increase the complexity of calculating.General to c maximal value, such as, order: c=3.

According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), construct following 9 kinds of multiplication window functions:

Work as c ₁=3, c ₂=0, c ₃when=0, construct and take advantage of 300 window functions mutually;

Work as c ₁=2, c ₂=1, c ₃when=0, construct and take advantage of 210 window functions mutually;

Work as c ₁=2, c ₂=0, c ₃when=1, construct and take advantage of 201 window functions mutually;

Work as c ₁=1, c ₂=2, c ₃when=0, construct and take advantage of 120 window functions mutually;

Work as c ₁=1, c ₂=0, c ₃when=2, construct and take advantage of 102 window functions mutually;

Work as c ₁=0, c ₂=3, c ₃when=0, construct and take advantage of 030 window function mutually;

Work as c ₁=1, c ₂=1, c ₃when=1, construct and take advantage of 111 window functions mutually;

Work as c ₁=0, c ₂=2, c ₃when=1, construct and take advantage of 021 window function mutually;

Work as c ₁=0, c ₂=0, c ₃when=3, construct and take advantage of 003 window function mutually;

As can be seen from above-mentioned 9 kinds of multiplication window functions: take advantage of 300 window functions mutually, take advantage of 030 window function mutually, take advantage of 003 window function mutually, in fact these three kinds of window functions belong to multiplication window function (the present invention will not discuss), all the other the 6 kinds cross multiplication window functions for constructing.

S2, Signal Pretreatment:

Mutual inductor gathers power network signal, and power network signal comprises current signal, voltage signal, and power network signal x (n) collected by mutual inductor, is transferred to host computer; Add cross multiplication window function w (n) to power network signal x (n) to block, obtain windowing signal x _w(n):

x _w(n)＝x(n)w(n) (2)

After FFT conversion is carried out to the windowing signal of formula (2), obtain windowing FFT frequency spectrum:

Wherein, the frequency spectrum that W () is window function, k is positive integer, and X (k) represents the frequency spectrum of kth subharmonic, and Ak is the amplitude of kth subharmonic, and j represents imaginary unit, and e is the truth of a matter of natural logarithm, for the initial phase of kth subharmonic, harmonic wave is first-harmonic for the first time, f _sfor sample frequency, f ₀for fundamental frequency, Δ f is discrete frequency intervals, and Δ f=f _s/ N.

Order: k ₀=f ₀/ Δ f, k ₀for the position of spectral line of real frequency spectrum, because non-synchronous sampling or non-integer-period block, k ₀being generally decimal, is not integer.

If ignore the impact of negative frequency point place secondary lobe, formula (3) becomes:

S3, determine three spectral lines:

At the windowing FFT spectrum peak near zone that S2 obtains, k ₀three spectral lines that place's Frequency point is larger are respectively: kth ₁, k ₂, k ₃root, k ₁, k ₂, k ₃be positive integer, k ₁~ k ₃pass be: k ₂=k ₁+ 1, k ₃=k ₂+ 1, the amplitude that these three spectral lines are corresponding is respectively y ₁, y ₂, y ₃.

For convenience of calculating, note variable α=k-k ₂, then-0.5≤α≤0.5;

Another note variable $\mathrm{\β}=\frac{y_3-y_1}{y_2}---\left(5\right)$

S4, calculate the correction formula of three spectral line interpolation algorithms:

Can obtain according to formula (4) and (5):

$\mathrm{\β}=\frac{\left|W(-\mathrm{\α}+1)\right|-\left|W(-\mathrm{\α}-1)\right|}{\left|W(-\mathrm{\α})\right|}---\left(6\right)$

When N value is larger, such as: during N > 1000, formula (6) can abbreviation be a function β=g (α), and its inverse function is α=g ^-1(β).Because adopted Cosine Window coefficient is real coefficient, its frequency response is even symmetry, thus g () and g ^-1() is odd function.

Polynomial Approximation Method is adopted to calculate odd function β=g ^-1(α), expression formula is:

α≈p ₁₁×β+p ₁₃×β ³+…+p _1pβ ^p(7)

P ₁₁, p ₁₃; P _1pfor the odd item coefficient of approximation by polynomi-als, p is odd number.

According to formula (4), try to achieve the amplitude Ai of power network signal i-th subharmonic:

$A_i=2y_i\left|W(k-\frac{f_0}{\mathrm{\Δf}})\right|---\left(8\right)$

Wherein, i is positive integer, y _ifor the amplitude of the i-th subharmonic after windowing FFT;

For first-harmonic, consider y ₂be from the nearest spectral line of true spectral line point, give larger weight, can obtain:

$A_1=\frac{2(y_3+2y_2+y_1)}{\left|W(-\mathrm{\α}+1)\right|+2\left|W(-\mathrm{\α})\right|+\left|W(-\mathrm{\α}-1)\right|}---\left(9\right)$

The approach method of similar formula (7), when N is larger, such as, during N > 1000, window function coefficient is real coefficient, and formula (9) can be expressed as:

A ₁＝N ^-1(y ₃+2y ₂+y ₁)u(α)

Wherein, u (α) is correction formula, and is even function, and approximating polynomial is not containing odd item.

Three spectral line correction approximating polynomials are as follows:

u(α)＝(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (10)

In formula (10), p ₂₀, p ₂₂p _2dfor the even term coefficient of approximation by polynomi-als, d is the highest order of matching, and d is even number.

S5, calculating first-harmonic parameter:

Consider y ₂be from two nearest spectral lines of true spectral line point, give larger weight, can fundamental frequency f be calculated ₀, fundamental voltage amplitude A ₁:

f ₀＝k·Δf＝(k ₂+α)Δf (11)

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (12)

According to formula (4), the phase place of first-harmonic can also be calculated:

Copy asking for of first-harmonic parameter, the analysis of each harmonic parameter can be carried out according to formula (6), (7), (9), (11), (12), (13).

S6, determine harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, step S3 ~ S5 is repeated, until all harmonic parameters calculate complete;

S7, carry out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compare with custom window function interpolation arithmetic accuracy.

Be described in detail below by a specific embodiment.The harmonic wave of a selected power network signal, the waveform of harmonic wave is shown in Figure 1.

The embodiment of the present invention provides a kind of three spectral line interpolation FFT harmonic analysis methods based on cross multiplication window function, specifically comprises the following steps:

The structure of step 1, cross multiplication window function:

For Hanning window, Hamming window, Blackman window, provide the structural model of cross multiplication window function, 6 kinds of multiplication window function performances of structure are as shown in table 2.Consider the problem of calculated amount, the order of multiplication window function is not suitable for too high, therefore limits order c=3.The sub-order of often kind of window function may get 0,1,2 three value.In order to the aspect of writing, simplify expression-form below: Hanning is abbreviated as Hn, Hamming is abbreviated as Hm, and Blackman is abbreviated as Bm.

Table 1, performance based on normal function structure cross multiplication window function

Hn(c 1)	Hm(c 2)	Bm(c 3)	Main lobe width (× π)	Sidelobe reduction (dB)
					2	1	0	0.0054	-62.8
2	0	1	0.0059	-84.5
					1	2	0	0.0051	-64.3
1	0	2	0.0061	-100.6
					1	1	1	0.0056	-84.9
0	2	1	0.0056	-86.1

The pre-service of step 2, signal:

Discrete power network signal x (n) collected by mutual inductor, signal model ginseng is shown in Table 2.In order to verify the precision of put forward algorithm, carry out 10 subharmonic simulation analysis.Signal model is:

Wherein: fundamental frequency f ₀for 50.5Hz, sample frequency f _sfor 5120Hz, sampling number N is 1024.

The harmonic parameters of table 2, power network signal

Signal spectrum is obtained after windowing FFT conversion:

Step 3, determine three spectral lines:

In 45 ~ 55Hz, seek maximum three spectral lines of mould, be respectively y ₁, y ₂, y ₃.

So $\mathrm{\β}=\frac{y_3-y_1}{y_2}$

Step 4, calculate the correction formula of three spectral line interpolation algorithms:

According to the derivation of step S4 in embodiment, the correction formula α=g of 4 kinds of window functions in table 2 ^-1(β), u (α) is as follows respectively:

(1) 210 window functions are taken advantage of mutually:

α＝1.11458715β-0.0908720β ³+0.01413650β ⁵

u(α)＝1.80783068+0.41014621α ²+0.05139583α ⁴

Take advantage of 210 window function amplitude versus frequency charactes shown in Figure 3 mutually.

(2) 201 window functions are taken advantage of mutually:

α＝1.28076537β-0.09843530β ³+0.01499840β ⁵

u(α)＝1.95369960+0.38423865α ²+0.04102745α ⁴

Take advantage of 201 window function amplitude versus frequency charactes shown in Figure 4 mutually.

(3) 120 window functions are taken advantage of mutually:

α＝1.08516300β-0.08842400β ³+0.0140181β ⁵

u(α)＝1.78667431+0.41627357α ²+0.05339022α ⁴

Take advantage of the amplitude versus frequency characte of 120 window functions shown in Figure 5 mutually.

(4) 102 window functions are taken advantage of mutually:

α＝1.42354957β-0.10399918β ³+0.01605665β ⁵

u(α)＝2.07275558+0.36590106α ²+0.03458435α ⁴

Take advantage of the amplitude versus frequency characte of 102 window functions shown in Figure 6 mutually.

(5) 111 window functions are taken advantage of mutually:

α＝1.2530640β-0.09572246β ³+0.01451709β ⁵

u(α)＝1.93442536+0.38883257α ²+0.04234880α ⁴

Take advantage of the amplitude versus frequency characte of 111 window functions shown in Figure 7 mutually.

(6) 021 window function is taken advantage of mutually:

α＝1.22447178β-0.09322158β ³+0.01445149β ⁵

u(α)＝1.91516084+0.39387119α ²+0.04376915α ⁴

Take advantage of the amplitude versus frequency characte of 021 window function shown in Figure 8 mutually.

Step 5, calculating first-harmonic parameter:

f ₀＝k·Δf＝(k ₂+α)Δf

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d)

Step 6, determine harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, step 3 ~ 5 are repeated, until all harmonic parameters calculate complete.

Step 7, carry out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compare with custom window function interpolation arithmetic accuracy.Error analysis comparative result see shown in table 3 ~ table 5, wherein, D _airepresent the relative error of first-harmonic and each harmonic amplitude measure; D _f0represent the relative error of fundamental wave frequency measurement value; represent the relative error of first-harmonic and each harmonic initial phase measured value, all represent with number percent.

Table 3, cross multiplication window amplitude, frequency measurement relative error

Relatively item

D A1

D A2

D A3

D A4

D A5

210

8.05E-7

1.25E-6

-4.83E-7

7.27E-7

-2.42E-6

201

6.54E-7

2.70E-8

-3.15E-7

6.81E-8

-1.54E-6

120

8.38E-7

-8.76E-7

-1.87E-7

-1.18E-7

-6.22E-6

102

5.58E-7

-6.91E-7

8.36E-7

8.71E-7

-3.98E-6

111

6.67E-7

-2.11E-7

1.57E-6

6.91E-9

-4.82E-6

021

6.96E-7

-8.10E-7

-1.60E-6

-1.23E-7

-5.02E-6

Relatively item

D A6

D A71

D A8

D A9

D A10

D f0

210

7.56E-7

-4.47E-7

-4.13E-7

8.38E-7

1.61E-6

-1.20E-8

201

-1.44E-8

-2.37E-7

-5.84E-8

7.02E-7

1.20E-6

-2.78E-9

120

-1.25E-6

-2.18E-6

-6.44E-8

9.73E-7

-1.80E-6

-3.14E-9

102

1.18E-5

1.40E-6

-3.24E-7

5.06E-7

8.57E-7

-1.82E-8

111

-1.15E-5

-2.05E-6

-2.66E-8

7.24E-7

1.30E-6

6.77E-9

021

-1.18E-5

-2.05E-6

1.33E-8

7.73E-7

1.36E-6

6.44E-9

Table 4, cross multiplication window phase measurement relative error

And adopt separately Hanning, Hamming and Blackman window function amplitude, frequency and phase measurement error as shown in table 5.

The amplitude of table 5, Hn, Hm, Bm window, frequency and phase place relative error

As can be seen here, the cross multiplication window function in the embodiment of the present invention is in basic frequency, and 1 ~ 3 order of magnitude higher than basic window function, count high 2 ~ 3 orders of magnitude at crest meter, in phase calculation, advantage is more outstanding, high 3 ~ 6 orders of magnitude.As can be seen from above simulation result, result of calculation is generally better than general window function spectral line interpolation algorithm.

In sum, in the embodiment of the present invention based on cross multiplication window function harmonic analysis method, realize structure cross multiplication window function first, cross multiplication window function be applied in three spectral line interpolation FFT algorithms, calculate harmonic parameters, the precision higher than custom window function can be obtained.

The embodiment of the present invention also provides a kind of three spectral line interpolation FFT frequency analysis systems based on cross multiplication window function, comprise cross multiplication window function tectonic element, Signal Pretreatment unit, spectral line determining unit, correction formula computing unit, first-harmonic parameter calculation unit, harmonic parameters determining unit, error analysis unit, wherein:

Cross multiplication construction of function unit, for constructing cross multiplication window function:

The general formula of multiplication window function is produced by multiple window function product, and the general formula of multiplication window function w (n) is:

Wherein, n is the ordinal number of sampled point, and n is natural number; M is the number of the basic window function participating in structure multiplication window, and m is positive integer; w _in () is i-th basic window function expression formula, i=1,2 ... m; c _ifor participating in the number of i-th basic window function, be called w _ithe sub-order of (n); c is total order of cross multiplication window function; According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), namely construct multi-form multiplication window function;

Work as w _iwhen () expression formula is identical n, when namely adopted basic window is identical, w (n) is multiplication window function, works as w _iwhen () expression formula is different n, when namely adopted basic window is different, w (n) is cross multiplication window function;

Below for three conventional window functions, expression formula is respectively w ₁(n), w ₂(n), w ₃(n).Total order c=c ₁+ c ₂+ c ₃.Because c is larger, mean that the basic window function participating in multiplication window is more, the multiplication window function item number of formation is more, can increase the complexity of calculating.General to c maximal value, such as, order: c=3.

According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), construct following 9 kinds of multiplication window functions:

Work as c ₁=3, c ₂=0, c ₃when=0, construct and take advantage of 300 window functions mutually;

Work as c ₁=2, c ₂=1, c ₃when=0, construct and take advantage of 210 window functions mutually;

Work as c ₁=2, c ₂=0, c ₃when=1, construct and take advantage of 201 window functions mutually;

Work as c ₁=1, c ₂=2, c ₃when=0, construct and take advantage of 120 window functions mutually;

Work as c ₁=1, c ₂=0, c ₃when=2, construct and take advantage of 102 window functions mutually;

Work as c ₁=0, c ₂=3, c ₃when=0, construct and take advantage of 030 window function mutually;

Work as c ₁=1, c ₂=1, c ₃when=1, construct and take advantage of 111 window functions mutually;

Work as c ₁=0, c ₂=2, c ₃when=1, construct and take advantage of 021 window function mutually;

Work as c ₁=0, c ₂=0, c ₃when=3, construct and take advantage of 003 window function mutually;

As can be seen from above-mentioned 9 kinds of multiplication window functions: take advantage of 300 window functions mutually, take advantage of 030 window function mutually, take advantage of 003 window function mutually, in fact these three kinds of window functions belong to multiplication window function (the present invention will not discuss), all the other the 6 kinds cross multiplication window functions for constructing.

Signal Pretreatment unit, for carrying out pre-service to signal:

Mutual inductor gathers power network signal, and power network signal comprises current signal, voltage signal; Power network signal x (n) collected by mutual inductor, is transferred to host computer; Add cross multiplication window function w (n) to power network signal x (n) to block, obtain windowing signal x _w(n):

x _w(n)＝x(n)w(n) (2)

After FFT conversion is carried out to the windowing signal of formula (2), obtain windowing FFT frequency spectrum:

Wherein, the frequency spectrum that W () is window function, k is positive integer, and X (k) represents the frequency spectrum of kth subharmonic, and Ak is the amplitude of kth subharmonic, and j represents imaginary unit, and e is the truth of a matter of natural logarithm, for the initial phase of kth subharmonic, harmonic wave is first-harmonic for the first time, f _sfor sample frequency, f ₀for fundamental frequency, Δ f is discrete frequency intervals, and Δ f=f _s/ N;

Order: k ₀=f ₀/ Δ f, k ₀for the position of spectral line of real frequency spectrum, because non-synchronous sampling or non-integer-period block, k ₀being generally decimal, is not integer.

Ignore the impact of negative frequency point place secondary lobe, formula (3) becomes:

Spectral line determining unit, for determining three spectral lines:

At the windowing FFT spectrum peak near zone that S2 obtains, k ₀three spectral lines that place's Frequency point is larger are respectively: kth ₁, k ₂, k ₃root, k ₁, k ₂, k ₃be positive integer, k ₁~ k ₃pass be: k ₂=k ₁+ 1, k ₃=k ₂+ 1, the amplitude that these three spectral lines are corresponding is respectively y ₁, y ₂, y ₃;

Note variable α=k-k ₂, then-0.5≤α≤0.5;

Another note variable $\mathrm{\β}=\frac{y_3-y_1}{y_2}---\left(5\right);$

Correction formula computing unit, for calculating the correction formula of three spectral line interpolation algorithms:

Obtain according to formula (4) and (5):

$\mathrm{\β}=\frac{\left|W(-\mathrm{\α}+1)\right|-\left|W(-\mathrm{\α}-1)\right|}{\left|W(-\mathrm{\α})\right|}---\left(6\right)$

When N value is larger, such as: during N > 1000, formula (6) can abbreviation be a function β=g (α), and its inverse function is α=g ^-1(β).Because adopted Cosine Window coefficient is real coefficient, its frequency response is even symmetry, thus g () and g ^-1() is odd function.

Polynomial Approximation Method is adopted to calculate odd function β=g ^-1(α), expression formula is:

α≈p ₁₁×β+p ₁₃×β ³+…+p _1pβ ^p(7)

P ₁₁, p ₁₃; P _1pfor the odd item coefficient of approximation by polynomi-als, p is odd number;

According to formula (4), try to achieve the amplitude Ai of power network signal i-th subharmonic:

$A_i=2y_i\left|W(k-\frac{f_0}{\mathrm{\Δf}})\right|---\left(8\right)$

Wherein, i is positive integer, y _ifor the amplitude of the i-th subharmonic after windowing FFT;

For first-harmonic, consider y ₂be from the nearest spectral line of true spectral line point, give larger weight, can obtain:

$A_1=\frac{2(y_3+2y_2+y_1)}{\left|W(-\mathrm{\α}+1)\right|+2\left|W(-\mathrm{\α})\right|+\left|W(-\mathrm{\α}-1)\right|}---\left(9\right)$

During N > 1000, window function coefficient is real coefficient, and formula (9) is expressed as:

A ₁＝N ^-1(y ₃+2y ₂+y ₁)u(α)

Wherein, u (α) is correction formula, and is even function, and approximating polynomial is not containing odd item;

Three spectral line correction approximating polynomials are as follows:

u(α)＝(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (10)

In formula (10), p ₂₀, p ₂₂p _2dfor the even term coefficient of approximation by polynomi-als, d is the highest order of matching, and d is even number;

First-harmonic parameter calculation unit, for calculating first-harmonic parameter:

Consider y ₂be from two nearest spectral lines of true spectral line point, give larger weight, can fundamental frequency f be calculated ₀, fundamental voltage amplitude A ₁:

f ₀＝k·Δf＝(k ₂+α)Δf (11)

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (12)

According to formula (4), calculate the phase place of first-harmonic:

Copy asking for of first-harmonic parameter, carry out the analysis of each harmonic parameter according to formula (6), (7), (9), (11), (12), (13);

Harmonic parameters determining unit, for determining harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, described spectral line determining unit, correction formula computing unit, first-harmonic parameter calculation unit repeat to calculate, until all harmonic parameters calculate complete;

Error analysis unit, for carrying out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compares with custom window function interpolation arithmetic accuracy.

Those skilled in the art can carry out various modifications and variations to the embodiment of the present invention, if these amendments and modification are within the scope of the claims in the present invention and equivalent technologies thereof, then these revise and modification also within protection scope of the present invention.

The prior art that the content do not described in detail in instructions is known to the skilled person.

Claims (6)

1., based on three spectral line interpolation FFT harmonic analysis methods of cross multiplication window function, it is characterized in that, comprise the following steps:

S1, structure cross multiplication window function:

The general formula of multiplication window function is produced by multiple window function product, and the general formula of multiplication window function w (n) is:

Wherein, n is the ordinal number of sampled point, and n is natural number; M is the number of the basic window function participating in structure multiplication window, and m is positive integer; w _in () is i-th basic window function expression formula, i=1,2 ... m; c _ifor participating in the number of i-th basic window function, be called w _ithe sub-order of (n); c is total order of cross multiplication window function; According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), namely construct multi-form multiplication window function;

Work as w _iwhen () expression formula is identical n, when namely adopted basic window is identical, w (n) is multiplication window function, works as w _iwhen () expression formula is different n, when namely adopted basic window is different, w (n) is cross multiplication window function;

The expression formula of three conventional window functions is respectively w ₁(n), w ₂(n), w ₃(n), total order c=c ₁+ c ₂+ c ₃, to c maximal value, order: c=3;

According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), construct following 9 kinds of multiplication window functions:

Work as c ₁=3, c ₂=0, c ₃when=0, construct and take advantage of 300 window functions mutually;

Work as c ₁=2, c ₂=1, c ₃when=0, construct and take advantage of 210 window functions mutually;

Work as c ₁=2, c ₂=0, c ₃when=1, construct and take advantage of 201 window functions mutually;

Work as c ₁=1, c ₂=2, c ₃when=0, construct and take advantage of 120 window functions mutually;

Work as c ₁=1, c ₂=0, c ₃when=2, construct and take advantage of 102 window functions mutually;

Work as c ₁=0, c ₂=3, c ₃when=0, construct and take advantage of 030 window function mutually;

Work as c ₁=1, c ₂=1, c ₃when=1, construct and take advantage of 111 window functions mutually;

Work as c ₁=0, c ₂=2, c ₃when=1, construct and take advantage of 021 window function mutually;

Work as c ₁=0, c ₂=0, c ₃when=3, construct and take advantage of 003 window function mutually;

Wherein, take advantage of 300 window functions mutually, take advantage of 030 window function mutually, take advantage of 003 window function mutually, these three kinds of window functions belong to multiplication window function, all the other the 6 kinds cross multiplication window functions for constructing;

S2, Signal Pretreatment:

Mutual inductor gathers power network signal, and power network signal x (n) collected by mutual inductor, is transferred to host computer; Add cross multiplication window function w (n) to power network signal x (n) to block, obtain windowing signal x _w(n):

x _w(n)＝x(n)w(n) (2)

After FFT conversion is carried out to the windowing signal of formula (2), obtain windowing FFT frequency spectrum:

Wherein, the frequency spectrum that W () is window function, k is positive integer, and X (k) represents the frequency spectrum of kth subharmonic, A _kfor the amplitude of kth subharmonic, j represents imaginary unit, and e is the truth of a matter of natural logarithm, for the initial phase of kth subharmonic, harmonic wave is first-harmonic for the first time, f _sfor sample frequency, f ₀for fundamental frequency, Δ f is discrete frequency intervals, and Δ f=f _s/ N;

Order: k ₀=f ₀/ Δ f, k ₀for the position of spectral line of real frequency spectrum,

Ignore the impact of negative frequency point place secondary lobe, formula (3) becomes:

S3, determine three spectral lines:

At the windowing FFT spectrum peak near zone that S2 obtains, k ₀three spectral lines that place's Frequency point is larger are respectively: kth ₁, k ₂, k ₃root, k ₁, k ₂, k ₃be positive integer, k ₁~ k ₃pass be: k ₂=k ₁+ 1, k ₃=k ₂+ 1, the amplitude that these three spectral lines are corresponding is respectively y ₁, y ₂, y ₃;

Note variable α=k-k ₂, then-0.5≤α≤0.5;

Another note variable $\mathrm{\β}=\frac{y_3-y_1}{y_2}---\left(5\right)$

S4, calculate the correction formula of three spectral line interpolation algorithms:

Obtain according to formula (4) and (5):

$\mathrm{\β}=\frac{\left|W(-\mathrm{\α}+1)\right|-\left|W(-\mathrm{\α}-1)\right|}{\left|W(-\mathrm{\α})\right|}---\left(6\right)$

Polynomial Approximation Method is adopted to calculate odd function β=g ^-1(α), expression formula is:

α≈p ₁₁×β+p ₁₃×β ³+…+p _1pβ ^p(7)

P ₁₁, p ₁₃; P _1pfor the odd item coefficient of approximation by polynomi-als, p is odd number;

According to formula (4), try to achieve the amplitude A of power network signal i-th subharmonic _i:

$A_i=2y_i\left|W(k-\frac{f_0}{\mathrm{\Δf}})\right|---\left(8\right)$

Wherein, i is positive integer, y _ifor the amplitude of the i-th subharmonic after windowing FFT;

Consider y ₂be from the nearest spectral line of true spectral line point, obtain:

$A_1=\frac{2(y_3+2y_2+y_1)}{\left|W(-\mathrm{\α}+1)\right|+2\left|W(-\mathrm{\α})\right|+\left|W(-\mathrm{\α}-1)\right|}---\left(9\right)$

During N > 1000, window function coefficient is real coefficient, and formula (9) is expressed as:

A ₁＝N ^-1(y ₃+2y ₂+y ₁)u(α)

Wherein, u (α) is correction formula, and is even function, and approximating polynomial is not containing odd item;

Three spectral line correction approximating polynomials are as follows:

u(α)＝(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (10)

In formula (10), p ₂₀, p ₂₂p _2dfor the even term coefficient of approximation by polynomi-als, d is the highest order of matching, and d is even number;

S5, calculating first-harmonic parameter:

Calculate fundamental frequency f ₀, fundamental voltage amplitude A ₁:

f ₀＝k·Δf＝(k ₂+α)Δf (11)

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (12)

According to formula (4), calculate the phase place of first-harmonic:

Copy asking for of first-harmonic parameter, carry out the analysis of each harmonic parameter according to formula (6), (7), (9), (11), (12), (13);

S6, determine harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, step S3 ~ S5 is repeated, until all harmonic parameters calculate complete;

S7, carry out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compare with custom window function interpolation arithmetic accuracy.

2., as claimed in claim 1 based on three spectral line interpolation FFT harmonic analysis methods of cross multiplication window function, it is characterized in that: described power network signal comprises current signal, voltage signal.

3., as claimed in claim 1 based on three spectral line interpolation FFT harmonic analysis methods of cross multiplication window function, it is characterized in that: the position of spectral line k of described real frequency spectrum ₀for decimal.

4. three spectral line interpolation FFT frequency analysis systems based on cross multiplication window function, it is characterized in that: comprise cross multiplication window function tectonic element, Signal Pretreatment unit, spectral line determining unit, correction formula computing unit, first-harmonic parameter calculation unit, harmonic parameters determining unit, error analysis unit, wherein:

Described cross multiplication construction of function unit, for constructing cross multiplication window function:

The general formula of multiplication window function is produced by multiple window function product, and the general formula of multiplication window function w (n) is:

Wherein, n is the ordinal number of sampled point, and n is natural number; M is the number of the basic window function participating in structure multiplication window, and m is positive integer; w _in () is i-th basic window function expression formula, i=1,2 ... m; c _ifor participating in the number of i-th basic window function, be called w _ithe sub-order of (n); c is total order of cross multiplication window function; According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), namely construct multi-form multiplication window function;

Work as w _iwhen () expression formula is identical n, when namely adopted basic window is identical, w (n) is multiplication window function, works as w _iwhen () expression formula is different n, when namely adopted basic window is different, w (n) is cross multiplication window function;

The expression formula of three conventional window functions is respectively w ₁(n), w ₂(n), w ₃(n), total order c=c ₁+ c ₂+ c ₃, to c maximal value, order: c=3;

According to c ₁, c ₂, c ₃different values, form different combination (c ₁c ₂c _m), construct following 9 kinds of multiplication window functions:

Work as c ₁=3, c ₂=0, c ₃when=0, construct and take advantage of 300 window functions mutually;

Work as c ₁=2, c ₂=1, c ₃when=0, construct and take advantage of 210 window functions mutually;

Work as c ₁=2, c ₂=0, c ₃when=1, construct and take advantage of 201 window functions mutually;

Work as c ₁=1, c ₂=2, c ₃when=0, construct and take advantage of 120 window functions mutually;

Work as c ₁=1, c ₂=0, c ₃when=2, construct and take advantage of 102 window functions mutually;

Work as c ₁=0, c ₂=3, c ₃when=0, construct and take advantage of 030 window function mutually;

Work as c ₁=1, c ₂=1, c ₃when=1, construct and take advantage of 111 window functions mutually;

Work as c ₁=0, c ₂=2, c ₃when=1, construct and take advantage of 021 window function mutually;

Work as c ₁=0, c ₂=0, c ₃when=3, construct and take advantage of 003 window function mutually;

Wherein, take advantage of 300 window functions mutually, take advantage of 030 window function mutually, take advantage of 003 window function mutually, these three kinds of window functions belong to multiplication window function, all the other the 6 kinds cross multiplication window functions for constructing;

Described Signal Pretreatment unit, for carrying out pre-service to signal:

Mutual inductor gathers power network signal, and power network signal x (n) collected by mutual inductor, is transferred to host computer; Add cross multiplication window function w (n) to power network signal x (n) to block, obtain windowing signal x _w(n):

x _w(n)＝x(n)w(n) (2)

After FFT conversion is carried out to the windowing signal of formula (2), obtain windowing FFT frequency spectrum:

Wherein, the frequency spectrum that W () is window function, k is positive integer, and X (k) represents the frequency spectrum of kth subharmonic, A _kfor the amplitude of kth subharmonic, j represents imaginary unit, and e is the truth of a matter of natural logarithm, for the initial phase of kth subharmonic, harmonic wave is first-harmonic for the first time, f _sfor sample frequency, f ₀for fundamental frequency, Δ f is discrete frequency intervals, and Δ f=f _s/ N;

Order: k ₀=f ₀/ Δ f, k ₀for the position of spectral line of real frequency spectrum,

Ignore the impact of negative frequency point place secondary lobe, formula (3) becomes:

Described spectral line determining unit, for determining three spectral lines:

At the windowing FFT spectrum peak near zone that S2 obtains, k ₀three spectral lines that place's Frequency point is larger are respectively: kth ₁, k ₂, k ₃root, k ₁, k ₂, k ₃be positive integer, k ₁~ k ₃pass be: k ₂=k ₁+ 1, k ₃=k ₂+ 1, the amplitude that these three spectral lines are corresponding is respectively y ₁, y ₂, y ₃;

Note variable α=k-k ₂, then-0.5≤α≤0.5;

Another note variable $\mathrm{\β}=\frac{y_3-y_1}{y_2}---\left(5\right);$

Described correction formula computing unit, for calculating the correction formula of three spectral line interpolation algorithms:

Obtain according to formula (4) and (5):

$\mathrm{\β}=\frac{\left|W(-\mathrm{\α}+1)\right|-\left|W(-\mathrm{\α}-1)\right|}{\left|W(-\mathrm{\α})\right|}---\left(6\right)$

Polynomial Approximation Method is adopted to calculate odd function β=g ^-1(α), expression formula is:

α≈p ₁₁×β+p ₁₃×β ³+…+p _1pβ ^p(7)

P ₁₁, p ₁₃; P _1pfor the odd item coefficient of approximation by polynomi-als, p is odd number;

According to formula (4), try to achieve the amplitude A of power network signal i-th subharmonic _i:

$A_i=2y_i\left|W(k-\frac{f_0}{\mathrm{\Δf}})\right|---\left(8\right)$

Wherein, i is positive integer, y _ifor the amplitude of the i-th subharmonic after windowing FFT;

Consider y ₂be from the nearest spectral line of true spectral line point, obtain:

$A_1=\frac{2(y_3+2y_2+y_1)}{\left|W(-\mathrm{\α}+1)\right|+2\left|W(-\mathrm{\α})\right|+\left|W(-\mathrm{\α}-1)\right|}---\left(9\right)$

During N > 1000, window function coefficient is real coefficient, and formula (9) is expressed as:

A ₁＝N ^-1(y ₃+2y ₂+y ₁)u(α)

Wherein, u (α) is correction formula, and is even function, and approximating polynomial is not containing odd item;

Three spectral line correction approximating polynomials are as follows:

u(α)＝(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (10)

In formula (10), p ₂₀, p ₂₂p _2dfor the even term coefficient of approximation by polynomi-als, d is the highest order of matching, and d is even number;

Described first-harmonic parameter calculation unit, for calculating first-harmonic parameter:

Calculate fundamental frequency f ₀, fundamental voltage amplitude A ₁:

f ₀＝k·Δf＝(k ₂+α)Δf (11)

A ₁＝N ^-1(y ₃+2y ₂+y ₁)(p ₂₀+p ₂₂α ²+…+p _2dα ^d) (12)

According to formula (4), calculate the phase place of first-harmonic:

Copy asking for of first-harmonic parameter, carry out the analysis of each harmonic parameter according to formula (6), (7), (9), (11), (12), (13);

Described harmonic parameters determining unit, for determining harmonic parameters: determine fundamental frequency f ₀after, at scope (kf ₀-5, kf ₀+ 5) in, described spectral line determining unit, correction formula computing unit, first-harmonic parameter calculation unit repeat to calculate, until all harmonic parameters calculate complete;

Described error analysis unit, for carrying out error analysis: the error analyzing the three spectral line interpolation FFT algorithms based on cross multiplication window function, and compares with custom window function interpolation arithmetic accuracy.

5., as claimed in claim 4 based on three spectral line interpolation FFT frequency analysis systems of cross multiplication window function, it is characterized in that: described power network signal comprises current signal, voltage signal.

6., as claimed in claim 4 based on three spectral line interpolation FFT frequency analysis systems of cross multiplication window function, it is characterized in that: the position of spectral line k of described real frequency spectrum ₀for decimal.