CN102072987A - 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

The 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.The wherein the easiest measurement of flip-flop B only needs the sampled point of all sinusoidal signals is asked on average and can be estimated, thereby is the simplified measurement method, does not consider the influence of flip-flop B in the elaboration of back of the present invention.

Usually sinusoidal signal is called " sine wave " again, and obviously, the implication of its " ripple " has reflected the periodicity and the oscillatory of signal.Yet this also illustrates from the another side, if periodically and under 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 influence.Wherein to phase theta ₀Have the greatest impact because for certain sampled point, phase place itself is a relative notion, just only between the enough wide area of observation coverage in, enough big difference is arranged on the amplitude between forward and backward sampled point, just can embody " phase place " feature.As shown in Figure 1, if will estimate the phase place that sinusoidal wave P is ordered, Fig. 1 (a) is easily more than Fig. 1 (b) certainly, because the sinusoidal signal among Fig. 1 has presented undulatory property in 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 under enough " undulatory property " situations having, the money sky then spells out collection in document [10] sample will comprise more than 20 cycle at least; On the other hand, a reason of various algorithms degradation when low frequency and observation time section are narrower is to have ignored the real signal sine wave and be actually complex exponential ripple by two conjugation more than

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 correspondence be k ∈ [0, the spectral line of left half axle N/2) bunch, and that negative frequency complex exponential composition correspondence is 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 this 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 notion,, at this moment just can ignore the influence of negative frequency if analyzing samples contains abundant number cycle of fluctuation (CiR:cycles in record) then can be considered high frequency.Yet under a lot of situations, the CiR that gathers sample may not satisfy this condition (as Fig. 1 (b)), at this moment must consider the influence that the negative frequency composition brings.So document [12] has proposed the spectrum correction method of low-frequency component, but also only studied the low frequency signal parameter problem of counting cycle of fluctuation under CiR＞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 under 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 this Several Parameters of sinusoidal signal, " phase place " weighed in time is first, back relation in the 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 the underwater submarine communication often is exactly " extremely low frequency frequency " (the Extremely low frequency between 3～30Hz, ELF), and because observation time is restricted, thereby in the very short observation time of this section, measure the phase place of signal, its measuring accuracy just unavoidably can be subjected to spectrum limit leakage shown in Figure 2, and then is difficult to estimate the distance of enemy warship from my warship.For another example, the seismic signal dominant frequency mainly changes between 20～50Hz, and reduce with the degree of depth, thereby when omen appears in earthquake if can estimate the phase place of signal, the phase place of overhead far bottom layer signal 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 quarter window weighting translation stack and Hilbert transform technology such as (Hilbert Transform), sine wave in proposing between the more perfect short area of observation coverage is measured estimation scheme, when counting CiR＜1 the cycle of fluctuation shown in Fig. 1 (b) when observation sample to solve, sinusoidal wave two complex exponential signals that comprised exist under the situation about disturbing between very serious spectrum, how to finish 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?improveelectrical?measurement?accuracy，”IEEE?Trans.Instrum.Meas.，vol.38，no.4，pp.856-863，Aug.1989.

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

[3] Qi is clear, Jia Xinle. and interpolation FFT is estimated precision analysis [J] electronic letters, vol of sinusoidal signal frequency, 2004,32 (4): 625-629.

[4] fourth health, Luojiang County is triumphant, Xie Ming. discrete spectrum time shift phase difference correction method [J] applied mathematics and mechanics, 2002,23 (7): 729-735.

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

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

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

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

[9] fourth health, Jiang Liqi. energy barycenter correction method [J] the vibration engineering journal of discrete spectrum, 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. the frequency spectrum of eliminating the negative frequency influence is proofreaied and correct [J] physical strength 2004,26 (1): 25-28.

[12] Chen Kuifu, Wang Jianli, Zhang Senwen. the frequency spectrum of low-frequency component is proofreaied and correct [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 lacking the phase parameter of interval sinusoidal signal, be implemented in and be observed sample cycle and count under the condition of CiR＜1, 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 the interval sinusoidal signal of above-mentioned weak point, be used for optimal parameter settings such as corresponding sample frequency of the present invention.

The present invention proposes a kind of phase estimation method of lacking interval sinusoidal signal, CiR＜1 is counted by tested sinusoidal signal sample cycle, and this method may further comprise the steps:

Step 1 is to simulating signal x (n)=Acos (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 the analytic trnasformation to x (n) ₁(n), the concrete processing of analytic trnasformation comprises:

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

Behind the Xi Er baud conversion 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 the formula _-(j ω) is negative frequency spectrum, X ₊(j ω) is positive frequency spectrum;

Will Multiply by j and promptly 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)/and N, 1, (N-1)/N ..., 2/N, 1/N] to x ₁(n) generation length in weighting and translation stack back is the sequences y (n) of N; Again with sequences y (n) and known array

Carry out obtaining complex values Q after the 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,HDA0000032153260000031.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 being provided with of described sample frequency _s＞(2N-1) f ₀/ CiR.

The invention allows for a kind of experimental provision of realizing the phase estimation method of the interval sinusoidal signal of weak point as claimed in claim 1, input signal is that CiR＜1 is counted by tested sinusoidal signal sample cycle, it is characterized in that, this experimental provision comprises that signal conditioning circuit, A/D change-over circuit, digital signal processor DSP and output drive and display unit, wherein, signal conditioning circuit is connected with the A/D change-over circuit, 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 again, sampling obtains sample sequence x (n), 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 sample rate f according to actual needs _s, finish the parameter estimation of received signal and handle; Demonstrate the estimated value of frequency modulation rate and centre frequency at last by output driving demonstration and display module thereof.

Wherein, the signal parameter of digital signal processor DSP is estimated to handle, and may further comprise the steps:

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

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 among the)+B;

Carry out phase measurement, obtain phase measurement.

Compared with prior art, the present invention can fill up the blank of sine wave phase measuring method at sinusoidal wave observation cycle number less than 1 cycle situation; Have higher phase estimation precision, its phase estimation precision can be higher than the existing measuring accuracy of method under CiR＜1 situation, thereby precision is higher; After quarter window weighting translation stack, obtain the sequences y that length is N (n), use complex exponential sequence and y (n) to carry out inner product operation again to replace the FFT computing, its computation complexity greatly reduces, calculated amount is little, the estimated efficiency height, also saved simultaneously hardware resources such as multiplier, resource cost is few, has saved hardware cost greatly.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 synoptic 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 the interval sinusoidal signal of weak point of the present invention;

Fig. 4 finds the solution synoptic diagram for signal resolution shift step in the method for estimation flow process of the phase parameter of the interval sinusoidal signal of weak point of the present invention;

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

Fig. 6 finds the solution synoptic diagram for the weighting translation stack of method of estimation flow process intermediate cam window and the inner product process of the phase parameter of the interval sinusoidal signal of weak point of the present invention;

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

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

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

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

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

Embodiment

As shown in Figure 3, at first to the input simulating signal x (t)=Acos (2 π f ₀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 the analytic trnasformation ₁(n), to eliminate the negative frequency sideband as far as possible; With length is that the quarter window of 2N-1 is to x ₁(n) generation length in weighting and translation stack back is the sequences y (n) of N; Again with y (n) and known array

(this sequence is corresponding to the first harmonic vector of N rank FFT) carries out obtaining complex values Q after the inner product and (promptly 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, then 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 between the area of observation coverage of the 2N-1 that an is adopted sampling point less than CiR signal period T=1/f ₀(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 being observed 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 to satisfy very easily in engineering.

Yet sample frequency is not high more good more, because under situation fixing between the area of observation coverage, sample frequency is high more, the sampling point of being adopted is many 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 that the spectrum at the left and right two ends of Fig. 2 in fact is close to, thereby this will inevitably make the phase mutual interference of spectral line bunch of left and right two sidebands become more violent, is unfavorable for extracting phase information on the contrary.This also is that existing each class methods can't be improved the basic reason of parameter estimation by improving sample frequency.

Below the concrete internal processes of each step of Fig. 3 and related theoretical principle are done explanation one by one.

(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 the interval sinusoidal signal of weak point of the present invention, analytic trnasformation specifically find the solution flow process as shown in Figure 4, promptly divide two-way as the discrete sample sequence x (n) of input signal, wherein one the tunnel obtain sequence through after the Hilbert transform

Will

Multiply by j and promptly 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 the 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 the Hilbert transform, positive frequency spectrum X ₊The phase place of (j ω) is changeed 90 ° clockwise, and negative frequency spectrum X _-The phase place of (j ω) is changeed 90 ° counterclockwise, thereby the frequency spectrum behind 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, the signal spectrum X after superposeing with x (n) again ₁(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 the 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 limited 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 the low frequency sinusoidal signal through the forward and backward FFT spectral amplitude of analytic trnasformation.Fig. 5 reflects that intuitively through after the analytic trnasformation, the negative frequency of right axle shaft is composed X ₁(k) be inhibited substantially, the spectrum of a left side half frequency axis disturbed thereby significantly reduced.Fig. 5 (b) has provided the partial enlarged drawing of left side band spectrum line bunch, and it is still relatively more serious as seen to compose leakage phenomenon.

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

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 the interval sinusoidal signal of weak point of the present invention, with N=6 is example, input data and length be the quarter window [12345654321]/6 of 2N-1 be weighted multiply each other after, be that the element in twos of N=6 sampling unit carries out the translation stack and generates data y (0)～y (5) again with the space, 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 promptly obtains phase estimation

Among 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 behind the 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 big, just be equivalent to introduce a very big waveform saltus step, big waveform saltus step occurs and mean that then the analysis of spectrum result can introduce a lot of leakage compositions.FFT stave after amplifying among Fig. 5 (b) is bright, for simple signal, it ideally should be a frequency spectrum, owing to do not resemble and utilize quarter window translation overlap-add procedure head and the tail waveform saltus step problem Fig. 6, thereby the spectral line after amplifying runs far deeper than one, but having leaked out many, spectrum leakage is serious more, must make that then the effect of phase estimation is poor more.

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

Here only with quarter window without the reason of other windows (as Hanning window, Hamming window) be, (1) than other windows, the main lobe width of the fourier spectra of quarter window correspondence is the narrowest, and this helps further suppressing the negative frequency sideband spectrum of positive frequency sideband is disturbed; (2) after the noise process quarter window weighted sum translation overlap-add procedure, the degree that its noise energy weakens is bigger than other windows, so 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 the analytic trnasformation and the quarter window weighting translation stack before handling with quarter window.As can be seen from Figure 7, the former sample waveform of Fig. 7 (a) has very big discontinuous from beginning to end, through after the analytic trnasformation, real part shown in Fig. 7 (b) and imaginary part still have very big 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 among the N-1).

FFT spectrum after the quarter window translation stack as shown in Figure 8 is an effectively left side band spectrum of intercepting part, from Fig. 8 (a) as can be seen, after quarter window translation stack, its FFT amplitude spectrum amplifies spectrum than the part among Fig. 5 (b), its spectral line bunch radical that is comprised 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 the position of peak value spectral line correspondence: can further find from Fig. 8 (b), at the peak value spectrum k=1 place of Fig. 8 (a), its pairing phase spectrum

(1) angle value is 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,HDA0000032153260000031.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}\mathrm{nk}},$ k＝0，1，...，N-1(7)

Less than a wave period situation, 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.So make k=1 in the formula (7), then have

$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,HDA0000032153260000031.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;}}{N}n}---\left(8\right)$ The multiplication computation amount of formula (8) only is formula (7) 1/N, so consumes resources greatly reduces.

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

Be experimental result explanation of the present invention below:

Table 1 does not have the phase measurement (20 ° of true phase values) of this method under the situation of making an uproar

Can find out from table 1, the phase measurement accuracy of this method is higher, under the situation of the not enough half period of observation cycle number (CiR:cycles in record), its measuring error is about 0.2 °, 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.

More than discuss and only limit to not have the situation of making an uproar.For existing under the noise situation, each the processing stage waveform as shown in Figure 9.Because the randomness of noise, 500 times Monte Carlo simulation 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)

Can find out that from table 2 this measuring method has certain anti-noise ability, its noiseproof feature can be improved by increasing spectrum exponent number N.Contrast table 1, table 2 can find out that table 1 data of noise-free case show: when spectrum exponent number N when 128 increase to 512, frequency is 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 when 128 increase to 512, frequency is that the phase measurement error of 1Hz signal then 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 the noise situations is arranged, the improvement degree of measuring accuracy is much obvious when making an uproar than nothing.

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 accumulate on one or two FFT spectral line, less than 1 situation, then only need to consider a spectral line at k=1 place for sampling period number of the present invention.And noise signal is a 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 subjected 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 big more, then σ ²/ N is more little, and it is also more little by the error that noise causes to survey phase time.

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

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

(1) calls core algorithm, finish the parameter estimation of received signal and handle;

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

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

Need point out, owing to adopted digitized method of estimation, thereby determined the complexity of Figure 10 system, in real time the principal element of degree and degree of stability is not that the periphery of DSP device among Figure 10 is connected, but the core algorithm for estimating 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 the phase measurement of interval sinusoidal signal " short " this core algorithm for estimating that is proposed in the DSP device, finishes high precision, low complex degree, phase estimation efficiently based on this.

Figure 11 flow process is divided into following several steps:

(1) at first needs according to concrete application requirements (as the concrete measurement requirement of seismic event and submarine), guestimate signal frequency f ₀Count CiR with the signal observation cycle, and set the phase accuracy requirement according to concrete needs.This step is to propose real needs from the engineering aspect, so that follow-up flow process is handled 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 influence among the)+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, move this algorithm after, 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) meet engine request until measurement result, then the output bus by DSP exports outside display drive device to, phase measurement is carried out number show.

Need point out, owing to adopted the DSP realization, make entire parameter estimate that operation becomes more flexible, can be according to the concrete condition of the various components that signal comprised, change the inner parameter setting of algorithm flexibly by programming, as the exponent number N of analysis of spectrum, the spectral line number that participation is proofreaied and correct etc.

Claims (3)

1. phase estimation method of lacking interval sinusoidal signal, CiR＜1 is counted by tested sinusoidal signal sample cycle, and this method may further comprise the steps:

Step 1 is to simulating signal x (n)=Acos (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 the analytic trnasformation to x (n) ₁(n), the concrete processing of analytic trnasformation comprises:

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

Behind the Xi Er baud conversion

Frequency spectrum

Be expressed as

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

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

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

$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)$

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

Carry out obtaining complex values Q after the x inner product

$Q=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}y\left(n\right)e^{-j\frac{2\mathrm{\&pi;}}{N}n}$

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

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

3. experimental provision of realizing the phase estimation method of the interval sinusoidal signal of weak point as claimed in claim 1, input signal is that CiR＜1 is counted by tested sinusoidal signal sample cycle, it is characterized in that, this experimental provision comprises that signal conditioning circuit, A/D change-over circuit, digital signal processor DSP and output drive and display unit, wherein, signal conditioning circuit is connected with the A/D change-over circuit, is 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 again, sampling obtains sample sequence x (n), 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 sample rate f according to actual needs _s, finish the parameter estimation of received signal and handle; Demonstrate the estimated value of frequency modulation rate and centre frequency at last by output driving demonstration and display module thereof.

Wherein, the signal parameter of digital signal processor DSP is estimated to handle, and may further comprise the steps:

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

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 among the)+B; Carry out phase measurement, obtain phase measurement.

Images