CN104931777B - A kind of signal frequency measuring method based on two DFT plural number spectral lines - Google Patents

Abstract

The invention relates to a signal frequency measurement method based on two DFT complex spectral lines, and belongs to the technical field of signal parameter measurement. The present invention is characterized in that its processing steps include: performing DFT transformation on the sampled signal after windowing processing, searching for the highest and second highest spectral lines by comparing the spectral line amplitude or the sum of the absolute values of the real and imaginary parts within the corresponding frequency range, based on The operation of the sum of the absolute values of the real and imaginary parts of two spectral lines is used to obtain the intermediate parameters, and then the frequency offset parameters are obtained by solving equations, or inverse function formulas, or approximating polynomials, and finally the signal frequency is measured. The invention avoids multiplication and square root calculation when searching for the highest and second highest spectral lines and calculating the intermediate parameter λ, simplifies the realization of signal frequency measurement algorithm, and shortens the calculation speed.

Description

A Signal Frequency Measurement Method Based on Two DFT Complex Lines

technical field

The invention relates to a signal frequency measurement method based on two DFT complex spectral lines, and belongs to the technical field of signal parameter measurement.

Background technique

At present, the method of analyzing frequency signals based on discrete Fourier transform DFT or its fast algorithm FFT has been widely used. However, DFT has a barrier effect, that is, the actual signal frequency may not fall on the discrete spectral line, so it is necessary to use an interpolation algorithm to estimate the actual signal frequency. In the article "Improved Algorithm for Harmonic Analysis of Power System Using FFT" published in "Journal of China Electrical Engineering Society" Volume 23, Issue 6 in 2003, it is proposed that after adding windowed Fourier transform to the input discrete signal, by selecting the highest amplitude and The second highest two spectral lines, the method of interpolating the frequency of the measured signal. If the discrete frequency numbers of the two spectral lines correspond to k ₁ and k ₂ =k ₁ +1 respectively, then the position k ₀ corresponding to the actual signal frequency satisfies k ₁ ≤k ₀ ≤k ₂ . Introducing an auxiliary parameter α=k ₀ -k ₁ -0.5, ignoring other signal interference, the value range of α is [-0.5,0.5]. Thus, the magnitudes |Y(k ₁ )| and |Y(k ₂ )| of the two discrete spectral lines satisfy:

When N is large, the above formula can be simplified as λ=g(α), and its inverse function is recorded as α=g ^-1 (λ). This paper further proposes to use polynomial approximation method to calculate

The shortcoming of the existing methods is that the highest and second highest spectral lines are searched based on the spectral line amplitude. However, obtaining the spectral line amplitude requires calculating the sum of the squares of the real part and the imaginary part, and then performing the square root, which requires a large amount of calculation.

Contents of the invention

The purpose of the present invention is to provide a method for measuring signal frequency based on two DFT complex spectral lines to solve the problem of large amount of calculation in the prior art.

To achieve the above object, the solution of the present invention includes:

A method for measuring signal frequency based on two DFT complex spectral lines, characterized in that the steps are as follows:

Step (1): Perform windowing processing on the sampling signal x(n) whose sampling rate is FS and sampling points are continuously intercepted to obtain a windowed signal y(n). The formula for windowing processing is:

y(n)=x(n)·w(n),

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

Step (2): Carry out discrete Fourier DFT transform to windowed signal y (n), obtain discrete frequency spectrum Y (k), wherein discrete frequency serial number k=0:(N-1);

Step (3): Within the range of discrete frequency numbers [ks,ke] corresponding to the set frequency range, search for the spectral line with the largest (|Re(Y(k))|+|Im(Y(k))|) As the highest spectral line, and compare the (|Re(Y(k))|+|Im(Y(k))|) values of the spectral lines on both sides of the spectral line, select the largest spectral line as the second highest spectral line, Record the discrete frequency numbers k1 and k2 of the two spectral lines, where 0<ks<ke<(N/2), k2=k1+1, Re(Z) is the real part of the complex number Z, and Im(Z) is the complex number Z the imaginary part of

Step (4): Calculate the intermediate parameter λ according to the two spectral lines Y(k1) and Y(k2) corresponding to k1 and k2:

Step (5): Solve the frequency offset parameter α in the following equation:

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

Step (6): Calculate the measured signal frequency fα according to the frequency offset parameter α, the calculation formula is:

f _α =(k ₁ +0.5+α)·Fs/N.

The process of searching the highest and second highest spectral lines in the step (3) is: within the range [ks, ke] of discrete frequency numbers corresponding to the set frequency range, search for the largest spectral line of |Y(k)| as the highest spectral line, and compare the spectral line |Y(k)| values on both sides of the largest spectral line, select the largest spectral line as the second highest spectral line, and record the discrete frequency numbers k1 and k2 of the two spectral lines, where 0<ks< ke<(N/2), k2=k1+1.

Described step (5) adopts the inverse function formula of intermediate parameter measured value λ and frequency offset parameter α relationship function λ=g(α) to calculate frequency offset parameter α, and this inverse function is:

α=g ⁻¹ (λ).

The step (5) adopts the approximate polynomial formula of the inverse function g ^-1 (λ) of the intermediate parameter measured value λ and the frequency offset parameter α relationship function λ=g(α) to calculate the frequency offset parameter α, and the approximate polynomial formula is :

Wherein, M is the highest degree of the approximation polynomial, and am(m=0:M) is the coefficient of the mth term λm of the polynomial.

The design principle of the signal frequency measurement method of the present invention is: assuming a single frequency signal x(t) whose frequency is f ₀ , amplitude is A, and initial phase is θ, after the analog-to-digital conversion of the sampling rate Fs, the following is obtained A discrete signal of the form:

If the time-domain form of the added window function is w(n), and the continuous frequency spectrum obtained by its discrete-time Fourier transform DTFT is W(ω), then the influence of the sidelobe of the frequency peak at the negative frequency point - f ₀ is ignored. The continuous spectrum function near the positive frequency point f ₀ can be expressed as:

The above formula performs discrete sampling, and the expression of the discrete Fourier transform DFT can be obtained as:

Wherein, the discrete frequency interval is Δf=F _S /N.

The cosine window function is the most commonly used type of window function in DFT. The unified time-domain form corresponding to the cosine window function is:

The discrete-time Fourier transform DTFT result of the cosine window w(n) is:

in:

In the main lobe of the DTFT spectrum curve of the signal, and when N is large, approximately:

when , the above formula takes an equal sign. According to the commonly used window function coefficients, within the main lobe -H<k<H, the adjacent two spectral lines W(ω) and The phase difference is approximately π; while the corresponding H<k<N/2 side lobe W(ω) and close to the same phase. Thus, when the phase difference between spectral lines Y(k ₁ ) and Y(k ₂ ) is 0 or π, there are:

Therefore, without calculating the spectral line amplitude, the search for the highest and second highest spectral lines and the calculation of the intermediate parameter λ can be realized directly through the sum of the absolute values of the real and imaginary parts of the spectral line. Furthermore, the actual signal frequency is further obtained. The method of the present invention is designed based on the above principle, omits some multiplication and square root operations, simplifies the realization of the signal frequency measurement algorithm, and shortens the calculation speed.

Description of drawings

FIG. 1 is a calculation process diagram of a signal frequency measurement method based on two DFT complex spectral lines according to an embodiment of the present invention.

Detailed ways

The present invention will be described in further detail below in conjunction with the accompanying drawings.

Two embodiments are provided below, and these two embodiments are applied to the measurement of frequency signals near 50 Hz.

Example 1

The specific steps are as follows:

Step (1): Sampling rate Fs=1500Hz, continuously intercepting signal x(n) of N=512 points, performing windowing processing to obtain windowing signal y(n), the formula of windowing processing is:

y(n)=x(n)·w(n),

Among them, w(n) selects the Hanning window function sequence of N=512 points, namely:

Step (2): Carry out discrete Fourier DFT transform to windowed signal y (n), obtain discrete frequency spectrum Y (k), wherein discrete frequency serial number k=0:(N-1);

Step (3): In the range of discrete frequency numbers [15, 23], that is, within the range of corresponding frequencies [43.945, 67.383] Hz, search for the largest spectral line of |Y(k)| as the highest spectral line, and compare the largest spectral line For the |Y(k)| values of the spectral lines on both sides, select the largest spectral line as the second highest spectral line, and record the discrete frequency numbers k ₁ and k ₂ =k ₁ +1 of the two spectral lines;

Step (4): Calculate the intermediate parameter λ according to the two spectral lines Y(k ₁ ) and Y(k ₂ ) corresponding to k ₁ and k ₂ :

Step (5): Solve the frequency offset parameter α in the following equation:

Therefore, the functional formula for calculating the frequency offset parameter α is:

α=1.5λ,

Step (6): Calculate the measured signal frequency f _α according to the frequency offset parameter α, the calculation formula is:

f _α =(k ₁ +0.5+α)·Fs/N.

Example 2

Step (1): Sampling rate Fs=1500Hz, continuously intercepting signal x(n) of N=512 points, performing windowing processing to obtain windowing signal y(n), the formula of windowing processing is:

y(n)=x(n)·w(n),

Among them, w(n) selects the Blackman (BlackMan) window function sequence of N=512 points, namely:

Step (2): Carry out discrete Fourier DFT transform to windowed signal y (n), obtain discrete frequency spectrum Y (k), wherein discrete frequency serial number k=0:(N-1);

Step (3): In the discrete frequency range [15,23], that is, within the range of the corresponding frequency [43.945,67.383]Hz, search for (|Re(Y(k))|+|Im(Y(k))| ) as the highest spectral line, and compare the (|Re(Y(k))|+|Im(Y(k))|) values of the spectral lines on both sides of the largest spectral line, and select the largest spectral line as Sub-high spectral line, record the discrete frequency numbers k ₁ and k ₂ =k ₁ +1 of the two spectral lines;

Step (4): Calculate the intermediate parameter λ according to the two spectral lines Y(k ₁ ) and Y(k ₂ ) corresponding to k ₁ and k ₂ :

Step (5): Using the approximate polynomial formula of the inverse function g ^-1 (λ) of the relationship function λ=g(α) between the measured value of the intermediate parameter λ and the frequency offset parameter α to calculate the frequency offset parameter α, using a maximum of 7 approximations A polynomial formula, of the form:

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

Step (6): Calculate the measured signal frequency f _α according to the frequency offset parameter α, the calculation formula is:

f _α =(k ₁ +0.5+α)·Fs/N.

According to the first and second implementation manners, the same set of simulation test data is respectively input to verify the calculation results of the two embodiments. The input signal x(n) is a signal whose fundamental frequency f ₁ is 50.1 Hz and contains 2nd to 9th harmonics, and the specific form is:

Among them, the amplitudes of the fundamental wave and each harmonic are: 1, 0.02, 0.1, 0.01, 0.05, 0.0, 0.02, 0.0, 0.01; the initial phases are -23.1°, 115.6°, 59.3°, 52.4°, 123.8°, 161.8°, -31.8°, 119.9°, -63.7°.

In the first embodiment, the calculation results of the nine spectral lines with discrete frequency numbers ranging from 15 to 23 are: 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. Regardless of the amplitude comparison or the sum of the absolute values of the real and imaginary parts, the highest and second highest spectral lines can be obtained corresponding to k ₁ =17 and k ₂ =18. Therefore, the calculation result of the intermediate parameter is λ=-0.26613833, the frequency offset parameter α=1.5λ=-0.39920750, the finally measured signal frequency is 50.099978 Hz, and the measurement error is -0.000044%.

In the second embodiment, the calculation results of the nine spectral lines with discrete frequency numbers ranging from 15 to 23 are: -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. Regardless of the amplitude comparison or the sum of the absolute values of the real and imaginary parts, the highest and second highest spectral lines can be obtained corresponding to k ₁ =17 and k ₂ =18. Therefore, the intermediate parameter calculation result is λ=-0.21160379. After being brought into the approximation polynomial, the frequency offset parameter α=-0.39909308 is obtained, and the finally measured signal frequency is 50.100313 Hz, and the measurement error is 0.00063%.

Specific embodiments have been given above, but the present invention is not limited to the described embodiments. The basic idea of the present invention lies in the above-mentioned basic scheme. For those of ordinary skill in the art, according to the teaching of the present invention, it does not need to spend creative labor to design various deformation models, formulas, and parameters. Changes, modifications, substitutions and variations to the implementations without departing from the principle and spirit of the present invention still fall within the protection scope of the present invention.

Claims (4)

1. A signal frequency measurement method based on two DFT complex spectral lines is characterized by comprising the following steps:

step (1): carrying out windowing processing on a sampling signal x (n) with a sampling rate of FS and a sampling point of continuous interception to obtain a windowed signal y (n), wherein the windowing processing formula is as follows:

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

wherein w (N) is a window function sequence of N points, N =0 (N-1);

step (2): performing Discrete Fourier Transform (DFT) on the windowed signal Y (N) to obtain a discrete spectrum Y (k), wherein the discrete frequency index k =0 (N-1);

and (3): searching a spectral line with the maximum (| Re (Y (k)) + | Im (Y (k)) |) in a discrete frequency number range [ ks, ke ] corresponding to a set frequency range as the highest spectral line, comparing (| Re (Y (k)) + | Im (Y (k)) |) values of spectral lines on two sides of the spectral line, selecting the largest spectral line as the next highest spectral line, and recording discrete frequency numbers k1 and k2 of the two spectral lines, wherein 0 & ltks ks & lt ke & lt (N/2), k2= k1+1, re (Z) is the real part of a complex number Z, and Im (Z) is the imaginary part of the complex number Z;

and (4): calculating an intermediate parameter lambda according to the two spectral lines Y (k 1) and Y (k 2) corresponding to k1 and k 2:

and (5): solving the frequency offset parameter α in the following equation:

wherein W (ω) is the result of a discrete-time fourier transform DTFT of the window function W (N) and the normalized angular frequency ω =2 π f/FS =2 π k/N;

and (6): and calculating the frequency f alpha of the measured signal according to the frequency offset parameter alpha, wherein the calculation formula is as follows:

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

2. a method for measuring signal frequency based on two DFT complex spectral lines as claimed in claim 1, wherein: the step (3) of searching the highest and second highest spectral lines comprises the following steps: in a discrete frequency number range [ ks, ke ] corresponding to a set frequency range, searching a spectral line with the maximum | Y (k) | as the highest spectral line, comparing the values of the spectral lines | Y (k) | on two sides of the maximum spectral line, selecting the maximum spectral line as the next highest spectral line, and recording discrete frequency numbers k1 and k2 of the two spectral lines, wherein 0 & lt ks & lt ke & lt (N/2), and k2= k1+1.

3. A method for measuring a signal frequency based on two DFT complex spectral lines as claimed in claim 1, wherein: the step (5) calculates the frequency offset parameter α by using an inverse function formula of the relation function λ = g (α) between the measured value λ of the intermediate parameter and the frequency offset parameter α, where the inverse function is:

α＝g ^-1 (λ)。

4. a method for measuring signal frequency based on two DFT complex spectral lines as claimed in claim 3, wherein: the step (5) adopts an inverse function g of a relation function lambda = g (alpha) of an intermediate parameter measured value lambda and a frequency deviation parameter alpha ^-1 (λ) calculating a frequency offset parameter α by an approximating polynomial equation:

where M is the highest degree of approximation to the polynomial and am (M = 0:M) is the coefficient of the mth term λ M of the polynomial.