CN106405234A - Digital harmonic analysis method - Google Patents
Digital harmonic analysis method Download PDFInfo
- Publication number
- CN106405234A CN106405234A CN201610772684.6A CN201610772684A CN106405234A CN 106405234 A CN106405234 A CN 106405234A CN 201610772684 A CN201610772684 A CN 201610772684A CN 106405234 A CN106405234 A CN 106405234A
- Authority
- CN
- China
- Prior art keywords
- sequence
- point
- points
- algorithm
- wfta
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R23/00—Arrangements for measuring frequencies; Arrangements for analysing frequency spectra
- G01R23/16—Spectrum analysis; Fourier analysis
Landscapes
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
The invention discloses a digital harmonic analysis method comprising the steps that objects requiring harmonic analysis are sampled so that a real number sequence is obtained; and the obtained real number sequence is segmented, and harmonic analysis is performed on the segmented sequence by adopting a mixed radix algorithm, a PFA algorithm and a WFTA algorithm so as to complete digital harmonic analysis. According to the digital harmonic analysis method, the algorithm levels are clear and orderly performed, and the advantages of the four algorithms of real number decomposition, mixed radix, PFA and WFTA can be fully exerted; the computational burden is low, and translation sequence arrays and twiddle factors required in the computational process are computed in advance so that the computational burden is reduced; meanwhile, the measurement precision is high and reaches the national requirement of harmonic A-class computation; besides, the method is fast in computational speed, computational time consumption is far less than that of other independent algorithms, the engineering application effect is great, the objective of real-time and rapid harmonic analysis can be achieved, and powerful technical guarantee is provided for harmonic analysis of digital transformer substations.
Description
Technical field
Present invention relates particularly to a kind of digitized harmonic analysis method.
Background technology
Development with national economy technology and the raising of people's living standard, digitizing technique has been rooted in the hearts of the people, and
And in the middle of also having been widely used for the production of people and living.
Digital transformer substation is the trend of intelligent grid development, and with respect to traditional transformer station, data acquisition digitized is several
One of mark of Zi Hua transformer station technology application.According to the regulation of IEC61850-9-2, the sample frequency of combining unit is fixed
Several data are become, according to domestic power frequency 50HZ, the points that every cycle is commonly used are converted to 40,80,200, are wherein up to 200
Point.Because the points that every cycle is commonly used are more, then the waveform that digitized sampling obtains is more accurate, and therefore 200 points have become as
Every cycle commonly uses the trend of points.According to power quality standard《IEC61000-4-30》Requirement, Harmonics Calculation window becomes
10 cycles, corresponding sampled point has just become at 2000 points.Due to 2000 points and be unsatisfactory for 2 integral number power it is impossible to directly using FFT
Carry out frequency analyses.In such a case, it is possible to the method adopting has DFT algorithm, it is interpolated to the interpolation side of 2 integral number power
Method, based on the mixed base algorithm of factor composition principle, prime factor algorithm（Abbreviation PFA）, Winograd algorithm（Abbreviation WFTA calculates
Method）Etc., but the arithmetic speed of these methods single is not very efficient, does not reach the mesh of rapid computations in practical engineering application
, so existing harmonic analysis method also has some problems of practical application.
Content of the invention
It is an object of the invention to provide a kind of high digitized harmonic analysis method of calculating speed.
This digitized harmonic analysis method that the present invention provides, comprises the steps：
S1. the object needing to carry out frequency analyses is sampled, obtain the sequence of real numbers [x of n point1,x2,x3,……xn]；
S2. to the X point sequence of real numbers [x obtaining1,x2,x3,……xn] split, and the sequence after segmentation is respectively adopted
Mixed base algorithm, PFA algorithm and WFTA algorithm carry out frequency analyses, complete digitized frequency analyses.
Described digitized harmonic analysis method, after step S1, also comprises the steps before step S2：
Sequence of real numbers [the x to the n point obtaining for the S1.51,x2,x3,……xn], according to following formula by sequence [x1,x2,
x3,……xn] be converted to sequence of complex numbers Y [y1,y2,y3,……,ym]：
y1=x1+j*x2, y2=x3+j*x4, y3=x5+j*x6... ..., yi=x2i-1+j*x2i... ..., ym=xn-1+j*xn；Wherein m=n/
2；
Under above-mentioned steps, described in step S2 to the n point sequence of real numbers [x obtaining1,x2,x3,……xn] split, as
The sequence of complex numbers of m point is split.
Described splits to the m point sequence of complex numbers obtaining, and specifically includes following steps：
I. the characteristic according to WFTA algorithm and the value of m, the optimal calculating points A of selected WFTA algorithm；
II. the characteristic according to PFA algorithm and the value of m, the optimal computed points B of selected PFA algorithm；
III. the characteristic according to mixed base computational algorithm and the value of m, the optimal computed points C of selected mixed base computational algorithm；
IV. count A, B and C to selected calculating, is tested using following formula：
B=N1*A；
C=N2*B；
m=N3*C；
N2=N4*A；
In formula, N1, N2, N3 and N4 are natural number；
If V. check the sequence successfully m point sequence of complex numbers being partitioned into N3 group C points, and the sequence minute by C point
It is slit into the sequence for N2 B points, then the sequence sequences segmentation of B point being become N1 A points.
Described mixed base algorithm, PFA algorithm and WFTA algorithm be respectively adopted to the sequence after segmentation carry out frequency analyses,
Specifically include following steps：
. the sequence of m points is decomposed into the sequence of N3 group C points, carries out N3 time, the mixed base of each C point calculates；
. the sequence of C points is resolved into the sequence of N2 group B points, carry out N2 time, the PFA of each B point calculates；
. the sequence of B points is decomposed into the sequence of N1 group A points, carries out N1 time, the WFTA of each A point calculates；
. B point-number sequence after calculating step is changed, and the sequence of B points after conversion is decomposed into A
The sequence of N1 point, then carry out A time, the WFTA of each N1 point calculates；
. B point-number sequence after calculating step carries out translating sequence, and the PFA completing B point-number sequence calculates；Again by step
The sequence of the C points that ~ step obtains is multiplied by the twiddle factor of C point；
. the sequence of the postrotational C point obtaining step is decomposed into the sequence of B group N2 points, then carries out B time, often
Secondary N2 point PFA calculates；
. the sequence of N2 point is decomposed into the sequence of N4 group A points, carries out N4 time, the WFTA of each A point calculates；
. the sequence conversion of N2 point after step is calculated, then the sequence of N2 point after conversion is decomposed into A group N4
The sequence of individual points, then carry out A time, the WFTA of each N4 point calculates；
. the sequence of calculated for step N2 points is carried out N2 point PFA and translates sequence, complete the PFA meter of N2 point sequence
Calculate；
. the sequence of calculated for step ~ step C points is carried out translating sequence, completes the mixed base meter of C points
Calculate；
Xi. the sequence of calculated for step ~ step m points is rotated；
X. the sequence of m points obtaining step xi is decomposed into every group N3 sequence counted of C group, then carries out C time, often
The WFTA of secondary N3 points calculates；
X. the sequence of m points calculated to step x carries out translating sequence, completes the sequence of complex numbers analysis of m points；
X. the sequence of complex numbers of m points obtaining step x is converted to the sequence of real numbers of n point, and according to this real number sequence
Row carry out the calculating of harmonic amplitude and phase place, complete digitized frequency analyses.
Described mixed base algorithm, PFA algorithm and WFTA algorithm be respectively adopted to the sequence after segmentation carry out frequency analyses,
Also comprised the steps before step：
O. the segmentation result according to m points sequence of complex numbers, by fft algorithm, mixed base algorithm, PFA algorithm and WFTA algorithm
Need the constant used and constant calculated in advance and store, to improve the calculating speed of the inventive method.
The constant used and constant is needed to include m in described fft algorithm, mixed base algorithm, PFA algorithm and WFTA algorithm
The FFT of point combination translates ordinal number group, and C point mixed base translates ordinal number group, and B point PFA translates ordinal number group, N2 point PFA translates ordinal number group, A point WFTA
Trigonometric function constant, N4 point WFTA trigonometric function constant, N2 point WFTA trigonometric function constant, the twiddle factor of C point, m point
Twiddle factor and m point plural number turn the trigonometric function constant of real number.
Described n value is 2000；M value is 1000；N3 value is 5；C value is 200；N2 value is 10；B value is
20；N1 value is 4；A value is 5；N4 value is 2.
This digitized harmonic analysis method that the present invention provides, algorithm clear layer：Combinational algorithm clear layer, sequentially
Execution, gives full play to real number decomposition, mixed base, the advantage of tetra- algorithms of PFA, WFTA, maximizes favourable factors and minimizes unfavourable ones；And operand is little, meter
Need during calculation translate ordinal number group and twiddle factor calculates in advance, reduce operand；Certainty of measurement is high simultaneously, reaches GB
Harmonic wave A level calculates and requires；Additionally, this method fast operation, according to dominant frequency 300MHZ of industrial computer MCU, calculating speed is permissible
Reach 10ms, well below the tens and hundreds of ms of other independent algorithms, engineer applied effect is good, reaches harmonic wave and divides real-time
The purpose of analysis, provides strong technology to ensure to the frequency analyses of digital transformer substation.
Brief description
Fig. 1 is embodied as schematic flow sheet for the inventive method.
Specific embodiment
Be illustrated in figure 1 the inventive method is embodied as schematic flow sheet：This digitized harmonic wave that the present invention provides
Analysis method, samples to the object needing to carry out frequency analyses first, obtains the sequence of real numbers [x of n point1,x2,x3,……
xn]；Then to the X point sequence of real numbers [x obtaining1,x2,x3,……xn] split, and the sequence after segmentation is respectively adopted mixed
Close base algorithm, PFA algorithm and WFTA algorithm and carry out frequency analyses, complete digitized frequency analyses.
Below in conjunction with accompanying drawing 1, the inventive method is described in detail：With according to power quality standard《IEC61000-4-
30》Requirement, Harmonics Calculation window becomes 10 cycles, and each cycle is sampled 200 points, and corresponding sampled point has just become 2000
Point；The optimal calculating points A of selected WFTA algorithm is at 5 points, and the optimal computed points B of selected PFA algorithm is at 20 points, selectes mixing
The optimal computed points C of base computational algorithm is at 200 points；
（1）By the constant needing in fft algorithm, mixed base algorithm, PFA algorithm and WFTA algorithm to use and constant calculated in advance simultaneously
Storage, to improve the calculating speed of the inventive method：Need in fft algorithm, mixed base algorithm, PFA algorithm and WFTA algorithm to use
The constant arriving and constant, the FFT including 1000 points of combinations translates ordinal number group, and 200 points of mixed bases translate ordinal number group, and 20 points of PFA translate ordinal number
Group, 10 points of PFA translate ordinal number group, 5 points of WFTA trigonometric function constants, and 4 points of WFTA trigonometric function constants, 2 points of WFTA trigonometric functions are normal
Amount, the twiddle factor of 200 points, the twiddle factor of 1000 points and 1000 points of plural numbers turn the trigonometric function constant of real number,；
（2）By 2000 points of sequence of real numbers [x1,x2,x3,……x2000] become 1000 points of sequence of complex numbers [y1,y2,y3,……,
y1000]；
y1=x1+j*x2, y2=x3+j*x4, y3=x5+j*x6... ..., yi=x2i-1+j*x2i... ..., y1000=x1999+j*x2000；
（3）1000 points of sequence of complex numbers are resolved into 5*200, M=5, L=200, M are cycle-indexes, be circulated 5 times, every time 200
The mixed base of point calculates；
（4）200 point sequences are resolved into 10*20, P=10, Q=20, P are cycle-indexes, be circulated 10 times, 20 point sequence every time
PFA calculate；
（5）20 point sequences are resolved into 4*5, is circulated 4 times, the WFTA calculating of 5 point sequences every time；
（6）Will（5）After 20 point sequence ranks conversions after calculating, then 20 point sequences obtaining are resolved into 5*4 again, be circulated
5 times, the WFTA calculating of 4 point sequences every time；
（7）Will（6）What the sequence after calculating carried out 20 point sequences translates sequence, translates ordinal number group using initialized, completes 20 points of sequences
The PFA of row calculates, then will be from step（4）200 point sequences obtaining to this are multiplied by the twiddle factor that 200 point sequences initialize；
（8）With step（4）Corresponding, will（7）200 point sequences calculating resolve into 20*10, and Q=20, P=10, Q are circulation time
Number, is circulated 20 times, the PFA calculating of 10 point sequences every time；
（9）10 point sequences are resolved into 2*5, is circulated 2 times, 5 points of WFTA calculating every time；
（10）Will（9）After 10 point sequence ranks conversions after calculating, then 10 point sequences obtaining are resolved into 5*2, be circulated 5
Secondary, 2 points of WFTA calculates every time；
（11）Will（10）What the sequence after calculating carried out 10 points of PFA translates sequence, translates ordinal number group using initialized, completes 10 points of sequences
The PFA of row calculates；
（12）To from（4）Arrive（11）200 point sequences calculating carry out translating sequence, translate ordinal number group using initialized, complete 200
The mixed base of point calculates；
（13）Will be from（3）1000 point sequences obtaining to this are multiplied by the 1000 points of twiddle factors having initialized；
（14）Will（13）1000 point sequences calculating resolve into 200*5 again, and L=200, M=5, L are cycle-indexes, are circulated
200 times, 5 points of WFTA calculating every time；
（15）Right（14）1000 sequences calculating are translated ordinal number group using 1000 points having initialized and are carried out translating sequence, complete 1000
The sequence of complex numbers analysis of point；
（16）Will（15）The 1000 points of sequence of complex numbers obtaining carry out 2000 points of sequence of real numbers using the trigonometric function having initialized
Conversion；
（17）Real part according to 2000 points decompositing and imaginary part, carry out the calculating of each harmonic amplitude and phase place, complete 2000
The frequency analyses of point sequence.
Claims (7)
1. a kind of digitized harmonic analysis method, comprises the steps：
S1. the object needing to carry out frequency analyses is sampled, obtain the sequence of real numbers [x of n point1,x2,x3,……xn]；
S2. to the n point sequence of real numbers [x obtaining1,x2,x3,……xn] split, and the sequence after segmentation is respectively adopted mixed
Close base algorithm, PFA algorithm and WFTA algorithm and carry out frequency analyses, complete digitized frequency analyses.
2. digitized harmonic analysis method according to claim 1 is it is characterised in that after step S1, go back before step S2
Comprise the steps：
Sequence of real numbers [the x to the X point obtaining for the S1.51,x2,x3,……xn], according to following formula by sequence [x1,x2,x3,……
xn] be converted to sequence of complex numbers Y [y1,y2,y3,……,ym]：
y1=x1+j*x2, y2=x3+j*x4, y3=x5+j*x6... ..., yi=x2i-1+j*x2i... ..., ym=xn-1+j*xn；Wherein m=n/
2；
Under above-mentioned steps, described in step S2 to the n point sequence of real numbers [x obtaining1,x2,x3,……xn] split, as
The sequence of complex numbers of m point is split.
3. digitized harmonic analysis method according to claim 2 it is characterised in that described to the m point obtaining plural number sequence
Row are split, and specifically include following steps：
I. according to the characteristic of WFTA algorithm, the optimal calculating points A of selected WFTA algorithm；
II. according to the characteristic of PFA algorithm, the optimal computed points B of selected PFA algorithm；
III. according to the characteristic of mixed base computational algorithm, the optimal computed points C of selected mixed base computational algorithm；
IV. count A, B and C to selected calculating, is tested using following formula：
B=N1*A；
C=N2*B；
m=N3*C；
N2=N4*A；
In formula, N1, N2, N3 and N4 are natural number；
If V. check the sequence successfully m point sequence of complex numbers being partitioned into N3 group C points, and the sequence minute by C point
It is slit into the sequence for N2 B points, then the sequence sequences segmentation of B point being become N1 A points.
4. digitized harmonic analysis method according to claim 3 is it is characterised in that described divides to the sequence after segmentation
Frequency analyses are not carried out using mixed base algorithm, PFA algorithm and WFTA algorithm, specifically include following steps：
. the sequence of m points is decomposed into the sequence of N3 group C points, carries out N3 time, the mixed base of each C point calculates；
. the sequence of C points is resolved into the sequence of N2 group B points, carry out N2 time, the PFA of each B point calculates；
. the sequence of B points is decomposed into the sequence of N1 group A points, carries out N1 time, the WFTA of each A point calculates；
. B point-number sequence after calculating step is changed, and the sequence of B points after conversion is decomposed into A
The sequence of N1 point, then carry out A time, the WFTA of each N1 point calculates；
. B point-number sequence after calculating step carries out translating sequence, and the PFA completing B point-number sequence calculates；Again by step
The sequence of the C points that ~ step obtains is multiplied by the twiddle factor of C point；
. the sequence of the postrotational C point obtaining step is decomposed into the sequence of B group N2 points, then carries out B time, often
Secondary N2 point PFA calculates；
. the sequence of N2 point is decomposed into the sequence of N4 group A points, carries out N4 time, the WFTA of each A point calculates；
. the sequence conversion of N2 point after step is calculated, then the sequence of N2 point after conversion is decomposed into A group N4
The sequence of individual points, then carry out A time, the WFTA of each N4 point calculates；
. the sequence of calculated for step N2 points is carried out N2 point PFA and translates sequence, complete the PFA meter of N2 point sequence
Calculate；
. the sequence of calculated for step ~ step C points is carried out translating sequence, completes the mixed base meter of C points
Calculate；
Xi. the sequence of calculated for step ~ step m points is rotated；
X. the sequence of m points obtaining step xi is decomposed into every group N3 sequence counted of C group, then carries out C time, often
The WFTA of secondary N3 points calculates；
X. the sequence of m points calculated to step x carries out translating sequence, completes the sequence of complex numbers analysis of m points；
X. the sequence of complex numbers of m points obtaining step x is converted to the sequence of real numbers of n point, and according to this real number sequence
Row carry out the calculating of harmonic amplitude and phase place, complete digitized frequency analyses.
5. digitized harmonic analysis method according to claim 4 is it is characterised in that described divides to the sequence after segmentation
Frequency analyses are not carried out using mixed base algorithm, PFA algorithm and WFTA algorithm, also comprised the steps before step：
O. the segmentation result according to m points sequence of complex numbers, by fft algorithm, mixed base algorithm, PFA algorithm and WFTA algorithm
Need the constant used and constant calculated in advance and store, to improve the calculating speed of the inventive method.
6. digitized harmonic analysis method according to claim 5 is it is characterised in that fft algorithm, mixed base algorithm, PFA
Need the constant used in algorithm and WFTA algorithm and constant includes the FFT of m point combination and translates ordinal number group, C point mixed base translates ordinal number
Group, B point PFA translates ordinal number group, N2 point PFA translates ordinal number group, A point WFTA trigonometric function constant, N4 point WFTA trigonometric function constant, N2
The trigonometric function that point WFTA trigonometric function constant, the twiddle factor of C point, the twiddle factor of m point and m point plural number turn real number is normal
Amount.
7. digitized harmonic analysis method according to claim 6 is it is characterised in that described n value is 2000；M value
For 1000；N3 value is 5；C value is 200；N2 value is 10；B value is 20；N1 value is 4；A value is 5；N4 value is
2.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610772684.6A CN106405234A (en) | 2016-08-31 | 2016-08-31 | Digital harmonic analysis method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610772684.6A CN106405234A (en) | 2016-08-31 | 2016-08-31 | Digital harmonic analysis method |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN106405234A true CN106405234A (en) | 2017-02-15 |
Family
ID=58002296
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610772684.6A Pending CN106405234A (en) | 2016-08-31 | 2016-08-31 | Digital harmonic analysis method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106405234A (en) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102214159A (en) * | 2010-11-11 | 2011-10-12 | 福州大学 | Method for realizing 3780-point fast Fourier transform/inverse fast Fourier transform (FFT/IFFT) and processor thereof |
-
2016
- 2016-08-31 CN CN201610772684.6A patent/CN106405234A/en active Pending
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102214159A (en) * | 2010-11-11 | 2011-10-12 | 福州大学 | Method for realizing 3780-point fast Fourier transform/inverse fast Fourier transform (FFT/IFFT) and processor thereof |
Non-Patent Citations (3)
| Title |
|---|
| PKOLBA等: "A Prime Faetor FFT algorithm using high一speed convolution", 《IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING》 * |
| 余飞: "DTMB系统中3780点FFT处理器的算法设计及FPGA实现", 《中国优秀硕士学位论文全文数据库 信息科技辑》 * |
| 段义隆: "适用于数字化变电站的电能质量检测算法的研究", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑（月刊）》 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102288807B (en) | Method for measuring electric network voltage flicker | |
| CN104897960B (en) | Harmonic wave rapid analysis method and system based on the spectral line interpolation FFT of adding window four | |
| CN105004939B (en) | A kind of complex electric energy quality disturbance signal quantitative analysis method | |
| CN103926462B (en) | Rapid harmonic wave analyzing method and device of power system | |
| CN109633262A (en) | Three phase harmonic electric energy gauging method, device based on composite window multiline FFT | |
| CN101900761B (en) | High-accuracy non-integer-period sampled harmonic analysis and measurement method | |
| CN109946512A (en) | A kind of dynamic power analysis method for improving frequency domain interpolation | |
| CN102520245A (en) | Micro-grid harmonic and inter-harmonic analysis method based on cubic spline interpolation waveform reconstruction | |
| CN104146709A (en) | Quick acquiring method for multi-frequency-point bioelectrical impedance | |
| CN103983849B (en) | A kind of Electric Power Harmonic Analysis method of real-time high-precision | |
| CN103529294A (en) | HHT (Hilbert-Huang Transform)-based harmonic detection system and method for grid-connected inverter of photovoltaic system | |
| CN106250904A (en) | Based on Power Disturbance analyser and the sorting technique of improving S-transformation | |
| CN109655665A (en) | All phase Fourier's harmonic analysis method based on Blackman window | |
| CN109446643B (en) | Method for establishing household appliance load harmonic model based on measured data | |
| CN110068729A (en) | A kind of signal phasor calculating method | |
| CN109581045A (en) | A kind of m-Acetyl chlorophosphonazo power measurement method meeting IEC standard frame | |
| Stanisavljević et al. | Reduced FFT algorithm for network voltage disturbances detection | |
| CN103728616A (en) | Field programmable gate array (FPGA) based inverse synthetic aperture radar (ISAP) imaging parallel envelope alignment method | |
| CN105548711B (en) | A kind of multifrequency information filter recursive demodulation method | |
| CN105372492A (en) | Signal frequency measurement method based on three DFT complex spectral lines | |
| CN104931777A (en) | Signal frequency measurement method based on two DFT complex spectral lines | |
| CN106405234A (en) | Digital harmonic analysis method | |
| CN102170318B (en) | Spectral analysis algorithm used for receiver performance test | |
| CN108982954B (en) | Method and system for calculating phase voltage amplitude and phase suitable for feeder line terminal | |
| CN107085133A (en) | Method and device for calculating single-phase virtual value |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20170215 |
|
| RJ01 | Rejection of invention patent application after publication |