CN102072987B - Phase estimation method for short-interval sinusoidal signal and experimental device thereof - Google Patents

Abstract

The invention discloses a phase estimation method for a short-interval sinusoidal signal and an experimental device thereof. The CiR (cycles in record) of the tested sinusoidal signal is less than 1. The method comprises the following steps: sampling an inputted sinusoidal signal x(n)=Acos(2pif0+theta0) at an equal interval to obtain 2N-1 discrete sample sequences x(n); carrying out analytic transformation on the x(n) to obtain a complex sequence x1(n), and carrying out triangular window weighting and shift superposition to generate a sequence y(n) of which the length is N; computing the inner product of the sequence y(n) and a known sequence to obtain a complex value Q; and calculating the phase angle value of the complex value Q to obtain a phase measurement estimated value. The experimental device comprises a signal conditioning circuit, an A/D (analog-to-digital) conversion circuit, a digital signal processor (DSP), and an output drive and display unit. Compared with the prior art, the invention fills the gap of the short-interval sinusoidal wave phase estimation method within one cycle, has the advantages of high phase estimation precision, low computation quantity, high estimation efficiency and low resource consumption, and greatly saves hardware cost.

Description

Phase estimation method and the experimental provision thereof of short interval sinusoidal signal

Technical field

The invention belongs to digital signal processing technique field, be specifically related to carry out the digital measuring method of phase estimation between the area of observation coverage less than the sinusoidal signal of one-period.

Background technology

Sinusoidal signal x (n)=Acos (2 π f ₀+ θ ₀)+B (or claiming cosine signal, difference only is θ ₀Differ 90 °) be common, simple and typical signal.Obviously, the feature of sinusoidal signal is by frequency f ₀, amplitude A, phase theta ₀And four parameters of flip-flop B are portrayed fully.Wherein flip-flop B the most easily measures, and only needs the sampled point of all sinusoidal signals is averaging and can estimates, thereby is the simplified measurement method, does not consider the impact of flip-flop B in the elaboration of back of the present invention.

Usually sinusoidal signal is called again " sine wave ", and obviously, the implication of its " ripple " has reflected periodicity and the oscillatory of signal.Yet this also illustrates from the another side, if periodically and in the inadequate situation of concussion property signal is being observed, and will inevitably be to amplitude A, frequency f ₀And phase theta ₀The estimated accuracy of these parameters is brought impact.Wherein to phase theta ₀Have the greatest impact because for certain sampled point, phase place itself is a relative concept, namely only between the enough wide area of observation coverage in, enough large difference is arranged on amplitude between forward and backward sampled point, just can embody " phase place " feature.As shown in Figure 1, if estimate the phase place that sinusoidal wave P is ordered, Fig. 1 (a) is certainly easily many than Fig. 1 (b), because the sinusoidal signal in Fig. 1 has presented undulatory property within 2 cycles, and the not enough one-period of Fig. 1 (b) sinusoidal signal, undulatory property and oscillation characteristics are all not obvious.

For above-mentioned reasons, the Parameter Estimation Method of existing sine wave such as amplitude ratioing technigue ^[1,2,3], phase difference method ^[4,5,6], the energy barycenter method ^[7,8,9]Deng, substantially all proposing in enough " undulatory property " situations having, the money sky explicitly points out collection in document [10] sample will comprise more than 20 cycle at least; On the other hand, above various algorithms low frequency and observation time section when narrower a reason of degradation be to have ignored the real signal sine wave and be actually complex exponential ripple by two conjugation With

Form, corresponding is positive and negative frequency respectively, and when frequency was low, the frequency spectrum of these two complex exponential ripples suffered very closely, thus its spectrum separately leak into branch and influence each other to cause between spectrum and disturb, this all can reduce the parameter estimation of signal.

Fig. 2 has provided FFT spectrum (the Fast Fourier Transform of Fig. 1 (b) medium and low frequency signal, Fast Fourier Transform (FFT), get its exponent number N=512), can find out positive frequency complex exponential composition corresponding be k ∈ [0, the spectral line of left half axle N/2) bunch, and negative frequency complex exponential composition corresponding be k ∈ [N/2, the spectral line of right axle shaft N) bunch, because FFT has ring shift, thereby left and right two bunches of spectral lines are actually and are close to, and these two bunches of spectral lines exist between very serious spectrum and disturb obviously.

In recent years, domestic scholars Chen Kui inspires confidence in and has recognized that in document [11] there are the problems referred to above in the low-frequency sine parameter measurement.Point out that in document [12] the high and low of signal frequency is a relative concept, if analyzing samples contains abundant number cycle of fluctuation (CiR:cycles in record) can be considered high frequency, at this moment just can ignore the impact of negative frequency.Yet under a lot of situations, the CiR of collecting sample may not satisfy this condition (as Fig. 1 (b)), at this moment must consider the impact that the negative frequency composition brings.So document [12] has proposed the spectrum correction method of low-frequency component, but also only studied and to have counted CiR cycle of fluctuation low frequency signal parameter problem in 1 situation, document [12] has also only been made CiR from 1 simulation study that changes at 8 o'clock low-frequency parameter is estimated, count in CiR＜1 situation such cycle of fluctuation for Fig. 1 (b), how accurately to estimate low-frequency parameter, do not relate to.

In fact, in these parameters of sinusoidal signal, " phase place " weighed in time is first, rear relation in signals transmission, thereby the sine wave phase between the short area of observation coverage is measured have very high effect in military affairs, seismic prospecting, radio application.For example the signal frequency in underwater submarine communication is exactly often " extremely low frequency frequency " (the Extremely low frequency between 3～30Hz, ELF), and because observation time is restricted, thereby to measure the phase place of signal in the very short observation time of this section, its measuring accuracy just unavoidably can be subject to spectrum limit leakage shown in Figure 2, and then is difficult to estimate that enemy warship is from the distance of my warship.For another example, the seismic signal dominant frequency mainly changes between 20～50Hz, and reduce with the degree of depth, if thereby can estimate the phase place of signal when omen appears in earthquake, the phase place of far bottom layer signal overhead particularly, for even accomplish to forecast earthquake accurately several seconds of trying to be the first, all has very high meaning to retrieving national economy loss and casualties.

This patent will be introduced the technology such as quarter window weighting translation stack and Hilbert transform (Hilbert Transform), sine wave in proposing between the more perfect short area of observation coverage is measured estimation scheme, when counting CiR＜1 cycle of fluctuation when observation sample shown in Fig. 1 (b) to solve, sinusoidal wave two complex exponential signals that comprise exist in situation about disturbing between very serious spectrum, how to complete the estimation problem of high-precision phase parameter.

The list of references of this patent is as follows:

[1]G.Andria,M.Savino,and A.Trotta.“Windows and interpolation algorithms to improve electrical measurement accuracy,”IEEE Trans.Instrum.Meas.,vol.38,no.4,pp.856–863,Aug.1989.

[2] Liu Min, Wang Keying. based on the high precision FFT frequency analysis [J] of windowing double peak of spectral lines interpolation. electrical measurement and instrument, 2006,43 (3): 112-116.

[3] Qi is clear, Jia Xinle. the precision analysis of Frequency Estimation of Sinusoid Based on Interpolated FFT [J]. and electronic letters, vol, 2004,32 (4): 625-629.

[4] fourth health, Luojiang County is triumphant, Xie Ming. Time-Shifting Correcting Method of Phase Difference on Discrete Spectrum [J]. and authorship, 2002,23 (7): 729-735.

[5] Qi is clear, Jia Xinle. based on the sine wave freuqency of DFT and the high precision method of estimation [J] of first phase. and electronic letters, vol, 2001,29 (9): 1164-116.

[6] Yang Zhijian, the fourth health. the Frequency Estimation precision analysis [J] of time shift phase difference correction method under the white Gaussian noise background. vibration engineering journal, 2007,20 (3): 274-279.

[7] Zhu Xiaoyong, the fourth health. the comprehensive comparison [J] of discrete spectrum correction method. signal is processed, and 2001,17 (1): 91-97.

[8] fourth health, Xie Ming. the error analysis [J] of discrete spectrum three-point convolution amplitude rectification method. vibration engineering journal, 1996,9 (1): 92 – 98.

[9] fourth health, Jiang Liqi. the energy barycenter correction method [J] of discrete spectrum. vibration engineering journal, 2001,14 (3): 354 – 359.

[10] money sky, Zhao Rongxiang. based on frequency analysis between the interpolation FFT algorithm [J]. Proceedings of the CSEE, 2005,25 (21): 87-91.

[11] Chen Kuifu, Zhang Senwen. eliminate the Spectrum Correction [J] of negative frequency impact. physical strength 2004,26 (1): 25-28.

[12] Chen Kuifu, Wang Jianli, Zhang Senwen. the Spectrum Correction of low-frequency component [J]. vibration engineering journal, 2008,21 (1): 38-42.

Summary of the invention

Based on above-mentioned prior art, the present invention proposes a kind of method of estimation of phase parameter of short interval sinusoidal signal, realization is counted under the condition of CiR＜1 being observed sample cycle, to the high-precision method of estimation of the phase parameter of low-frequency sine signal, to the various sine wave parameter estimations technique the time carry out phase estimation.

The invention allows for a kind of experimental provision of method of estimation of the phase parameter of realizing above-mentioned short interval sinusoidal signal, be used for the optimal parameter settings such as corresponding sample frequency of the present invention.

The present invention proposes a kind of phase estimation method of short interval sinusoidal signal, CiR＜1 is counted by tested sinusoidal signal sample cycle, and the method comprises the following steps:

Step 1 is to the simulating signal x (n)=A cos (2 π f of input ₀+ θ ₀) carry out equal interval sampling, obtain 2N-1 discrete sample sequence x (n);

Step 2 carries out obtaining complex sequences x after analytic trnasformation to x (n) ₁(n), the concrete processing of analytic trnasformation comprises:

Will be as the discrete sample sequence x (n) of input signal minute two-way, wherein the first via as shown in Figure 4, need not done any processing, and the second the tunnel obtains sequence after need doing Hilbert transform

After Hilbert transform

Frequency spectrum

Be expressed as

$\hat{X}\left(\mathrm{j\ω}\right)=-j{X}_{+}\left(\mathrm{j\ω}\right)+j{X}_{-}\left(\mathrm{j\ω}\right)$

X in formula _-(j ω) is negative frequency spectrum, X ₊(j ω) is positive frequency spectrum;

Will Multiply by j and namely get analytic trnasformation output signal x again with after original signal x (n) stack ₁(n) frequency spectrum X ₁(j ω)

$X_1\left(\mathrm{j\ω}\right)=X\left(\mathrm{j\ω}\right)+j\hat{X}\left(\mathrm{j\ω}\right)=2X_{+}\left(\mathrm{j\ω}\right)$

Step 3 is the quarter window w of 2N-1 with length _c=[1/N, 2/N ..., (N-1)/N, 1, (N-1)/N ..., 2/N, 1/N] and to x ₁(n) generating length after weighting and translation stack is the sequences y (n) of N; Again with sequences y (n) and known array

Carry out obtaining complex values Q after inner product

$Q=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}y\left(n\right)e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HDA0000032153260000011.png,HDA0000032153260000012.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}n}$

Step 4, the angle values of getting complex values Q obtains the phase measurement estimated value

F is satisfied in arranging of described sample frequency _s(2N-1) f ₀/ CiR.

The invention allows for a kind of experimental provision of phase estimation method of short interval sinusoidal signal, input signal is tested sinusoidal signal, CiR＜1 is counted by its sample cycle, this experimental provision comprises that signal conditioning circuit, A/D converter, digital signal processor DSP and output drive and display unit, wherein, signal conditioning circuit is connected with A/D converter, then is connected with digital signal processor DSP, be connected at last that output drives and display unit as output:

Short interval sinusoidal input signal is simulated pre-service through signal conditioning circuit, so that the signal amplitude scope is adjusted, and removes outer interference noise; Pass through A/D converter, sampling obtains sample sequence x (n) again, with sample sequence x (n) and the V-CLK CP that obtains ₂As parallel digital signal supplied with digital signal processor DSP, call core algorithm through digital signal processor DSP, adjust according to actual needs sample rate f _s, complete the parameter estimation that receives signal and process; Drive by output at last and show and display module demonstrates the phase measurement estimated value;

Wherein, the signal parameter of digital signal processor DSP is estimated to process, and comprises the following steps:

According to concrete application requirements, guestimate signal frequency f ₀Count CiR with the signal observation cycle, and set according to specific needs the phase accuracy requirement;

The digital signal processor DSP internal processes is according to f ₀Count the CiR value with the signal observation cycle, set corresponding frequency spectrum exponent number N and fully collect the data sampling frequency f _s

The CPU primary controller of digital signal processor DSP enters internal RAM from IO port reads sampled data.

Realize follow-up " DC processing ": the mean value that calculates sample point obtains flip-flop B, obtains needing to eliminate sinusoidal signal x (n)=Acos (2 π f ₀+ θ ₀Flip-flop B in)+B;

Carry out phase measurement, obtain phase measurement.

Compared with prior art, the present invention can for sinusoidal wave observation cycle number less than 1 cycle situation, fill up the blank of sine wave phase measuring method; Have higher phase estimation precision, its phase estimation precision can be higher than the existing measuring accuracy of method in CiR＜1 situation, thereby precision is higher; Obtain the sequences y that length is N (n) after quarter window weighting translation stack, use again complex exponential sequence and y (n) to carry out inner product operation to replace the FFT computing, its computation complexity greatly reduces, calculated amount is little, estimated efficiency is high, also saved simultaneously the hardware resources such as multiplier, resource cost is few, has greatly saved hardware cost.Have stronger noise resisting ability, its noiseproof feature can improve by increasing the analysis of spectrum exponent number.

Description of drawings

Fig. 1 is sine waveform and the sampling point schematic diagram between the different areas of observation coverage;

Fig. 2 is the FFT spectrogram of low frequency cosine signal;

Fig. 3 is the schematic flow sheet of method of estimation of the phase parameter of short interval sinusoidal signal of the present invention;

Fig. 4 is that in the method for estimation flow process of phase parameter of short interval sinusoidal signal of the present invention, the signal resolution shift step is found the solution schematic diagram;

Fig. 5 is the forward and backward signal FFT spectral amplitude schematic diagram of analytic trnasformation in the method for estimation flow process of phase parameter of short interval sinusoidal signal of the present invention;

Fig. 6 is that the weighting translation stack of method of estimation flow process intermediate cam window and the inner product process of the phase parameter of short interval sinusoidal signal of the present invention found the solution schematic diagram;

Fig. 7 be in the method for estimation flow process of phase parameter of short interval sinusoidal signal of the present invention each the processing stage waveform schematic diagram relatively;

Fig. 8 is the FFT spectrum schematic diagram after quarter window translation stack;

Fig. 9 be under the noise of experimental result of the present invention each the processing stage waveform schematic diagram relatively;

Figure 10 is the hardware implementation figure for the experiment dress of the method for estimation of the phase parameter of the short interval sinusoidal signal of the present invention;

Figure 11 is the hardware DSP internal processes flow graph for the experimental provision of the method for estimation of the phase parameter of the short interval sinusoidal signal of the present invention.

Embodiment

As shown in Figure 3, at first to the simulating signal x (t)=Acos (2 π f of input ₀T+ θ ₀) carry out equal interval sampling and obtain 2N-1 discrete sample sequence x (n); X (n) is carried out obtaining discrete sample complex sequences x after analytic trnasformation ₁(n), to eliminate as far as possible the negative frequency sideband; Be that the quarter window of 2N-1 is to x with length ₁(n) generating length after weighting and translation stack is the sequences y (n) of N; Again with y (n) and known array

(this sequence is corresponding to the first harmonic of N rank FFT vector) carries out obtaining complex values Q(after inner product and namely calculates the summation of series

); The angle values of getting Q (is

) can obtain the phase measurement estimated value

Suppose that sample frequency is f _s, its sampling interval is 1/f _s, for satisfying the restriction (observation cycle is counted CiR value＜1) between the short area of observation coverage, should be less than CiR signal period T=1/f between the area of observation coverage of the 2N-1 that an is adopted sampling point ₀(f ₀Be the measured signal frequency), i.e. sample frequency f _s, analysis of spectrum exponent number N and signal frequency thereof the condition that should satisfy is

$(2N-1)\frac{1}{f_s}<\mathrm{CiR}\&CenterDot;\frac{1}{{\mathrm{<figure-callout\; id="0"\; label="f"\; filenames="HDA0000032153260000011.png,HDA0000032153260000012.png"\; state="\{\{state\}\}">f</figure-callout>}}_0}\&DoubleRightArrow;f_s>\frac{(2N-1){\mathrm{<figure-callout\; id="0"\; label="f"\; filenames="HDA0000032153260000011.png,HDA0000032153260000012.png"\; state="\{\{state\}\}">f</figure-callout>}}_0}{\mathrm{CiR}}---\left(1\right)$

Because the time of observing is short, and signal frequency f ₀Be low frequency, and existing analog to digital converter can reach very high switching rate, thereby formula (1) condition is very easily to satisfy in engineering.

Yet sample frequency is not more high better, because in the situation that fixing between the area of observation coverage, sample frequency is higher, the sampling point of being adopted is more, and the left and right FFT spectral line of low frequency cosine signal shown in Figure 2 bunch certainty be drawn close towards left and right both sides more.The front points out, and the spectrum at Fig. 2 left and right two ends in fact is close to, thereby this will inevitably make the phase mutual interference of the spectral line bunch of left and right two sidebands become more violent, is unfavorable on the contrary extracting phase information.This is also that existing each class methods can't be improved by improving sample frequency the basic reason of parameter estimation.

The below explains one by one to the concrete internal processes of each step of Fig. 3 and related theoretical principle.

(1) analytic trnasformation

Because signal comprises positive and negative frequency content, therefore in the present invention, the purpose of introducing analytic trnasformation is: the negative frequency composition of removing signal.

In the method for estimation flow process of the phase parameter of short interval sinusoidal signal of the present invention, analytic trnasformation specifically find the solution flow process as shown in Figure 4, namely as the discrete sample sequence x (n) of input signal minute two-way, wherein one the tunnel obtain sequence through after Hilbert transform

Will

Multiply by j and namely get analytic trnasformation output x again with after original signal x (n) stack ₁(n).

As everyone knows, the transport function of Hilbert transform in frequency domain is

The frequency spectrum X (j ω) that supposes x (n) can be expressed as suc as formula the positive frequency spectrum X shown in (3) and formula (4) ₊(j ω) and negative frequency spectrum X _-The form of (j ω) stack

X(jω)=X ₊(jω)+X _-(jω) (3)

And $X\left(\mathrm{j\&omega;}\right)=\left\{\begin{array}{cc}{X}_{+}\left(\mathrm{j\&omega;}\right)& \mathrm{\&omega;}>0\\ X_{-}\left(\mathrm{j\&omega;}\right)& \mathrm{\&omega;}<0\end{array}\right.---\left(4\right)$

According to formula (2), through after Hilbert transform, positive frequency spectrum X ₊The phase place of (j ω) turns 90 ° clockwise, and negative frequency spectrum X _-The phase place of (j ω) turns 90 ° counterclockwise, thereby the frequency spectrum after the Xi Er baud conversion can be expressed as

$\hat{X}\left(\mathrm{j\&omega;}\right)=-j{X}_{+}\left(\mathrm{j\&omega;}\right)+j{X}_{-}\left(\mathrm{j\&omega;}\right)---\left(5\right)$

Will

Multiply by j, then the signal spectrum X after superposeing with x (n) ₁(j ω) is

$X_1\left(\mathrm{j\&omega;}\right)=X\left(\mathrm{j\&omega;}\right)+j\hat{X}\left(\mathrm{j\&omega;}\right)=2X_{+}\left(\mathrm{j\&omega;}\right)---\left(6\right)$

Obviously, the spectrum X after analytic trnasformation ₁The positive spectral magnitude of (j ω) is strengthened, and becomes original 2 times, and negative spectrum component is by filtering.

But, limited by finite sample, when reality realizes with FFT, negative frequency spectrum X _-(j ω) can not be inhibited fully, and always some left behind.Fig. 5 has provided Low Frequency Sine Signals through the forward and backward FFT spectral amplitude of analytic trnasformation.Fig. 5 reflects intuitively, through after analytic trnasformation, and the negative frequency of right axle shaft spectrum X ₁(k) substantially be inhibited, thereby greatly reduced, the spectrum of a left side half frequency axis disturbed.Fig. 5 (b) has provided the partial enlarged drawing of left side band spectral line bunch, as seen composes leakage phenomenon still more serious.

(2) quarter window translation stack and inner product thereof are processed

As shown in Figure 6, for adopting the process of quarter window translation overlap-add procedure in the method for estimation flow process of the phase parameter of short interval sinusoidal signal of the present invention, take N=6 as example, input data and length be the quarter window [1 23 4565432 1]/6 of 2N-1 be weighted multiply each other after, the element in twos that with the space is again N=6 sampling unit carries out the translation stack and generated data y (0) ~ y (5), at last with y (0) ~ y (5) and sequence

Carry out obtaining plural Q after the inner product summation, the phase angle of getting the Q value namely obtains phase estimation

In Fig. 6, use quarter window w _c=[1/N, 2/N ..., (N-1)/N, 1, (N-1)/N ... 2/N, 1/N] purpose that the data of analytic trnasformation are weighted the translation stack is that the waveform that causes in order to eliminate sample waveform to be interrupted from beginning to end is discontinuous, this can play the effect that reduces to compose leakage.This is because based on the analysis of spectrum of Fourier transform, be actually to have given tacit consent to input signal is carried out analysis of spectrum after periodic extension, thereby the Wave data that is used for carrying out analysis of spectrum in fact should be regarded as to join end to end and connects.If waveform head and the tail data differences is too large, just be equivalent to introduce a very large waveform saltus step, large waveform saltus step occurs and mean that the analysis of spectrum result can introduce a lot of leakage compositions.FFT stave after amplifying in Fig. 5 (b) is bright, for simple signal, it should be ideally a frequency spectrum, utilize quarter window translation overlap-add procedure head and the tail waveform saltus step problem Fig. 6 owing to not resembling, thereby the spectral line after amplifying runs far deeper than one, but having leaked out many, spectrum leakage is more serious, must make the effect of phase estimation poorer.

Select different window function settings in Fig. 6, its final phase measurement can differ greatly.Reason is: (1) than other windows, the main lobe width of the fourier spectra that quarter window is corresponding is the narrowest, and this is conducive to further suppress the negative frequency sideband spectrum of positive frequency sideband is disturbed; (2) after noise process quarter window weighted sum translation overlap-add procedure, the degree that its noise energy weakens is larger than other windows, therefore noiseproof feature is better than selecting other windows.

Here with quarter window without the reason of other windows (as Hanning window, Hamming window) be only, (1) than other windows, the main lobe width of the fourier spectra that quarter window is corresponding is the narrowest, and this is conducive to further suppress the negative frequency sideband spectrum of positive frequency sideband is disturbed; (2) after noise process quarter window weighted sum translation overlap-add procedure, the degree that its noise energy weakens is larger than other windows, therefore noiseproof feature is better than selecting other windows.

As shown in Figure 7, provided waveform after original signal waveform, the waveform (comprising real part and imaginary part) after analytic trnasformation and quarter window weighting translation stack before processing with quarter window.As can be seen from Figure 7, the former sample waveform of Fig. 7 (a) has very large discontinuous from beginning to end, through after analytic trnasformation, real part shown in Fig. 7 (b) and imaginary part still have very large discontinuous from beginning to end, and after quarter window weighting translation stack pre-service shown in Figure 6, the real part of Fig. 7 (c) waveform and the head and the tail waveform non-continuous event of imaginary part are eliminated, it is that the analysis of spectrum result of N=512 is (from k=0 that Fig. 8 has provided corresponding to the analysis of spectrum exponent number, 1 ..., intercept out a section in k=0 ~ 8 in N-1).

FFT spectrum after quarter window translation stack as shown in Figure 8 is effectively left side band spectrum of intercepting part, can find out from Fig. 8 (a), after quarter window translation stack, its FFT amplitude spectrum amplifies spectrum than the part in Fig. 5 (b), its spectral line bunch radical that comprises drops to 1,2 from 6,7, thereby its spectrum leakiness greatly reduces.A benefit that reduces to compose leakiness is exactly can be directly to read out phase information from position corresponding to peak value spectral line: can further find from Fig. 8 (b), at the peak value spectrum k=1 place of Fig. 8 (a), its corresponding phase spectrum

Angle value substantially near 20 ° of phase values (measured value is 19.8105 °, only differs from 0.1895 °).

Because the computing formula of the N rank FFT of sequences y (n) is

$Y\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}y\left(n\right)e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HDA0000032153260000011.png,HDA0000032153260000012.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}\mathrm{nk}},k=\mathrm{0,1},...,N-1---\left(7\right)$

Less than a wave period situation, as shown in Fig. 8 (b), the peak value spectral line only drops on the position of k=1 for the sampling time width, thereby when asking phase value, Y (k) value of calculating other k ≠ 1 does not have much effects.Therefore k=1 in the formula of order (7) has

$Q=Y\left(1\right)=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}y\left(n\right)e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HDA0000032153260000011.png,HDA0000032153260000012.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}}---\left(8\right)$

The multiplication computation amount of formula (8) is only formula (7) 1/N, therefore consumes resources greatly reduces.

As shown in Figure 6, directly get the phase angle of Q value, namely obtain final phase estimation

The below is experimental result explanation of the present invention:

Table 1 is without the phase measurement (20 ° of true phase values) of this method in the situation of making an uproar

As can be seen from Table 1, the phase measurement accuracy of this method is higher, in the situation that the not enough half period of observation cycle number (CiR:cycles in record), its measuring error is 0.2 ° of left and right, and the measuring error of document [12] is in 0.3 °, and is to record under a plurality of CiR are far longer than 1 condition.

Above discussion only limits to without making an uproar situation.For existing in the noise situation, each the processing stage waveform as shown in Figure 9.Due to the randomness of noise, the Monte Carlo simulation of 500 times has been carried out in experiment, and counts the survey phase root-mean-square error of test of many times, and its statistics is as shown in table 2.

The phase measurement root-mean-square error of this method under table 2 noise situations

(20 ° of true phase values, signal to noise ratio snr=20dB)

As can be seen from Table 2, this measuring method has certain anti-noise ability, and its noiseproof feature can be composed exponent number N by increase and be improved.Contrast table 1, table 2 can find out, table 1 data of noise-free case show: when spectrum exponent number N increased to 512 from 128, frequency was that the phase measurement error of 1Hz signal only is reduced to-0.1895 ° from-0.2440 °; And have table 2 data of noise situations to show: when spectrum exponent number N increased to 512 from 128, frequency was that the phase measurement error of 1Hz signal is reduced to 2.6542 ° from 5.4279 °.This explanation: the increase of spectrum exponent number N (observe for maintaining in half signal period, corresponding sample frequency also needs to increase pro rata), for the improvement that the measuring accuracy under noise situations is arranged, when making an uproar than nothing, the improvement degree of measuring accuracy is much obvious.

Increasing analysis of spectrum N can alleviate the phase measurement error that noise causes and can be illustrated theoretically, this is because sinusoidal signal is the single-frequency narrow band signal, its energy can only be gathered on one or two FFT spectral line, less than 1 situation, only need to consider a spectral line at k=1 place for sampling period number of the present invention.And noise signal is broadband signal, and when making N rank FFT analysis of spectrum, its noise energy will be spread out on all N root spectral lines, therefore supposes that the power of noise is σ ², for a single spectral line, the energy of the noise effect that it is subject to is about σ ²/ N.And in fact the sine wave phase measurement is locked on a single spectral line at k=1 place, thereby obviously analysis of spectrum exponent number N is larger, σ ²/ N is less, surveys phase time also less by the error that noise causes.

The below gives simple declaration to implementing hardware of the present invention.Referring to Figure 10, for accurately estimating the phase parameter of short interval sinusoidal signal x (t), need simulate pre-service by signal conditioning circuit to input signal, so that the signal amplitude scope is carried out necessary adjustment, and remove outer interference noise etc.; Pass through again A/D(analog-to-digital conversion device) sampling obtains sample sequence x (n), and the form of inputting with Parallel Digital enters the DSP device, processes through the internal algorithm of DSP device, obtains the parameter estimation of signal; Drive by output at last and show and display module demonstrates the estimated value of frequency modulation rate and centre frequency.

Wherein, the DSP(Digital Signal Processor of Figure 10, digital signal processor) be core devices, in the signal parameter estimation procedure, complete following major function:

(1) call core algorithm, complete the parameter estimation that receives signal and process;

(2) adjust according to actual needs sample rate f _s(change the CP of Figure 10 ₂Clock frequency), make it under this sampling rate condition, estimate accurately signal parameter as far as possible;

Export to during (3) with the phase estimation fructufy and drive and display module.

Need point out, owing to having adopted digitized method of estimation, thereby determined the complexity of Figure 10 system, in real time degree be connected principal element with degree of stability be not that the periphery of DSP device in Figure 10 connects, but the kernel estimation algorithm that the DSP internal program memory is stored.

The internal processes flow process of DSP device as shown in figure 11.

The present invention implants " phase measurement of short interval sinusoidal signal " this kernel estimation algorithm that proposes in the DSP device, completes high precision, low complex degree, efficient phase estimation based on this.

Figure 11 flow process is divided into following several step:

(1) at first need according to concrete application requirements (the concrete measurement as seismic event and submarine requires), guestimate signal frequency f ₀Count CiR with the signal observation cycle, and set according to specific needs the phase accuracy requirement.This step is from engineering aspect proposition real needs, so that follow-up flow process is processed targetedly.

(2) then, the DSP internal processes is according to f ₀Count the CiR value with the signal observation cycle, set corresponding spectrum exponent number N and sample frequency f _s, purpose is to guarantee can fully collect data with current sample frequency in short observation time interval.

(3) then, the CPU primary controller enters internal RAM from IO port reads sampled data.

(4) follow-up " DC processing " is in order to eliminate sinusoidal signal x (n)=Acos (2 π f ₀+ θ ₀Flip-flop B impact in)+B.Otherwise the existence of flip-flop can reduce measuring accuracy.Flip-flop B is easy to measure, and the mean value that only need calculate sampling point can obtain.

(5) carrying out phase measurement by Fig. 3 processing procedure of the present invention is the most crucial part of DSP algorithm, after moving this algorithm, can obtain phase measurement.

(6) judge whether the inventive method satisfies engineering demand, if do not satisfy, program is returned, again according to f _s(2N-1) f ₀/ CiR requires to set sample frequency (as previously mentioned, can and correspondingly improve sample frequency and improve measuring accuracy by increase analysis of spectrum exponent number).

(7) until measurement result meets engine request, the output bus by DSP exports outside display drive device to, phase measurement is carried out number show.

Need point out, owing to having adopted the DSP realization, make the operation of whole parameter estimation become more flexible, the concrete condition of the various components that can comprise according to signal, change the inner parameter setting of algorithm by flexible in programming, as the exponent number N of analysis of spectrum, the spectral line number that participation is proofreaied and correct etc.

Claims (3)

1. the phase estimation method of a short interval sinusoidal signal, CiR＜1 is counted by tested sinusoidal signal sample cycle, and the method comprises the following steps:

Step 1 is to the simulating signal x (n)=A cos (2 π f of input ₀+ θ ₀) carry out equal interval sampling, obtain 2N-1 discrete sample sequence x (n);

Step 2 carries out obtaining complex sequences x after analytic trnasformation to x (n) ₁(n), the concrete processing of analytic trnasformation comprises:

Will be as the discrete sample sequence x (n) of input signal minute two-way, wherein the first via need not done any processing, and the second the tunnel obtains sequence after need doing Hilbert transform

After Hilbert transform

Frequency spectrum

Be expressed as

X in formula _-(j ω) is negative frequency spectrum, X ₊(J ω) is positive frequency spectrum;

Will Multiply by j and namely get analytic trnasformation output signal x again with after first via signal x (n) stack ₁(n) frequency spectrum X ₁(j ω)

Step 3 is the quarter window of 2N-1 with length

w _c=[1/N，2/N，…,(N-1)/N，1,(N-1)/N，…,2/N，1/N]

To x ₁(n) generating length after weighting and translation stack is the sequences y (n) of N; Again with sequences y (n) and known array

Carry out obtaining complex values Q after inner product

Step 4, the angle values of getting complex values Q obtains the phase measurement estimated value

2. the phase estimation method of short interval sinusoidal signal as claimed in claim 1, is characterized in that, f is satisfied in arranging of described sample frequency _s(2N-1) f ₀/ CiR.

3. experimental provision of realizing the phase estimation method of short interval sinusoidal signal as claimed in claim 1, input signal is tested sinusoidal signal, CiR＜1 is counted by its sample cycle, it is characterized in that, this experimental provision comprises that signal conditioning circuit, A/D converter, digital signal processor DSP and output drive and display unit, and wherein, signal conditioning circuit is connected with A/D converter, be connected with digital signal processor DSP again, be connected at last that output drives and display unit as output:

Short interval sinusoidal input signal is simulated pre-service through signal conditioning circuit, so that the signal amplitude scope is adjusted, and removes outer interference noise; Pass through A/D converter, sampling obtains sample sequence x (n) again, with sample sequence x (n) and the V-CLK CP that obtains ₂As parallel digital signal supplied with digital signal processor DSP, call core algorithm through digital signal processor DSP, adjust according to actual needs sample rate f _s, complete the parameter estimation that receives signal and process; Drive by output at last and show and display module demonstrates the phase measurement estimated value;

Wherein, the signal parameter of digital signal processor DSP is estimated to process, and comprises the following steps:

According to concrete application requirements, guestimate signal frequency f ₀Count CiR with the signal observation cycle, and set according to specific needs the phase accuracy requirement;

The digital signal processor DSP internal processes is according to f ₀Count the CiR value with the signal observation cycle, set corresponding frequency spectrum exponent number N and fully collect the data sampling frequency f _s

The CPU primary controller of digital signal processor DSP enters internal RAM from IO port reads sampled data;

Realize follow-up " DC processing ": the mean value that calculates sample point obtains flip-flop B, obtains needing to eliminate sinusoidal signal x (n)=Acos (2 π f ₀+ θ ₀Flip-flop B in)+B; Carry out phase measurement, obtain phase measurement.

Images