CN104931777A - Signal frequency measurement method based on two DFT complex spectral lines - Google Patents

Abstract

The present invention relates to a signal frequency measurement method based on two DFT complex spectral lines, belonging to the field of signal parameter measurement technology. The method is characterized in that the processing procedure comprises that: a sampling signal is subjected to windowing processing and then is subjected to DFT transform, a highest spectral line and a second highest spectral line are searched through comparing spectral line amplitude values or the sum of real and imaginary part absolute values in a corresponding frequency search range, based on the operation of the sum of real and imaginary part absolute values of two spectral lines, offset parameters are obtained through a solution equation, or an inverse function formula or an approximation polynomial, and finally a signal frequency is measured. According to the method, multiplication and square root calculation are avoided in searching the highest spectral line and the second highest spectral line and calculating a middle parameter lambda, the signal frequency measurement algorithm realization is simplified, and the calculation speed is shortened.

Description

A kind of signal frequency measuring method based on two DFT plural number spectral lines

Technical field

The present invention relates to a kind of signal frequency measuring method based on two DFT plural number spectral lines, belong to signal parameter field of measuring technique.

Background technology

Current, the method based on discrete Fourier transformation DFT or its fast algorithm fft analysis frequency signal widely uses.But DFT has hurdle effect, namely actual signal frequency may not drop on discrete spectral line, needs thus to adopt interpolation algorithm to estimate the frequency of actual signal.Propose after input discrete signal enlarge section in " application FFT carries out the innovatory algorithm of harmonic analysis in power system " article that 2003 " Proceedings of the CSEE " 23 is delivered by volume 6 phase, the highest and secondary high two spectral lines by selection amplitude, the method for interpolation measuring-signal frequency.If two the discrete frequency sequence number of spectral line corresponding k respectively ₁and k ₂=k ₁+ 1, then the position k that actual signal frequency is corresponding ₀meet k ₁≤ k ₀≤ k ₂.Introduce an auxiliary parameter α=k ₀-k ₁-0.5, ignore other signal disturbing, then the numerical range of α is [-0.5,0.5].Thus, the amplitude of two discrete spectral lines | Y (k ₁) | with | Y (k ₂) | meet:

$\mathrm{\λ}=\frac{\left|Y\left(k_2\right)\right|-\left|Y\left(k_1\right)\right|}{\left|Y\left(k_2\right)\right|+\left|Y\left(k_1\right)\right|}=\frac{\left|W(2\mathrm{\π}\·(-\mathrm{\α}+0.5)/N)\right|-\left|W(2\mathrm{\π}\·(-\mathrm{\α}-0.5)/N)\right|}{\left|W(2\mathrm{\π}\·(-\mathrm{\α}+0.5)/N)\right|+\left|W(2\mathrm{\π}\·(-\mathrm{\α}-0.5)/N)\right|}.$

When n is large, above formula can be simplified shown as λ=g (α), and its inverse function is designated as α=g ^-1(λ).This article proposes to adopt Polynomial Approximation Method to calculate further

Existing methods deficiency is to carry out the highest and secondary high spectral line search based on spectral line amplitude.And obtain spectral line amplitude need to calculate real part and imaginary part quadratic sum, then carry out evolution, this calculated amount is large.

Summary of the invention

The object of this invention is to provide a kind of signal frequency measuring method based on two DFT plural number spectral lines, in order to solve the large problem of prior art calculated amount.

For achieving the above object, the solution of the present invention comprises:

Based on a signal frequency measuring method for two DFT plural number spectral lines, it is characterized in that, step is as follows:

Step (1): be FS by sampling rate, sampled point is the sampled signal x (n) intercepted continuously, carry out windowing process and obtain windowing signal y (n), windowing process formula is:

y(n)＝x(n)·w(n)，

Wherein w (n) is the window function sequence of N point, n=0:(N-1);

Step (2): discrete fourier DFT conversion is carried out to windowing signal y (n), obtains discrete spectrum Y (k), wherein discrete frequency sequence number k=0:(N-1);

Step (3): the discrete frequency serial number range [ks corresponding to setpoint frequency scope, ke] in, search (| Re (Y (k)) |+| Im (Y (k)) |) maximum spectral line is as the highest spectral line, and compare this spectral line both sides spectral line (| Re (Y (k)) |+| Im (Y (k)) |) value, select wherein maximum spectral line as secondary high spectral line, the discrete frequency sequence number k1 of record two spectral lines and k2, wherein 0<ks<ke< (N/2), k2=k1+1, the real part that Re (Z) is plural Z, the imaginary part that Im (Z) is plural Z,

Step (4): two the spectral line Ys (k1) corresponding according to k1 and k2 and Y (k2) calculate intermediate parameters λ:

$\mathrm{\&lambda;}=\frac{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)-\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)}{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)+\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)};$

Step (5): solve the frequency deviation parameter alpha in following equation:

$\mathrm{\&lambda;}=\frac{\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}+0.5)/N)\right|-\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}-0.5)/N)\right|}{\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}+0.5)/N)\right|+\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}-0.5)/N)\right|},$

Wherein, W (ω) is the result of the discrete time Fourier transform DTFT of window function w (n), and normalizing angular frequency=2 π f/FS=2 π k/N;

Step (6): calculate measured signal frequency f α according to frequency deviation parameter alpha, computing formula is:

f _α＝(k ₁+0.5+α)·Fs/N。

The process of the highest and secondary high spectral line of described step (3) search is the discrete frequency serial number range [ks corresponding to setpoint frequency scope, ke] in, search | Y (k) | maximum spectral line is as the highest spectral line, and more maximum spectral line both sides spectral line | Y (k) | value, select wherein maximum spectral line as secondary high spectral line, the discrete frequency sequence number k1 of record two spectral lines and k2, wherein 0<ks<ke< (N/2), k2=k1+1.

Described step (5) adopts the inverse function formulae discovery frequency deviation parameter alpha of intermediate parameters measured value λ and frequency deviation parameter alpha relation function λ=g (α), and this inverse function is:

α＝g ^-1(λ)。

Described step (5) adopts the inverse function g of intermediate parameters measured value λ and frequency deviation parameter alpha relation function λ=g (α) ^-1(λ) approximating polynomial formulae discovery frequency deviation parameter alpha, this approximating polynomial formula is:

$\mathrm{\&alpha;}\&ap;\underset{m=0}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\&CenterDot;{\mathrm{\&lambda;}}^m,$

Wherein, M is the most high reps of approximating polynomial, and am (m=0:M) is the coefficient of polynomial expression the m time item λ m.

The design concept of signal frequency measuring method of the present invention is: suppose that a frequency is f ₀, amplitude is A, initial phase is single-frequency signals x (t) of θ, have passed through the discrete signal obtaining following form after sampling rate is the analog to digital conversion of Fs:

$x\left(n\right)=A\&CenterDot;\sin (2\mathrm{\&pi;}\frac{f_0}{F_S}n+\mathrm{\&theta;})$

If the forms of time and space of institute's windowed function is w (n), the continuous frequency spectrum that its discrete time Fourier transform DTFT obtains is W (ω), then ignore negative frequency-f ₀the secondary lobe impact at frequency peak, place, at positive frequency f ₀neighbouring continuous frequency spectrum function can be expressed as:

$\stackrel{\&OverBar;}{X}\left(f\right)=\frac{A}{2j}{e}^{\mathrm{j\&theta;}}W\left(\frac{2\mathrm{\&pi;}(f-f_0)}{f_s}\right)$

Above formula carries out discrete sampling, and the expression formula that can obtain discrete Fourier transform (DFT) DFT is:

$Y\left(k\right)=\stackrel{\&OverBar;}{X}(k\&CenterDot;\mathrm{\&Delta;f})=\frac{A}{2j}{e}^{\mathrm{j\&theta;}}W\left(\frac{2\mathrm{\&pi;}(k\&CenterDot;\mathrm{\&Delta;f}-f_0)}{F_S}\right)$

Wherein, discrete frequency intervals is Δ f=F _s/ N.

Cosine window function is the class window function that DFT is the most conventional.The unified forms of time and space of corresponding cosine window function is:

$w\left(n\right)=\underset{h=0}{\overset{H}{\mathrm{\&Sigma;}}}{g}_h\&CenterDot;\cos (\frac{2\mathrm{\&pi;h}}{N}\&CenterDot;n),n=\mathrm{0,1},...,N-1$

The discrete time Fourier transform DTFT result of Cosine Window w (n) is:

$W\left(\mathrm{\&omega;}\right)=N\&CenterDot;e^{-j\frac{\mathrm{N\&omega;}}{2}}\&CenterDot;\sin \frac{\mathrm{N\&omega;}}{2}\&CenterDot;F\left(\mathrm{\&omega;}\right)$

Wherein:

$F\left(\mathrm{\&omega;}\right)=\underset{h=0}{\overset{H}{\mathrm{\&Sigma;}}}\frac{g_h}{2N}\&CenterDot;\frac{\sin \mathrm{\&omega;}}{\sin (\frac{\mathrm{\&omega;}}{2}-\frac{\mathrm{\&pi;h}}{N})\sin (\frac{\mathrm{\&omega;}}{2}+\frac{\mathrm{\&pi;h}}{N})}+j\&CenterDot;\underset{h=0}{\overset{H}{\mathrm{\&Sigma;}}}{g}_h$

In the main lobe of signal DTFT spectrum curve, and when n is large, approximate have:

$F\left(\frac{2\mathrm{\&pi;k}}{N}\right)\&ap;\frac{1}{\mathrm{\&pi;}}\underset{h=0}{\overset{H}{\mathrm{\&Sigma;}}}\frac{g_h\&CenterDot;k}{(k^2-h^2)}$

When time, above formula gets equal sign.According to conventional window function coefficient, in main lobe-H<k<H, its adjacent two spectral line W (ω) and phase be approximately π; And in the secondary lobe of corresponding H<k<N/2 W (ω) and close to same-phase.Thus, as spectral line Y (k ₁) and Y (k ₂) phase 0 or π time, have:

$\frac{\left|Y\left(k_1\right)\right|}{\left|Y\left(k_2\right)\right|}=\frac{\left|\mathrm{Re}\left(Y\left(k_1\right)\right)\right|+\left|\mathrm{Im}\left(Y\left(k_1\right)\right)\right|}{\left|\mathrm{Re}\left(Y\left(k_2\right)\right)\right|+\left|\mathrm{Im}\left(Y\left(k_2\right)\right)\right|}$

Thus, spectral line amplitude need not be calculated, directly pass through the absolute value sum of spectral line real part and imaginary part, the highest and secondary search of high spectral line and the calculating of intermediate parameters λ can be realized.And then, ask for actual signal frequency further.The inventive method, based on above-mentioned principle design, eliminates some multiplication and extracting operation, simplifies signal frequency Measurement Algorithm and realizes, shorten computing velocity.

Accompanying drawing explanation

Fig. 1 is the computation process figure of the signal frequency measuring method based on two DFT plural number spectral lines of the embodiment of the present invention.

Embodiment

Below in conjunction with accompanying drawing, the present invention will be further described in detail.

Below provide two embodiments, these two embodiments are applied to be measured frequency signal near 50Hz.

Embodiment 1

Its concrete steps are as follows:

Step (1): by sampling rate Fs=1500Hz, signal x (n) intercepting N=512 point continuously, carry out windowing process and obtain windowing signal y (n), windowing process formula is:

y(n)＝x(n)·w(n)，

Wherein w (n) selects the Hanning window function sequence of N=512 point, that is:

$w\left(n\right)=0.5-0.5\&CenterDot;\cos \left(\frac{2\mathrm{\&pi;n}}{N}\right),n=0:(N-1);$

Step (2): discrete fourier DFT conversion is carried out to windowing signal y (n), obtains discrete spectrum Y (k), wherein discrete frequency sequence number k=0:(N-1);

Step (3): at discrete frequency serial number range [15,23], i.e. respective frequencies [43.945, in the scope of 67.383] Hz, search | Y (k) | maximum spectral line is as the highest spectral line, and more maximum spectral line both sides spectral line | Y (k) | value, select wherein maximum spectral line as secondary high spectral line, the discrete frequency sequence number k of record two spectral lines ₁and k ₂=k ₁+ 1;

Step (4): according to k ₁and k ₂two corresponding spectral line Y (k ₁) and Y (k ₂) calculate intermediate parameters λ:

$\mathrm{\&lambda;}=\frac{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)-\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)}{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)+\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)};$

Step (5): solve the frequency deviation parameter alpha in following equation:

$\mathrm{\&lambda;}=\frac{\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}+0.5)/N)\right|-\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}-0.5)/N)\right|}{\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}+0.5)/N)\right|+\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}-0.5)/N)\right|}=\frac{2}{3}\mathrm{\&alpha;},$

Thus, the function formula calculating frequency deviation parameter alpha is:

α＝1.5λ，

Step (6): calculate measured signal frequency f according to frequency deviation parameter alpha _α, computing formula is:

f _α＝(k ₁+0.5+α)·Fs/N。

Embodiment 2

Step (1): by sampling rate Fs=1500Hz, signal x (n) intercepting N=512 point continuously, carry out windowing process and obtain windowing signal y (n), windowing process formula is:

y(n)＝x(n)·w(n)，

Wherein w (n) selects Blacknam (BlackMan) the window function sequence of N=512 point, that is:

$w\left(n\right)=0.42659-0.49656\&CenterDot;\cos \left(\frac{2\mathrm{\&pi;n}}{N}\right)+0.076849\&CenterDot;\cos \left(\frac{4\mathrm{\&pi;n}}{N}\right),n=0:(N-1);$

Step (2): discrete fourier DFT conversion is carried out to windowing signal y (n), obtains discrete spectrum Y (k), wherein discrete frequency sequence number k=0:(N-1);

Step (3): at discrete frequency serial number range [15,23], i.e. respective frequencies [43.945, in the scope of 67.383] Hz, search (| Re (Y (k)) |+| Im (Y (k)) |) maximum spectral line is as the highest spectral line, and more maximum spectral line both sides spectral line (| Re (Y (k)) |+| Im (Y (k)) |) value, select wherein maximum spectral line as secondary high spectral line, the discrete frequency sequence number k of record two spectral lines ₁and k ₂=k ₁+ 1;

Step (4): according to k ₁and k ₂two corresponding spectral line Y (k ₁) and Y (k ₂) calculate intermediate parameters λ:

$\mathrm{\&lambda;}=\frac{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)-\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)}{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)+\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)};$

Step (5): the inverse function g adopting intermediate parameters measured value λ and frequency deviation parameter alpha relation function λ=g (α) ^-1(λ) approximating polynomial formulae discovery frequency deviation parameter alpha, adopt the approximating polynomial formula being up to 7 times, form is as follows:

α≈1.87500108·λ+0.24171091·λ ³+0.10443757·λ ⁵+0.06718760·λ ⁷，

Step (6): calculate measured signal frequency f according to frequency deviation parameter alpha _α, computing formula is:

f _α＝(k ₁+0.5+α)·Fs/N。

According to first and second embodiment, input one group of identical emulation testing data respectively, to verify the result of calculation of two embodiments.This input signal x (n) is fundamental frequency f ₁for 50.1Hz, the signal comprising 2 to 9 subharmonic, concrete form is:

$x\left(n\right)=\underset{i=1}{\overset{9}{\mathrm{\&Sigma;}}}{A}_i\sin (2\mathrm{\&pi;}\frac{i\&CenterDot;f_1}{F_S}n+{\mathrm{\&theta;}}_i)$

Wherein, the amplitude of first-harmonic and each harmonic respectively: 1,0.02,0.1,0.01,0.05,0.0,0.02,0.0,0.01; Initial phase is-23.1 ° respectively, 115.6 °, 59.3 °, 52.4 °, 123.8 °, 161.8 ° ,-31.8 °, 119.9 ° ,-63.7 °.

In first embodiment, 9 spectral line result of calculations of discrete frequency serial number range 15 to 23 are successively: 0.152783905+j1.76229599,4.70210906+j54.2254920,-10.9858392-j126.688437,6.36734959+j73.4295187 ,-0.221526634-j2.55345712 ,-0.0512202109-j0.589197696,-0.0199885230-j0.228698046 ,-0.00996303987-j0.112669252.No matter amplitude com parison or adopt real part imaginary part absolute value sum, the corresponding k of the highest and secondary high spectral line all can be obtained ₁=17, k ₂=18.So intermediate parameters result of calculation is λ=-0.26613833, λ=-0.39920750, frequency deviation parameter alpha=1.5, finally measuring the signal frequency obtained is 50.099978Hz, measuring error-0.000044%.

In second embodiment, 9 spectral line result of calculations of discrete frequency serial number range 15 to 23 are successively :-0.63151373-j7.27542732,4.87561219+j56.23242377,-9.38422261-j108.21320679,6.10560557+j70.41558734 ,-1.11464837-j12.84931264,0.00734999+j0.0889938,-0.00022265+j0.00101658,0.00048733+j0.00853355.No matter amplitude com parison or adopt real part imaginary part absolute value sum, the corresponding k of the highest and secondary high spectral line all can be obtained ₁=17, k ₂=18.So intermediate parameters result of calculation is λ=-0.21160379.After bringing approximating polynomial into, obtain frequency deviation parameter alpha=-0.39909308, finally measuring the signal frequency obtained is 50.100313Hz, measuring error 0.00063%.

Be presented above concrete embodiment, but the present invention is not limited to described embodiment.Basic ideas of the present invention are above-mentioned basic scheme, and for those of ordinary skill in the art, according to instruction of the present invention, designing the model of various distortion, formula, parameter does not need to spend creative work.The change carried out embodiment without departing from the principles and spirit of the present invention, amendment, replacement and modification still fall within the scope of protection of the present invention.

Claims (4)

1., based on a signal frequency measuring method for two DFT plural number spectral lines, it is characterized in that, step is as follows:

Step (1): be FS by sampling rate, sampled point is the sampled signal x (n) intercepted continuously, carry out windowing process and obtain windowing signal y (n), windowing process formula is:

y(n)＝x(n)·w(n)，

Wherein w (n) is the window function sequence of N point, n=0:(N-1);

Step (2): discrete fourier DFT conversion is carried out to windowing signal y (n), obtains discrete spectrum Y (k), wherein discrete frequency sequence number k=0:(N-1);

Step (3): the discrete frequency serial number range [ks corresponding to setpoint frequency scope, ke] in, search (| Re (Y (k)) |+| Im (Y (k)) |) maximum spectral line is as the highest spectral line, and compare this spectral line both sides spectral line (| Re (Y (k)) |+| Im (Y (k)) |) value, select wherein maximum spectral line as secondary high spectral line, the discrete frequency sequence number k1 of record two spectral lines and k2, wherein 0<ks<ke< (N/2), k2=k1+1, the real part that Re (Z) is plural Z, the imaginary part that Im (Z) is plural Z,

Step (4): two the spectral line Ys (k1) corresponding according to k1 and k2 and Y (k2) calculate intermediate parameters λ:

$\mathrm{\&lambda;}=\frac{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)-\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)}{\left(\right|\mathrm{Re}\left(Y\left(k_2\right)\right)|+|\mathrm{Im}\left(Y\left(k_2\right)\right)\left|\right)+\left(\right|\mathrm{Re}\left(Y\left(k_1\right)\right)|+|\mathrm{Im}\left(Y\left(k_1\right)\right)\left|\right)};$

Step (5): solve the frequency deviation parameter alpha in following equation:

$\mathrm{\&lambda;}=\frac{\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}+0.5)/N)\right|-\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}-0.5)/N)\right|}{\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}+0.5)/N)\right|+\left|W(2\mathrm{\&pi;}\&CenterDot;(-\mathrm{\&alpha;}-0.5)/N)\right|},$

Wherein, W (ω) is the result of the discrete time Fourier transform DTFT of window function w (n), and normalizing angular frequency=2 π f/FS=2 π k/N;

Step (6): calculate measured signal frequency f α according to frequency deviation parameter alpha, computing formula is:

f _α＝(k ₁+0.5+α)·Fs/N。

2. a kind of signal frequency measuring method based on two DFT plural number spectral lines according to claim 1, it is characterized in that: the process of the highest and secondary high spectral line of described step (3) search is the discrete frequency serial number range [ks corresponding to setpoint frequency scope, ke] in, search | Y (k) | maximum spectral line is as the highest spectral line, and more maximum spectral line both sides spectral line | Y (k) | value, select wherein maximum spectral line as secondary high spectral line, the discrete frequency sequence number k1 of record two spectral lines and k2, wherein 0<ks<ke< (N/2), k2=k1+1.

3. a kind of signal frequency measuring method based on two DFT plural number spectral lines according to claim 1, it is characterized in that: described step (5) adopts the inverse function formulae discovery frequency deviation parameter alpha of intermediate parameters measured value λ and frequency deviation parameter alpha relation function λ=g (α), and this inverse function is:

α＝g ^-1(λ)。

4. a kind of signal frequency measuring method based on two DFT plural number spectral lines according to claim 3, is characterized in that: described step (5) adopts the inverse function g of intermediate parameters measured value λ and frequency deviation parameter alpha relation function λ=g (α) ^-1(λ) approximating polynomial formulae discovery frequency deviation parameter alpha, this approximating polynomial formula is:

$\mathrm{\&alpha;}\&ap;\underset{m=0}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m\&CenterDot;{\mathrm{\&lambda;}}^m,$

Wherein, M is the most high reps of approximating polynomial, and am (m=0:M) is the coefficient of polynomial expression the m time item λ m.