CN104459315A - Inter-harmonic detection method based on non-base 2FFT transformation - Google Patents

Abstract

The invention discloses an inter-harmonic detection method based on non-base 2FFT transformation. When the number N of FFT points is equal to 2<M>*N', N1 is made to be equal to 2<M>, N2 is made to be equal to N', the FFT operation data are divided into N2 groups, each group has N1 point data, the data serial number n extracted from original FFT operation data is equal to N2n1+n2, 2FFT transformation is carried out on the N1 point data in each group, the k is equal to k1+N1k2, and the kth FFT transformation result in the original FFT operation data is obtained according to the calculation formula. The new non-base 2FFT algorithm is provided, FFT transformation of data with the point number N equal to 2<M>*N' is achieved, no interpolation zero-padding is needed, inter-harmonic detection precision can be improved, the frequency resolution meets national standard requirements, a transformation result needing to be calculated can be selected according to actual requirements, and operation time is saved.

Description

Based on the harmonic detection method that non-base 2FFT converts

Technical field

The invention belongs to m-Acetyl chlorophosphonazo detection technique field, more specifically say, relate to a kind of harmonic detection method converted based on non-base 2FFT.

Background technology

Along with the fast development of science and technology, the widespread use of the power equipments such as a large amount of impact load, the change of current, cause harmonic pollution in electric power net day by day serious, greatly have impact on the safety of electric system, especially along with various complex precise, constantly universal to the consumer of quality of power supply sensitivity, fractional harmoni (m-Acetyl chlorophosphonazo) also becomes the object that engineering and research field are constantly paid close attention to day by day.M-Acetyl chlorophosphonazo can cause the noise of the rolling of voltage flicker, TV images and radio receiver or other audio frequency apparatuses, also motor can be made to produce very strong noise and vibration, cause waveform distortion of the power supply network, thus reduce the power factor of load, increase various energy loss, for ensureing power grid quality and safety, need to strengthen the accurate detection to m-Acetyl chlorophosphonazo, so that the improvement of the electrical network quality of power supply.

Current most widely used harmonic measuring method is Fourier transform, and Fourier transform can obtain amplitude and the phase place of harmonic wave simultaneously, and function is more, and hardware implementing is simple.The Fourier transformation method of m-Acetyl chlorophosphonazo is: sample sequence extended to contiguous 2 by interpolation zero padding ^mpoint, then uses 2-base algorithm to carry out 2 ^mpoint FFT (Fast Fourier Transform, fast fourier transform).Below the base 2FFT algorithm that frequency domain extracts is introduced.

Defined by DFT (Discrete Fourier Transform, discrete Fourier transformation):

$F\left[k\right]=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}f\left[n\right]W_N^{\mathrm{kn}},0\&le;k\&le;N-1---\left(1\right)$

Wherein, N represents and counts, and f [n] represents discrete signal, represent twiddle factor.

F [n] is split into two halves, i.e. f [n] and f [n+N/2] (0≤n≤N/2-1), (1) formula of substitution obtains

$F\left[k\right]=\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}f\left[n\right]W_N^{\mathrm{kn}}+\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}f[n+N/2]W_N^{k(n+N/2)},0\&le;k\&le;N-1---\left(2\right)$

Due to $W_N^{k(n+N/2)}=W_N^{\mathrm{kn}}{W}_N^{\mathrm{kN}/2}={(-1)}^k{W}_N^{\mathrm{kn}},$ (2) can merge into again for two in formula:

$F\left[k\right]=\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}\left\{f\right[n]+{(-1)}^{k}f[n+N/2\left]\right\}{W}_N^{\mathrm{kn}},0\&le;k\&le;N-1---\left(3\right)$

When k=2r is even number, notice (-1) ^k=1, $W_N^{\mathrm{kn}}=W_N^{2\mathrm{rn}}=e^{-i2\mathrm{\&pi;}2\mathrm{rn}/N}=W_{N/2}^{\mathrm{rn}},$ (3) formula becomes:

$\begin{array}{c}F\left[2r\right]=\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}\left(f\right[n]+f[n+N/2\left]\right)W_{N/2}^{\mathrm{rn}}\\ =\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}g\left(n\right)W_{N/2}^{\mathrm{rn}}\\ =G\left(r\right),0\&le;r\&le;N/2-1\end{array}---\left(4\right)$

When k=2r+1 is odd number, $W_N^{\mathrm{kn}}=W_N^{(2r+1)n}=e^{-i2\mathrm{\&pi;}(2r+1)n/N}=W_N^n{W}_{N/2}^{\mathrm{rn}},$ (3) formula becomes:

$\begin{array}{c}F[2r+1]=\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}\left\{\left(f\right[n]-f[n+N/2\left]\right)W_N^n\right\}{W}_{N/2}^{\mathrm{rn}}\\ =\underset{n=0}{\overset{N/2-1}{\mathrm{\&Sigma;}}}p\left(n\right)W_{N/2}^{\mathrm{rn}}=p\left(r\right),0\&le;r\&le;N/2-1\end{array}---\left(5\right)$

So just the N point DFT (F [k]) of a N point sequence (f [n]) is calculated the N/2 point DFT (G [r] and P [r]) having changed into two N/2 point sequences (g [n] and p [n]) to calculate.The calculated amount being divided into g [n] and p [n] by f [n] be N number of add (i.e. f [n]+f [n+N/2] and f [n]-f [n+N/2], 0≤n≤N/2-1) and N/2 take advantage of (namely ( 0≤n≤N/2-1).

Due to the N/2 point DFT that g [n] calculates, that in the N point DFT (F [k]) of f [n], k is that half of even number, that half of to be then k the be even number calculated by p [n], therefore need the F of even number k [k] DFT (G (r)) put together as g [n] that releases to export, while, exports the F of odd number k [k] DFT (P (r)) put together as p [n] that releases.Because k is frequency domain variable, therefore this algorithm is called the fft algorithm that frequency domain extracts.

Then, two N/2 point DFT still respectively can take advantage of N/2 to add through N/4 with said method and be divided into two N/4 point DFT (simultaneously also will do corresponding frequency domain to extract), thus are divided into 4 N/2 point DFT altogether, and total division calculated amount is still and N number ofly adds take advantage of individual with N/2.Such division can be gone on doing step by step, is not difficult to find out, the total division calculated amount often walked is all take advantage of with N/2 N number of adding.

Work as N=2 ^m, after the division of M-1 step, be just divided into the computational problem of N/2 following 2 DFT

$\left.\begin{array}{c}A=a{W}_2^{0\&CenterDot;0}+b{W}_2^{0\&CenterDot;1}=a+b\\ B=a{W}_2^{1\&CenterDot;0}+b{W}_2^{1\&CenterDot;1}=(a-b)\&CenterDot;W_2^0\end{array}\right\}---\left(6\right)$

Fig. 1 is that 8 FFT that frequency domain extracts calculate schematic diagram.According to Fig. 1,8 FFT are divided into 42 DFT through 2 steps.(6) calculated amount needed for formula is 2 and adds and take advantage of with 1, takes advantage of with N/2 so the amount of calculation completing N/2 2 DFT is still N number of adding.Thus complete the amount of calculation M × N=Nlog of whole N point DFT ₂n number ofly to add and M × N/2=(N/2) log ₂n number ofly take advantage of, this is than the N needed for directly calculating by definition ²individual take advantage of and add want much less.Such as, N=2 ¹⁰=1024, then M=10, the multiplication number needed for calculating with base 2FFT algorithm is M × N/2=5 × 1024, and the multiplication number needed for directly calculating by definition is N ²=1024 × 1024, the two differs 1024/5 ≈ 200 times, and N is larger, and the efficiency of FFT improves more obvious.

Generally, owing to having done the extraction of M-1 time point of odd even, the F [k] that the last N/2 of this algorithm 2 DFT calculate has not been that order extracts.The change of order can illustrate with binary code: first time extract the odd even of dividing be 1 or 0 to distinguish by binary code the 1st, it is even number when this position is 0, it is odd number when this position is 1, it is 1 or 0 to distinguish that second time smokes odd even by binary code the 2nd ... each extraction is all (left side) limit before even item is placed on, (right side) limit after odd term is placed on, thus the binary code of number is arranged in order from left to right according to binary digit after extracting, just in time contrary with the rule that ordinary binary number is arranged in order from right to left, so be called bit reversed order.After calculating F [k], bit reversed order to be become order.

In m-Acetyl chlorophosphonazo detects, owing to will calculate the amplitude of m-Acetyl chlorophosphonazo, necessary several cycles of data processing, because base 2FFT algorithm will reach 2 ^mpoint, during general sampling, the sampling number of each signal period is the power of 2, but data processing cycle number differ be decided to be 2 power, therefore some value places of m-Acetyl chlorophosphonazo may obtain for interpolation, cause measuring accuracy inaccurate.

Summary of the invention

The object of the invention is to overcome the deficiencies in the prior art, a kind of harmonic detection method converted based on non-base 2FFT is provided, realize the FFT conversion of non-base 2, do not need zero padding, improve the accuracy that m-Acetyl chlorophosphonazo detects, make testing result can meet IEC standard completely, when first-harmonic be 50Hz, sample window width makes frequency resolution can obtain exact value at 5Hz place when being 10 cycle.

For achieving the above object, the present invention is based on the harmonic detection method that non-base 2FFT converts, comprise the following steps:

S1: sample to signal, arranges sampling number N=2 ^m× N ', makes N ₁=2 ^m, N ₂=N ', is divided into N by sampled data ₂group, often organizes N ₁point data, n-th ₂in group n-th ₁point data is designated as f [n ₁, n ₂], the data sequence number n=N extracted in sampled data ₂n ₁+ n ₂, wherein n ₁span be 0≤n ₁≤ N ₁-1, n ₂span be 0≤n ₂≤ N ₂-1;

S2: respectively to the N in every group ₁point data carries out base 2FFT conversion, by n-th ₂kth in group ₁individual transformation results is designated as F ' [k ₁, n ₂], k ₁span be 0≤k ₁≤ N ₁-1;

S3: make k=k ₁+ N ₁k ₂, k ₂span be 0≤k ₂≤ N ₂-1, a kth Fourier transform result F [k of calculating sampling data ₁, k ₂] and export as testing result, computing formula is:

$F[k_1,k_2]=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{W}_{N_2}^{\&prime;n_2{k}_2}{F}^{\&prime;}[k_1,n_2]$

Wherein, $W_{N_2}^{\&prime;n_2{k}_2}=W_{n_2}^{n_2{k}_2}{W}_N^{n_2{k}_1}.$

The present invention is based on the harmonic detection method that non-base 2FFT converts, when FFT points N=2 ^m× N ', makes N ₁=2 ^m, N ₂=N ', is divided into N by FFT operational data ₂group, often organizes N ₁point data, the data sequence number n=N extracted in former FFT operational data ₂n ₁+ n ₂, respectively to the N in every group ₁point data carries out base 2FFT conversion, then k=k ₁+ N ₁k ₂, obtain a kth FFT transformation results in former FFT operational data according to computing formula.The present invention proposes a kind of non-base 2FFT algorithm newly, realize points N=2 ^mthe FFT conversion of the data of × N ', does not need interpolation zero padding, can improve the accuracy of detection of m-Acetyl chlorophosphonazo, frequency resolution is up to state standards requirement, and can selects calculative transformation results according to actual needs, save operation time.

Accompanying drawing explanation

Fig. 1 is that 8 FFT that frequency domain extracts calculate schematic diagram;

Fig. 2 is the embodiment process flow diagram that the present invention is based on the harmonic detection method that non-base 2FFT converts;

Fig. 3 is the result exemplary plot adopting the present invention to carry out m-Acetyl chlorophosphonazo detection.

Embodiment

Below in conjunction with accompanying drawing, the specific embodiment of the present invention is described, so that those skilled in the art understands the present invention better.Requiring particular attention is that, in the following description, when perhaps the detailed description of known function and design can desalinate main contents of the present invention, these are described in and will be left in the basket here.

Fig. 2 is the embodiment process flow diagram that the present invention is based on the harmonic detection method that non-base 2FFT converts.As shown in Figure 2, the harmonic detection method that the present invention is based on non-base 2FFT conversion comprises the following steps:

S201: sampled data is divided into groups:

Signal is sampled, sampling number N=2 is set ^m× N ', makes N ₁=2 ^m, N ₂=N ', is divided into N by FFT operational data ₂group, often organizes N ₁point data, n-th ₂in group n-th ₁point data is designated as f [n ₁, n ₂], the data sequence number n=N extracted in sampled data ₂n ₁+ n ₂, wherein n ₁span be 0≤n ₁≤ N ₁-1, n ₂span be 0≤n ₂≤ N ₂-1.

S202: often organize data and carry out base 2FFT conversion respectively:

To the N in every group ₁point data carries out base 2FFT conversion, by n-th ₂kth in group ₁individual transformation results is designated as F ' [k ₁, n ₂], k ₁span be 0≤k ₁≤ N ₁-1.

S203: calculate Fourier transform result:

Make k=k ₁+ N ₁k ₂, k ₂span be 0≤k ₂≤ N ₂-1, a kth Fourier transform result F [k of calculating sampling data ₁, k ₂] export as testing result, computing formula is:

$F[k_1,k_2]=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{W}_{N_2}^{\&prime;n_2{k}_2}{F}^{\&prime;}[k_1,n_2]---\left(7\right)$

Wherein, $W_{N_2}^{\&prime;n_2{k}_2}=W_{n_2}^{n_2{k}_2}{W}_N^{n_2{k}_1}.$

Below principle of the present invention and derivation are described.

In the present invention, when making N point FFT, N=2 ^m× N ', makes N ₁=2 ^m, N ₂=N ', so by the DFT of f [n]:

$F\left[k\right]\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}f\left[n\right]W_N^{\mathrm{nk}},0\&le;k\&le;N-1---\left(8\right)$

By mapping:

$\begin{array}{cc}n=N_2{n}_1+n_2& 0\&le;n_1\&le;N_1-1,0\&le;n_2\&le;N_2-1\\ k=k_1+N_1{k}_2& 0\&le;k_1\&le;N_1-\mathrm{1,0}\&le;k_2\&le;N_2-1\end{array}---\left(9\right)$

Can obtain: $W_N^{\mathrm{nk}}=W_N^{(N_2{n}_1+n_2)(k_1+N_1{k}_2)}=W_N^{(N_2{n}_1{k}_1+N_1{N}_2{n}_1{k}_2+n_2{k}_1+N_1{n}_2{k}_2)},$ And N=N ₁n ₂, $W_N^{N_1}=W_{N_2},$ $W_N^{N_2}=W_{N_1},$ can abbreviation be:

$W_N^{\mathrm{nk}}=W_{N_1}^{n_1{k}_1}{W}_N^{n_2{k}_1}{W}_{N_2}^{n_2{k}_2}---\left(10\right)$

Thus formula (8) is converted into:

$F[k_1,k_2]=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{W}_{N_2}^{n_2{k}_2}\left(W_N^{n_2{k}_1}\underset{n_1=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}f\right[n_1,n_2\left]W_{N_1}^{n_1{k}_1}\right)---\left(11\right)$

Wherein, 0≤k ₁≤ N ₁-1,0≤k ₂≤ N ₂-1.

Front two twiddle factors in (11) formula are merged, according to W _nreducibility can obtain

$W_{N_2}^{n_2{k}_2}{W}_N^{n_2{k}_1}=W_{N_1{N}_2}^{N_1{n}_2{k}_2}{W}_N^{n_2{k}_1}=W_N^{N_1{n}_2{k}_2+n_2{k}_1}=W_N^{n_2(N_1{k}_2+k_1)}---\left(12\right)$

(9) formula is substituted into (12) Shi Ke get

$W_{N_2}^{n_2{k}_2}\&times;W_N^{n_2{k}_1}=W_N^{n_2(N_1{k}_2+k_1)}=W_N^{n_{2}k}={W^{\&prime;}}_{N_2}^{n_2{k}_2}---\left(13\right)$

Thus formula (11) can turn to:

$F[k_1,k_2]=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{W}_{N_2}^{\&prime;n_2{k}_2}\left(\underset{n_1=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}f\right[n_1,n_2\left]W_{N_1}^{n_1{k}_1}\right)---\left(14\right)$

Visible, according to n ₂value f [n] is divided into N ₂group data f [n ₁, n ₂], often organize in data and have N ₁point, in formula (14) represent n-th ₂the kth obtained after group data transformation ₁individual transformation results, is designated as F ' [k ₁, n ₂].Due to N ₁=2 ^m, therefore N ₁the FFT conversion of point data can be obtained by base 2FFT algorithm.This pattern (14) just can turn to:

$F[k_1,k_2]=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{W}_{N_2}^{\&prime;n_2{k}_2}{F}^{\&prime;}[k_1,n_2]---\left(15\right)$

Formula (15) is similar to DFT and calculates, wherein not proper twiddle factor, but $W_{N_2}^{\&prime;n_2{k}_2}=W_{N_2}^{n_2{k}_2}{W}_N^{n_2{k}_1}.$

Embodiment

In order to technique effect of the present invention is described, implement invention has been emulation.Signal expression is y=0.8sin (10 π t)+0.6sin (80 π t)+sin (100 π t).According to IEC standard, the fundamental frequency of signal is 50Hz.In the middle of the present embodiment, harmonic frequency is 5Hz, 40Hz.Window width is 10 cycles.Sample to signal, each periodic sampling 1024 point, amounts to 10240 points.Fig. 3 is the result exemplary plot adopting the present invention to carry out m-Acetyl chlorophosphonazo detection.As shown in Figure 3, the m-Acetyl chlorophosphonazo testing result adopting the present invention to obtain is consistent with signal.

In actual applications, according to the requirement of signal characteristic and IEC standard, whole transformation results can not be needed, and only need to select a part of result, if adopt traditional FFT transform method, need first calculate whole result and then select, and in the present invention, owing to using formula k=k ₁+ N ₁k ₂carry out transformation results mapping, therefore can the sequence number of result as required select corresponding data to calculate, for the present embodiment, data are divided into 5 groups, first calculate the FFT transformation results of 5 groups of 2048 data, assuming that the transformation results sequence number of current needs is 511, so according to formula k=k ₁+ N ₁k ₂k can be tried to achieve ₁and k ₂, and then the data selecting 5 groups of FFT transformation results corresponding carry out calculating.Adopt in this way, computing time can be saved further.By MATLAB simulating, verifying, adopting the present invention all to calculate required time by 10240 is 0.062962 second.According to the requirement of IEC standard to survey frequency scope, only can calculate limited point, if for 2000, required time is 0.021317 second.Visible arithmetic speed of the present invention is very fast, can reach the measuring speed requirement that m-Acetyl chlorophosphonazo detects.

According to national act.std, signal fundamental frequency f is 50Hz, and sample window width D is 10 cycles, requires that frequency resolution is 5Hz.Remember that the sampling number in each cycle is X, so sample frequency f _s=Xf, total sampling number N=DX, so frequency resolution f ₀=Xf DX=f D.If adopt base 2FFT algorithm, so just need interpolation zero padding, extend that the sampling period reaches nearest 2 ^{m '}point, namely D ' is 16, and therefore frequency resolution f D ' is 3.125Hz, and be interpolation gained estimated value at 5Hz multiple frequency values, accuracy of detection can not ensure; And adopting the harmonic detection method that the present invention is based on non-base 2FFT algorithm, frequency resolution f D=5Hz, can meet m-Acetyl chlorophosphonazo accuracy of detection, and be up to state standards requirement.

The present invention is directed to points N=2 ^m× N '=N ₁× N ₂data, due to N ₁=2 ^mbe the whole power number of 2, the FFT built-in function of DSP can be called, calculate N ₂group N ₁point base 2FFT conversion, arithmetic speed is very fast, then calculates N ₁group N ₂the similar computing of some DFT, completes the Fourier transform of non-base 2.Compared with classic method, adopt the present invention to carry out m-Acetyl chlorophosphonazo detection, frequency resolution can be up to state standards, and accuracy of detection is higher, and arithmetic speed also can reach existing checkout equipment requirement.In addition, the present invention can calculate N=2 ^m× N ' counts, and range of application is more extensive.

Although be described the illustrative embodiment of the present invention above; so that those skilled in the art understand the present invention; but should be clear; the invention is not restricted to the scope of embodiment; to those skilled in the art; as long as various change to limit and in the spirit and scope of the present invention determined, these changes are apparent, and all innovation and creation utilizing the present invention to conceive are all at the row of protection in appended claim.

Claims (1)

1., based on the harmonic detection method that non-base 2FFT converts, it is characterized in that comprising the following steps:

S1: sample to signal, arranges sampling number N=2 ^m× N ', makes N ₁=2 ^m, N ₂=N ', is divided into N by sampled data ₂group, often organizes N ₁point data, n-th ₂in group n-th ₁point data is designated as f [n ₁, n ₂], the data sequence number n=N extracted in sampled data ₂n ₁+ n ₂, wherein n ₁span be 0≤n ₁≤ N ₁-1, n ₂span be 0≤n ₂≤ N ₂-1;

S2: to the N in every group ₁point data is carried out non-base 2 base 2FFT and is converted, by n-th ₂kth in group ₁individual transformation results is designated as F ' [k ₁, n ₂], k ₁span be 0≤k ₁≤ N ₁-1;

S3: make k=k ₁+ N ₁k ₂, k ₂span be 0≤k ₂≤ N ₂-1, a kth Fourier transform result F [k of calculating sampling data ₁, k ₂] and export as testing result, computing formula is:

$F[k_1,k_2]=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{W}_{N_2}^{\&prime;n_2{k}_2}{F}^{\&prime;}[k_1,n_2]$

Wherein, $W_{N_2}^{\&prime;n_2{k}_2}=W_{N_2}^{n_2{k}_2}{W}_N^{n_2{k}_1}.$