CN102175916B - Short sample dense frequency signal parameter measurement method - Google Patents

Abstract

The invention discloses a short sample dense frequency signal parameter measurement method, which comprises the following steps of: sampling an input analogue signal x(t) at an equal interval with sampling frequency fs to obtain 2N-1 discrete samples x(n); performing all-phase discrete Fourier transform (apFFT) on the x(n) to obtain a complex sequence Xa(k), obtaining amplitude spectrum by using the modulus of the Xa(k), and finding a spectral peak position k=k0; extracting a frequency range omega of [(k0-1)Deltaomega, (k0+2)Deltaomega] taking k0Deltaomega as a center, performing all-phase discrete Fourier transform (DTFT) on x(n) on the frequency range, and taking the phase spectrum of the transform, wherein Deltaomega is equal to 2phi/N; observing the characteristics of a phase spectrum curve, and determining single-frequency condition if the phase spectrum curve presents straight distribution characteristics; if the phase spectrum curve presents continuous oscillation, determining dual-frequency distribution condition, and further estimating frequency, amplitude and phase values of two dense components; in addition, determining the multi-frequency conditions of component numbers more than or equal to 3. The method has the characteristic of accuracy controllability.

Description

The measurement method of parameters of short sample dense frequencies signal

Technical field

The present invention relates to a kind of parameter measurement of frequency signal.Particularly relate to a kind ofly when collecting short sample, the sine wave signal that comprises intensive contiguous composition is carried out the measurement method of parameters of short sample dense frequencies signal of the digital measurement of frequency, initial phase and amplitude.

Background technology

The parameter measurement of multifrequency sine wave frequency, amplitude and first phase is the problem that engineering circle usually runs into.At present, the digital photogrammetry of sinusoidal signal with its precision height, scheme flexibly and need not to introduce advantage such as complicated circuit and become widely used measurement means.Digital measuring method need be sampled to signal with certain sampling rate, relends to help digital signal processing algorithm that sampled data is carried out analytical calculation, and obtains the estimated value of frequency, amplitude and first phase.Therefore, the quality of digital signal processing method has directly determined final argument estimated accuracy and measurement range.

Yet, in the engineering practice in each field,, have a very thorny difficult problem for the parameter measurement of multifrequency sinusoid.That is exactly in the multiple frequency content of signal; run into situation that two adjacent frequency contents suffer very closely (as in military radio environment through regular meeting; reception carrier is subjected to nearly enemy's strong jamming of position frequently through regular meeting); and limit by the sampling time, we usually gather less than sufficiently long sample signal.How to differentiate and measure nearly frequency composition under the short sample situation and frequency thereof, amplitude and just phase parameter be a challenging difficult problem in the engineering practice always, press for solution.

In the last few years, people also were devoted to the signal parameter Measurement Study always, and achievement in research emerges in an endless stream, such as:

(1) the counting additive process (Da Shi and bright, stone Tian Xiushu. frequency measurement circuit [P]. Japan: 008144012000.08.18.), proposed frequency measurement unit and in the cycle of displacement mutually, reference clock has been counted.Totalizer is set comes the counting addition of frequency measurement unit.By the such addition of count value that in the cycle of displacement, records, obtain frequency measurement.

(2) adaptive frequency measuring appliance method (Hu Qiannan, Cao Quanxi, Gong Hongjin, Liu Yawei. adaptive frequency measuring appliance [P] China: 2008202215832008.12.04.).This device comprises the sample conversion circuit, and the sample conversion circuit connects filtering circuit, filtering circuit linking number analog conversion circuit, and D/A converting circuit connects the digital quantity processing single chip, and the digital quantity processing single chip connects man-machine interface.Finish frequency measurement by hardware circuit than relative broad range.

(3) sampling statistic law (Yan Keguo, Jin Cunshan. a kind of frequency measurement method and device [P]. China: 2008100571922008.01.30.).This measurement mechanism comprises: receive measured signal, the edge of described measured signal triggers and interrupts; When sampling is finished, obtain according to current interruption times and to require sampling number; According to timing interruption times and timer time, obtain the sampling time; Sampling number and sampling time obtain the frequency values of described measured signal according to described the requirement.

Above various methodologies all realizes by hardware circuit, too relies on special-purpose circuit, thereby can expend great amount of hardware resources and measurement range is limited, and precision is difficult to be held, and seems powerless when running into the situation said method that two frequency contents suffer very closely.Therefore, recent people turn to software to the signal parameter measuring method gradually by hardware, realize the measurement and the analysis of parameter by software programming with mathematical method.

DFT (Discrete Fourier Transform, discrete Fourier transformation) and DTFT (Discrete Time FourierTransform, the discrete time Fourier transform) be the parameter estimation method of the most frequently used sinusoidal signal, DFT is the discrete sampling result of DTFT in frequency domain.The both can directly carry out conversion and obtain spectrogram sample sequence, directly searches out the peak value spectral position on spectrogram, distinguishes by the identification spectrum peak position and closes on composition.

Viewpoint of measures from analysis of spectrum, on the one hand, sample length has directly determined the exponent number of analysis of spectrum, short sample means that then analysis of spectrum exponent number N can't obtain too high, thereby make the spacing of two neighbour's frequency contents that signal comprises probably less than frequency resolution Δ ω=2 π/N of DFT, promptly belong to the intensive spectrum situation of short sample, obviously the intensive spectrum that is formed by short sample can increase the identification difficulty of frequency content; On the other hand, the common defects of DTFT and DFT is that two kinds of analysis of spectrums all exist very serious spectrum leakage effect.As Sanjit K.Mitra.DIGITAL SIGNAL PROCESSING:A Computer Based Approach, 3rdEdition[M] .NewYork:McGraw-Hill, 2005, Gao Xiquan, Ding Yumei. digital signal processing-third edition [M]. Xi'an: publishing house of Xian Electronics Science and Technology University, 2008 and D.Agrez.Improving phase estimation with leakageminimization[J] IEEE Trans.Instrum.Meas., 2005,54 (4): 1347-1353 is described.

And DFT is than DTFT, though discretely realize relatively easily, has fence effect and reduced the analysis of spectrum performance.For the intensive spectrum situation of short sample, the DFT spectrum of each frequency content and DTFT spectrum (comprising spectral amplitude and phase spectrum) are all suffered very closely and are caused interference (Xie Ming between very large spectrum, the fourth health. the bearing calibration [J] of two overlapping frequency spectrums of dense frequencies composition. the vibration engineering journal, 1999,12 (3): 109-114. side body lotus, flood one. utilize FFT to proofread and correct the frequency and the phase place [J] of two intensive signals. radar science and technology, 2005,3 (6): 378-382.), this directly causes DFT spectrum and DTFT spectrum to estimate almost to lose efficacy when identification dense frequencies composition.

For intensive spectrum parameter recognition problem, modern spectrum is estimated (Petre Stoica, Randolph Moses.Sepctral Analysis ofSignals[M] .Upper Saddle River, Prentice Hall, 2005.) method (and as the AR model method (Kay S M.ModernSpcetral Estimation:Theroy and Application[M] .Englewood Cliffs.NJ:Prentice-Hall, 1988.), the Maximum Entropy Spectral Estimation method (Burg J P.Maximum entropy spectral analysis[D] .Ph.D dissertation, Dept.ofGeophysics, Stanford Univ., CA, 1975.)) can't be competent at equally.Reason is that modern spectrum estimates it is the spectral analysis method of constructing on statistical significance based on model parameter, must carry out statistical study (as calculating autocorrelation function) to great amount of samples and just can draw the precise analytic model parameter; In addition, modern spectral estimation method all be meant power Spectral Estimation (Petre Stoica, Randolph Moses.Sepctral Analysis of Signals[M] .Upper Saddle River, Prentice Hall, 2005.), abandoned the phase information of signal fully.

So above the whole bag of tricks has excessively been emphasized spectral amplitude (or power spectrum), rely on the way at search spectral amplitude peak to know frequency location singlely.And for the situation of intensive composition sample deficiency, because each composition exists mutually and disturb between serious spectrum, cause spectral amplitude unimodal shape to occur possibly and can't distinguish each son spectrum.Thereby intensive spectrum estimates it is the difficult problem of signal Processing always.

For explain the difference of intensive spectrum under short sample and the long sample from spectrogram, suppose and to estimate x (t)=cos (2 π f ₁T+ θ ₁)+2cos (2 π f ₂T+ θ ₂) each frequency content, f wherein ₁=3.3Hz, f ₂=3.5Hz, its x (t) at the waveform in the 10s scope shown in Fig. 1 (a), with f _sThe sampling rate of=32Hz is respectively to x (t) sampling, gathers waveform that N=32 sample obtain shown in Fig. 1 (b), gathers waveform that N=256 sample obtain shown in Fig. 1 (c).These two sections samples are carried out DFT and DTFT respectively, and the spectrogram that obtains is respectively shown in Fig. 1 (d) and Fig. 1 (e).

DFT and DTFT spectrogram when observing cosine signal x (t) sampling number N=32 among Fig. 1 (d): because sample is short, two frequency distribution are in a resolution, spectrum leakage is serious, has only a spectrum peak, obviously we can not read the frequency content of wanting from such collection of illustrative plates, and phase information is then covered fully.DFT and the intensive spectrogram of DTFT when Fig. 1 (e) is cosine signal x (t) sampling number N=256, pairing two the spectrum peaks of two intensive compositions are distinguished in big activation among the figure, and cost is that the sampling number that expends is the former 8 times.

As can be seen from Figure 1, we can prolong the sampling time, perhaps collect more sample and improve resolution, thereby weaken the accuracy rate that frequency leakage improves parameter estimation.But, usually be to can not get abundant sample in engineering practice: such as, militarily, time is life, in the process of data processing speed and accuracy rate had high requirement again, is difficult to the sufficient enough samples of time collection of assurance and does parameter estimation; And for example on current clinical medicine, more and more will consider standing the time limit of patient, pumping signal can not continue long usually, thereby gathers less than sufficiently long useful response data, so clinical analysis of diagnosis also must be started with from a small amount of only sample usually.Therefore, we have to select other approach for the intensive spectrum of short sample, open up new observation angle and estimate intensive spectrum parameter value.

In fact, the feature of signal both had been embodied in the spectral amplitude and had been also contained in the phase spectrum.For intensive spectrum distribution situation, under the situation that the spectral amplitude analysis was lost efficacy, be necessary to consider information extraction from phase spectrum.Whole phase FFT (all-phase DFT, apFFT ^[12,13,14]) be a kind of novel spectral analysis method that proposes in recent years, document [10] points out that whole phase FFT has than the better inhibition spectrum of conventional DFT leakage characteristics, and can directly extract phase information.

Summary of the invention

Technical matters to be solved by this invention is to provide a kind of calculated amount less, estimated efficiency height, the measurement method of parameters of the short sample dense frequencies signal that resource cost is few.

The technical solution adopted in the present invention is: a kind of measurement method of parameters of lacking sample dense frequencies signal comprises the steps:

1) at first with sample frequency f _sSimulating signal x (t) to input carries out equal interval sampling acquisition 2N-1 discrete sample x (n);

1) x (n) is carried out obtaining complex sequences X after the apFFT conversion _a(k), getting its mould is worth searching out spectrum peak position k=k behind the spectral amplitude ₀

3) for shortening the frequency domain region of search, get with k ₀Δ ω is frequency range ω the ∈ [(k at center ₀-1) Δ ω, (k ₀+ 2) Δ ω], Δ ω=2 π/N does the apDTFT conversion with x (n) on this frequency range, get the phase spectrum of this conversion;

4) observe the phase spectrum curvilinear characteristic,, then be judged as the single-frequency situation, at this moment do not have intensive spectrum discrimination problem if the phase spectrum curve presents straight distribution characteristics; If the phase spectrum curve presents continuous oscillation, and have wide transitional zone feature, then be judged as the double frequency distribution situation; In addition then be into the multifrequency situation of mark more than or equal to 3 compositions;

5) if be judged as the double frequency situation, then further estimate the frequency, amplitude and the phase value that two intensive compositions.

Described apFFT conversion is with long convolution window w for (2N-1) _cTo the input data weighting, the data that will be spaced apart N then superpose in twos, except the neutral element, and form data y (0), y (1) ..., y (N-1) carries out FFT to this N new data again and promptly obtains apFFT X as a result _a(0), X _a(1) ..., X _a(N-1), w wherein _cForm for the rectangular window convolution of N by two.

Described apDTFT conversion is that 2N-1 input data are divided into N sub-segmentation x ₀, x ₁, x ₂.., x _N-1Its neutron segmentation x _m=[x (and m), x (m+1) ..., x (m+N-1)] ^T, m=0,1 ..., N-1; And then the DTFT that asks for each sub-segmentation respectively composes

M=0,1 ..., N-1; Again to each son spectrum X _m(j ω) summation that adds up promptly obtains apDTFT phase spectrum X _a(j ω).

Described apFFT conversion, be with input signal x (n) through after the quarter window weighting respectively with e ^{-j ω n}, n ∈ [N+1, N-1] multiplies each other, and the accumulation summation promptly obtains the apDTFT phase spectrum X of equivalence again _a(j ω).

To the double frequency situation, the process of then further estimating the frequency, amplitude and the phase value that two intensive compositions is as follows: will further be reduced into ω ∈ [(k between the frequency area of observation coverage ₀+ 1) Δ ω, (k ₀+ 2) Δ ω], in this interval, continue to do the extreme point search, the 1st the extreme point ω that searches out ₁The Frequency Estimation of tie element 1 then

The extreme value phase place of this frequency then is the phase estimation of composition 2

The 2nd the extreme point ω that searches out ₂The Frequency Estimation of tie element 2 then

The extreme value phase place of this frequency then is the phase estimation of composition 1

Then with the digital angular frequency value ω of two frequency contents estimating to obtain ₁, ω ₂And Integer N value substitution Amplitude Estimation formula

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;})\left]\right|}{\frac{{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;})\right]}{\frac{{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}}$

Promptly obtain the Amplitude Estimation of these two frequency contents respectively.

The beneficial effect of the measurement method of parameters of short sample dense frequencies signal of the present invention is as follows:

1, filled up specially at the intensive spectrum situation of short sample, parameters such as signal initial phase, frequency have been done the blank of effective precise Estimation Method.

The spectrum estimation of original digital signal processing is to have under the adequate sample situation mostly, does spectrum estimation by the mode of observation amplitude spectrum or spectrum peak.The signal parameter information that they are not sure and comprise in the phase spectrum, the parameter estimation problem under the intensive spectrum situation of short sample exactly can be solved by phase spectrum, and through consulting, existing document does not relate to, and the present invention has filled up this blank.

2, the inventive method has the controlled characteristics of precision.

In the practical application, can adjust resolution by appropriate change frequency sweeping step-length and come control accuracy according to the precision needs in the actual engineering.

3, calculated amount is less, the estimated efficiency height, and resource cost is few, has saved hardware cost greatly.

In parameter estimation procedure, do the peak dot search, computer capacity is narrowed down to peak dot k ₀Neighbouring is a bit of, and remainder is not processed, and has reduced a lot of unnecessary calculating, has accelerated computing velocity so greatly.

4, for intensive bispectrum parameter measurement situation, higher Detection of Weak Signals ability is arranged.

This is because the front analyzes, knowing that intensive spectrum belongs under the situation of bispectrum, for accurately drawing the frequency values of two intensive compositions, only need search for two extreme point positions of the phase spectrum of an apDTFT spectrum in the Δ ω interval and can realize, it doesn't matter with the amplitude of two compositions.Under the high s/n ratio situation, the situation for there being a weak composition around the strong composition is very beneficial for the detection to weak composition like this.

5, certain noise resisting ability is arranged.

Can show that from table 4 even for the serious interference signals of noise that is subjected to shown in Figure 11, the inventive method still can be finished accurate Frequency Estimation and more accurate phase estimation.As long as accurate centre of location frequency location and extreme point position can be made estimates of parameters fast and accurately in full phase place DTFT phase place spectrogram.

Description of drawings

Fig. 1 is the time domain waveform figure and the spectrogram of the different sample lengths of certain double frequency cosine signal;

Fig. 2 is a technical scheme process flow diagram of the present invention;

Fig. 3 is that the whole phase FFT analysis of spectrum is simplified (N=4) synoptic diagram;

Fig. 4 is the apFFT amplitude spectrum that double frequency refers to signal x (n) again;

Fig. 5 is deriving from traditional DTFT to apDTFT (N=4) synoptic diagram;

Fig. 6 is that full phase place DTFT analysis of spectrum is simplified (N=4) synoptic diagram;

Fig. 7 is direct DTFT spectrum and full phase place DTFT spectrum (N=8) oscillogram;

Wherein:

(a) the DTFT spectral amplitude that directly obtains by Fig. 3 | Y (j ω) |

(b) the DTFT phase spectrum that directly obtains by Fig. 3

(c) spectral amplitude that obtains of the apDTFT flow process of simplifying by Fig. 6 | X _a(j ω) |

(d) phase spectrum that obtains of the apDTFT flow process of simplifying by Fig. 6

Fig. 8 is that the full phase place DTFT spectrum of single-frequency, double frequency, multiple-frequency signal compares (N=16) oscillogram;

Wherein:

(a) the apDTFT phase spectrum of apDTFT spectral amplitude (b) simple signal of simple signal

(c) the apDTFT phase spectrum of apDTFT spectral amplitude (d) two-frequency signal of two-frequency signal

(e) the apDTFT phase spectrum of apDTFT spectral amplitude (f) multiple-frequency signal of multiple-frequency signal

Fig. 9 is the apDTFT spectrum of two-frequency signal, N rank DTFT spectrum and 2N rank DTFT spectrum (N=8) oscillogram;

Wherein:

(a) apDTFT spectral amplitude (b) apDTFT phase spectrum

(c) N rank DTFT spectral amplitude (d) N rank DTFT phase spectrum

(e) 2N rank DTFT spectral amplitude (f) 2N rank DTFT spectral amplitude

Figure 10 is the extreme point and the transitional zone oscillogram of apDTFT phase spectrum;

Wherein:

(a) near the apDTFT phase spectrum the wide transitional zone

(b) near the partial enlarged drawing the 1st extreme point of phase spectrum transitional zone right side

(c) near the partial enlarged drawing the 2nd extreme point of phase spectrum transitional zone right side

Figure 11 is the (signal to noise ratio snr=24dB) of phase place spectrogram under the noise;

Wherein:

(a) near the apDTFT phase spectrum the wide transitional zone

(b) the apDTFT phase spectrum after near the amplification the wide transitional zone

(c) the apDTFT phase spectrum after near the further amplification the wide transitional zone

Figure 12 is a hardware embodiment schematic diagram of the present invention;

Figure 13 is a DSP internal processes process flow diagram.

Embodiment

Make a detailed description below in conjunction with embodiment and accompanying drawing measurement method of parameters short sample dense frequencies signal of the present invention.

As shown in Figure 2, the measurement method of parameters of short sample dense frequencies signal of the present invention is at first with sample frequency f _sSimulating signal x (t) to input carries out equal interval sampling acquisition 2N-1 discrete sample x (n); X (n) is carried out obtaining complex sequences X after the apFFT conversion _a(k), getting its mould is worth searching out spectrum peak position k=k behind the spectral amplitude ₀For shortening the frequency domain region of search, get with k ₀Δ ω is frequency range ω the ∈ [(k at center ₀-1) Δ ω, (k ₀+ 2) Δ ω], Δ ω 2 π/N do the apDTFT conversion with x (n) on this frequency range, get the phase spectrum of this conversion, observe its phase spectrum curvilinear characteristic: if the phase spectrum curve presents straight distribution characteristics, then be judged as the single-frequency situation, at this moment do not have intensive spectrum discrimination problem; If the phase spectrum curve presents continuous oscillation, and have wide transitional zone feature, then be judged as the double frequency distribution situation; In addition then be into the multifrequency situation of mark more than or equal to 3 compositions.

If be judged as the double frequency situation, also can further estimate frequency, amplitude and the phase value that two intensive compositions.The estimation process is as follows: will further be reduced into ω ∈ [(k between the frequency area of observation coverage ₀+ 1) Δ ω, (k ₀+ 2) Δ ω], in this interval, continue to be extreme point search, the 1st the extreme point ω that searches out at this ₁The Frequency Estimation of tie element 1 then

Its extreme value phase place then is the phase estimation of composition 2

The 2nd the extreme point ω that searches out ₂The Frequency Estimation of tie element 2 then

Its extreme value phase place then is the phase estimation of composition 1

Then with the digital angular frequency value ω of two frequency contents estimating to obtain ₁, ω ₂And Integer N value substitution Amplitude Estimation formula

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;})\left]\right|}{\frac{{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;})\right]}{\frac{{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}}$

Promptly obtain the Amplitude Estimation of these two frequency contents respectively.

Down the internal processes and the related theoretical principle of step shown in Figure 2 are done explanation one by one, and experimental result is analyzed.

(1) whole phase FFT

Its process is as shown in Figure 3: only need with long convolution window w for (2N-1) _cTo the input data weighting, the data that will be spaced apart N then superpose in twos (except the neutral element) and form data y (0), y (1) ..., y (N-1) carries out FFT to this N new data again and promptly obtains apFFT X as a result _a(0), X _a(1) ..., C _a(N-1).W wherein _cForm for the rectangular window convolution of N by two.

The present invention needs by whole phase FFT (apFFT) determining spectrum peak position, and then the frequency sweeping scope of dwindling full phase place DTFT, reduces calculated amount.

Yan Keguo, Jin Cunshan are at " a kind of frequency measurement method and device " [P]. China: point out among 200810057192 2008.01.30, whole phase FFT has very good inhibition spectrum leakage characteristics, the former very convenient spectrum peak position search that is used for determines that according to its apFFT amplitude spectrum the approximate range of frequency is believable.Just utilize apFFT spectrum to leak little good characteristic in the technical scheme of the present invention, to x (n) do after the conversion its frequency spectrum sequence, through the peak dot search, determine integer frequency position k then ₀, then x (n) comprises k ₀Near the Δ ω dense frequencies is provided with frequency range ω ∈ [(k ₀-1) Δ ω, (k ₀+ 2) Δ ω] be the computation interval of following apDTFT conversion.We are locked in frequency range in the very little scope with the peak dot search, can significantly reduce calculated amount like this.For example: with N=8 is example, right ω ₁=3.2 Δ ω, ω ₂3.5 Δ ω, θ ₁=20 °, θ ₂Do the peak dot search for=70 °.Because | ω ₁-ω ₂|＜Δ ω, so belong to intensive spectrum distribution situation.X (n) is done the normalized amplitude spectrum that apFFT obtains | X _a(k) | as shown in Figure 4, the spectral sequence of Fig. 4 is peak dot search, spectrum peak position k as can be known ₀=3, corresponding frequency is ω ₀=3 Δ ω, thereby the frequency content of x (n) is at ω ₀Near=3 Δ ω.Thus, the interval ω ∈ of frequency range [(k is set ₀-1) Δ ω, (k ₀+ 2) Δ ω]=[2 Δ ω, 5 Δ ω], be the analystal section of next step apDTFT phase spectral analysis.

(2) full phase place DTFT analysis of spectrum process

As shown in Figure 5, (all-phase Discrete Time Fourier Transform, process apDTFT) is full phase place discrete time Fourier transform: 2N-1 input data are divided into N sub-segmentation x ₀, x ₁, x ₂..., x _N-1Its neutron segmentation x _m=[x (and m), x (m+1) ..., x (m+N-1)] ^T, m=0,1 ..., N-1; And then the DTFT that asks for each sub-segmentation respectively composes M=0,1 ..., N-1; Again to each son spectrum X _m(j ω) summation that adds up promptly obtains apDTFT spectrum X _a(j ω).

Fig. 5 process can be reduced to shown in Fig. 6 flow process: with input signal x (n) through after the quarter window weighting respectively with e ^{-j ω n}, n ∈ [N+1, N-1] multiplies each other, and the accumulation summation promptly obtains the apDTFT spectrum X of equivalence again _a(j ω).

Fig. 5 is higher than its computation complexity of Fig. 6, and this N of access sub-DTFT spectrum needs bigger memory headroom.And through after the simplification of Fig. 6 process, N DTFT computing reduced to 1 time, so only need a small amount of multiplication and addition, thus calculated amount and storage space saved greatly.Be pointed out that because what generate is continuous spectrum, so Fig. 5 and Fig. 6 in frequency variable ω all need scan and obtain with certain frequency step.

The reason that above process is introduced apDTFT is: whole phase FFT is the same with traditional FFT, is discrete spectrum, has fence effect, can't observe the spectrum distribution situation (comprising spectral amplitude and phase spectrum) in the single frequency resolution.Yet for intensive spectrum distribution situation, spectrum distributed intelligence that frequency resolution is interior exactly of being concerned about, so apDTFT spectral analysis method that the present invention proposes, this novel spectral analysis method can obtain very little spectral amplitude of leakiness and phase spectrum intuitively on full frequency band, be highly susceptible to extracting intensive spectrum information.

Attention apDTFT analysis of spectrum can not directly be drawn the output X that Fig. 3 obtains from Fig. 3 _a(k) disperse, promptly can only ω=k Δ ω=k2 π/N (k=0,1 ..., observe the spectrum value on N-1).For obtaining continuous spectrum, way of directly expecting easily be exactly among Fig. 3 after phase place pre-service full length be that data y (0)～y (N-1) of N directly carries out traditional DTFT, that is:

Yet so direct way is wrong.Below illustrate its mistake.

With signal

N ∈ [N+1, N-1], ω ₀=3.3 Δ ω, Δ ω=2 π/N, θ ₀=20 °, A=1 is an example.The spectral amplitude of the Y (j ω) that N data y (0)～y (N-1) directly calculates that generates by Fig. 3 | Y (j ω) | shown in Fig. 7 (a), its phase spectrum

Shown in Fig. 7 (b).

(b) can find out from Fig. 7 (a), though right respectively on ω=k Δ ω | Y (j ω) | and Can obtain the apDFT spectral amplitude after the sampling | X _a(k) | and phase spectrum

But on other frequencies of ω ≠ k Δ ω, | Y (j ω) | spectrum is leaked very big,

Can not directly reflect 20 ° of true first phase values.Thereby Fig. 7 (a) DTFT spectrum (b) is not the apDTFT spectrum.

After Fig. 6 flow processing, can obtain the apDTFT spectral amplitude after the normalization shown in Fig. 7 (c) | X _a(j ω) | and the phase spectrum shown in Fig. 7 (d)

(b) the same with Fig. 7 (a), right respectively | X _a(j ω) | and phase spectrum

After carrying out equal interval sampling, promptly get the apDFT spectral amplitude | X _a(k) | and phase spectrum What Fig. 7 (c) was different with Fig. 7 (a) is, because each height spectrum positive and negative cancelling out each other in additive process, on other all frequencies of ω ≠ k Δ ω, | X _a(j ω) | the spectrum leakage part very little, in addition on full frequency band, Fig. 7's (d)

All equal 20 ° of true first phase values.

(3) based on the intensive spectrum type identification of full phase place DTFT phase spectrum

This intensive spectrum discrimination process is as follows: observe the phase spectrum curvilinear characteristic that is obtained by apDTFT: if the phase spectrum curve presents straight distribution characteristics, then be judged as the single-frequency situation, at this moment do not have intensive spectrum discrimination problem; If the phase spectrum curve presents continuous oscillation, and have wide transitional zone feature, then be judged as the double frequency distribution situation; In addition then be into the multifrequency situation of mark more than or equal to 3 compositions.

Give an example comparison signal to comprise the difference of the full phase place DTFT phase spectrum under single-frequency, double frequency, three kinds of situations of multifrequency.Suppose that single-frequency, double frequency, three signals of multifrequency are respectively x ₁(n), x ₂(n), x ₃(n):

$x_1\left(n\right)=A_1{e}^{j({\mathrm{\&omega;}}_{1}n+{\mathrm{\&theta;}}_1)},$ ω ₁＝(k ₀+Δk ₁)Δω；

$x_2\left(n\right)=A_1{e}^{j({\mathrm{\&omega;}}_{1}n+{\mathrm{\&theta;}}_1)}+A_2{e}^{j({\mathrm{\&omega;}}_{2}n+{\mathrm{\&theta;}}_2)},$ ω ₁＝(k ₀+Δk ₁)Δω，ω ₂＝(k ₀+Δk ₂)Δω；

$x_3\left(n\right)=A_1{e}^{j({\mathrm{\&omega;}}_{1}n+{\mathrm{\&theta;}}_1)}+A_2{e}^{j({\mathrm{\&omega;}}_{2}n+{\mathrm{\&theta;}}_2)}+,A_3{e}^{j({\mathrm{\&omega;}}_{3}n+{\mathrm{\&theta;}}_3),}$ ω ₁＝(k ₀+Δk ₁)Δω，ω ₂＝(k ₀+Δk ₂)Δω，

ω ₃＝(k ₀+Δk ₃)Δω；

N=16 wherein, Δ ω=2 π/16; A ₁=1, A ₂=2, A ₃=2; ω ₁=3.2 Δ ω, ω ₂=3.5 Δ ω, ω ₃=3.8 Δ ω.

Simplification analysis process according to Fig. 6 is made apDTFT and apFFT to these three signals respectively, and spectral amplitude that obtains and phase spectrum are as shown in Figure 8.

Here signal branch single-frequency composition, double frequency composition and three kinds of situations of multifrequency composition are discussed.

Because each intensive composition frequency gap is too small, all squeeze in size is the minizone of Δ ω=2 π/N, and the intensity of each composition each is variant, so the apDTFT spectrum is the same with traditional DTFT spectrum, only have a spectrum peak (as Fig. 8 (a) (c) shown in each spectral amplitude of (d)).The simple spectral amplitude that relies on goes to distinguish then complete failure of these three kinds of situations.Thereby being necessary to consider phase spectrum, the phase spectrum feature to three kinds of frequency distribution is described below.

(i) single-frequency situation

The phase spectrum feature of apDTFT clearly, shown in Fig. 8 (b): phase spectrum

On the full rate axle, present a straight line.On Fig. 8 (b), can directly read the phase value θ of simple signal ₁=20 °, by the readable frequency content ω that gets in simple signal spectral amplitude spectrum peak ₁=3.2 Δ ω.Just obtain amplitude parameter A=1 of signal at last through simple computation.

The (ii) intensive spectrum distribution situation of double frequency

Can find out that from the phase place spectrogram of Fig. 8 (d) phase spectrum of double frequency dense frequencies signal has it self characteristics: wide transitional zone, extreme point period profile, have only two different extreme values (maximal value, a minimum value), shown in phase spectrum among the figure.We also judge whether double frequency of a signal with these characteristics.

Because the phase place spectrogram presents period profile and local wide transitional zone feature, thereby when intensive spectrum type identification, need not frequency sweeping in whole frequency separation, only need at ω ∈ [(k ₀-1) Δ ω, (k ₀+ 2) Δ ω] in carry out frequency sweeping and get final product.Can find out from Fig. 8 (d), at this interval ω ∈ [(k ₀-1) Δ ω, (k ₀+ 2) Δ ω]=[2 Δ ω, 5 Δ ω] in, phase spectrum possesses following feature simultaneously: (1) is at ω=k ₀Present wide transitional zone characteristic near Δ ω=3 Δ ω; The periodicity of (2) blocking, promptly the width at ω ∈ [2 Δ ω, 5 Δ ω] is 3 extreme values should occur in the 3 Δ ω, because the influence of middle wide transitional zone is kept to 2 maximum value and 2 minimal values, two very big or minimizing 2 Δ ω that are spaced apart.If meet above feature, can predicate the intensive spectrum of double frequency.

(iii) multifrequency situation

If the apDTFT phase spectrum that obtains Do not meet the phase spectrum feature of single-frequency described above and double frequency composition, shown in Fig. 8 (f), just can assertive signal comprise the dense frequencies composition more than 3 or 3.

(4) parameter in assessing of two of the bispectrum situation composition frequencies, amplitude and phase places

Its process is as follows: after judging that intensive spectrum is distributed as the bispectrum situation, will further be reduced into ω ∈ [(k between the frequency area of observation coverage ₀+ 1) Δ ω, (k ₀+ 2) Δ ω], in this interval, continue to be extreme point search, the 1st the extreme point ω that searches out at this ₁The Frequency Estimation of tie element 1 then Its extreme value phase place then is the phase estimation of composition 2

The 2nd the extreme point ω that searches out ₁The Frequency Estimation of tie element 2 then

Its extreme value phase place then is the phase estimation of composition 1

Then with the digital angular frequency value ω of two frequency contents estimating to obtain ₁, ω ₂And Integer N value substitution Amplitude Estimation formula

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;})\left]\right|}{\frac{{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;})\right]}{\frac{{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}}$

Promptly obtain the Amplitude Estimation of these two frequency contents respectively.

Suppose that two are positioned at ω=k ₀Near the Δ ω dense frequencies is respectively ω ₁=(k ₀+ Δ k ₁) Δ ω, ω ₂=(k ₀+ Δ k ₂) Δ ω, k ∈ z, ω might as well be established in k ≠ 0 ₁＜ω ₂, 0＜Δ k is then arranged ₁＜Δ k ₂＜1.Can prove, at ω=(k+k ₀+ Δ k ₁) Δ ω, k ≠ 0 place, the phase spectrum value of apDTFT is the first phase value θ of second frequency content ₂Similarly, at ω=(k+k ₀+ Δ k ₂) Δ ω, k ≠ 0 place, the phase spectrum value of apDTFT is the first phase value θ of first frequency content ₁

For measuring conveniently, knowing amplitude peak position k ₀After, make k=1, then at ω=(1+k ₀+ Δ k ₂) the Δ ω place phase spectrum value of reading promptly can be used as ω ₁The first phase value θ of composition ₁Estimate, at ω=(1+k ₀+ Δ k ₁) the Δ ω place phase spectrum value of reading promptly can be used as ω ₁The first phase value θ of composition ₂Estimate.

Thereby knowing under the intensive spectrum distribution situation of double frequency, and to have judged the peak value position of spectral line be k=k ₀Situation under, for further estimating two frequencies omega of bispectrum composition exactly ₁And ω ₂Value can further be dwindled the frequency sweeping interval, only needs at ω ∈ [(k ₀+ 1) Δ ω, (k ₀+ 2) Δ ω] scanning in get final product.

Be distributed as after bispectrum distributes identifying spectrum,, need be further analyzed for estimating intensive spectrum parameter.For its principle is described, we do the analysis of spectrum contrast with full phase place DTFT and traditional DTFT, to two-frequency signal x ₂(n) carry out spectral amplitude that 8 rank apDTFT obtain shown in Fig. 9 (a), phase spectrum is shown in Fig. 9 (b); To two-frequency signal x ₂(n) carry out spectral amplitude that 8 rank DTFT obtain shown in Fig. 9 (c), phase spectrum is shown in Fig. 9 (d); To two-frequency signal x ₂(n) carry out spectral amplitude that 8 rank DTFT obtain shown in Fig. 9 (e), phase spectrum is shown in Fig. 9 (f);

Can find out from the amplitude spectrogram of Fig. 9, because the frequency interval of two compositions has only 0.3 Δ ω, so the apDTFT of Fig. 9 (a) and 8 rank of Fig. 9 (c) tradition DTFT spectral amplitude all have only 1 spectrum peak, even traditional DTFT spectrum exponent number is increased 1 times (corresponding number of samples that expends will more than apDTFT), the DTFT spectral amplitude peak, 16 rank of Fig. 9 (e) still has only 1, and spectrum leak do not have be improved significantly.So rely on the spectral amplitude analysis result, be to distinguish signal to comprise one, two or a plurality of frequency content on earth.

Thereby have to discern intensive spectrum composition by phase spectrum.Phase spectrum Fig. 9 (d) of tradition DTFT and Fig. 9 (f) are relatively more disorderly, do not have obvious characteristic.Yet we can find out that from the apDTFT phase place spectrogram of Fig. 9 (b) it has periodically, at the centre frequency place tangible transitional zone is arranged, and has only two extreme values (maximal value, a minimum value).It is bispectrum for these character representations of phase spectrum, therefore is fit to use the method for parameter estimation among the present invention to estimate.The apDTFT phase place spectrogram of Fig. 9 (b) shows: near ω=3 Δ ω

The transitional zone that has a broad is so can judge that the integer frequency position of intensive spectrum is k ₀=3.This point can more clearly be observed from Figure 10.

Figure 10 (a) is the transitional zone of two-frequency signal apDTFT phase place spectrogram, amplifies two extreme points on transitional zone right side are local, gets Figure 10 (b) (c).Figure 10 (b) observes first extreme point on transitional zone right side at ω=4.2 Δ ω places, and to have analyzed the theoretical value of this position be ω=(1+k in the front ₀+ Δ k ₁) Δ ω, thereby can determine Δ k ₁=0.2, also can read this frequency from Figure 10 (b)

Further can know θ by inference ₂=70 °; In addition Figure 10 (c) also can observe the transitional zone right side second extreme point position at ω=4.5 Δ ω places, and to have analyzed the theoretical value of this position be ω=(1+k in the front ₀+ Δ k ₂) Δ ω, thereby can determine Δ k ₂=0.5; Also can read this frequency from Figure 10 (c) So can know θ by inference equally ₁=20 °.At last, determine θ ₁=20 °, ω ₁=(k ₀+ Δ k ₁) Δ ω=3.2 Δ ω, θ ₂=70 °, ω ₂=(k ₀+ Δ k ₂) Δ ω=3.5 Δ ω.Can find that from above analytic process in the process of the frequency values of determining each composition, only needing to seek the extreme point position can realize, that is to say that the definite of frequency values of two intensive compositions do not influenced by the amplitude of two compositions fully.

Measure two frequency values ω according to above process ₁And ω ₂After, the value substitution following formula of exponent number N and Δ ω=2 π/16 just can be calculated the amplitude of two frequency contents, promptly

${\hat{A}}_1=\frac{\left|X_a\right[j({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;})\left]\right|}{\frac{{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_2+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_1)N/2]}},{\hat{A}}_2=\frac{X_a\left[j({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;})\right]}{\frac{{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}{{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000046274090000032.png,HDA0000046274090000041.png"\; state="\{\{state\}\}">N</figure-callout>}}^2{\sin }^2[({\mathrm{\&omega;}}_1+\mathrm{\&Delta;\&omega;}-{\mathrm{\&omega;}}_2)N/2]}}$

From Fig. 9 (a) amplitude spectrogram, can read | X _a[j (ω ₂+ Δ ω)] |=0.0262, | X _a[j (ω ₁+ Δ ω)] |=0.1131, thus can calculate A according to following formula ₁=1, A ₂=2.

So far, finish intensive two-frequency signal x ₂(n) the isoparametric estimation of frequency, initial phase and amplitude.From whole process method for parameter estimation of the present invention as can be seen the intensive spectrum of short sample there is very strong specific aim, can makes precise parameters to the single-frequency two-frequency signal rapidly and accurately and estimate, and whether judge its frequency content more than 2.Simultaneously, in different occasions, can pass through comprehensive adjustment frequency range ω length and scanning step control survey precision according to different requirements to degree of accuracy.

Therefore the present invention is a saving calculated amount, the controlled good method for parameter estimation of precision, has very strong engineering application.

(5) experimental result

During engineering was used, signal always was subjected to noise, also may be that the energy difference of two adjacent intensive compositions is very big.In addition, the signal that runs in the actual engineering can not be desirable complex exponential signal usually, but real cosine signal (the complex exponential signal composition that includes two conjugation).Thereby be necessary to study by experiment to compare noiseless respectively and single-frequency complex exponential signal under the noise conditions, the frequency of intensive complex exponential signal of double frequency and the intensive cosine signal of double frequency and the difference of first phase estimated accuracy are arranged.

Experiment 1:

To signal

$x\left(n\right)=e^{j({\mathrm{\&omega;}}_{0}n+{\mathrm{\&theta;}}_0)},$ ω ₀＝3.3Δω，Δω＝2π/N，θ ₀＝20°； $y\left(n\right)=e^{j({\mathrm{\&omega;}}_{1}n+{\mathrm{\&theta;}}_1)}+2e^{j({\mathrm{\&omega;}}_{2}n+{\mathrm{\&theta;}}_2)};$

z(n)＝cos(ω ₁n+θ ₁)+2cos(ω ₂n+θ ₂)，ω ₁＝3.2Δω，ω ₂3.5Δω，θ ₁＝20°，θ ₂＝70°。

Carry out the estimation of frequency and phase place by the inventive method under nothing is made an uproar situation, table 1 has provided different spectral factorization exponent number N situation corresponding parameters estimated results.

Table 1 does not have the parameter measurements (scanning step: Δ ω/10 of this method under the situation of making an uproar ³, Δ ω=2 π/N)

From table 1, can find out, with the parameter measurement that method of the present invention is carried out,, composing exponent number N=8 for single-frequency complex exponential signal x (n), when promptly only having used 2N-1=15 sample, just can obtain 20 ° of accurate Frequency Estimation result 3.3 Δ ω and phase estimation results; For double frequency complex exponential signal y (n),, when promptly only having used 2N-1=15 sample, also can obtain 20 °, 70 ° of accurate Frequency Estimation result 3.2 Δ ω, 3.5 Δ ω and initial phase estimated results at spectrum exponent number N=8; For double frequency cosine signal z (n), because there is the phase mutual interference of two sidebands in cosine signal, thereby its precision reduces (the frequency estimation error only is 0.03 Δ ω, and the phase place estimation error only is about 2 °) slightly than the complex exponential situation, but the measuring error of its frequency content is in 10 ^-2Resolution Δ ω also usually is satisfactory in engineering practice.

Experiment 2:

Scanning step is changed into: Δ ω/1007, Δ ω=2 π/N, repeated experiments 1, its experimental data that obtains is as shown in table 2.

Table 2 does not have the parameter measurements (scanning step: Δ ω/(1007), Δ ω=2 π/N) of this method under the situation of making an uproar

The measurement result that table 2 obtains has certain deviation with respect to table 1, and this is that sweep frequency ω does not traverse the cause of desirable frequency location because scanning step limits.Comparison sheet 1 can be reached a conclusion with the experimental result in the table 2, and behind definite Δ ω, the frequency sweeping step-length preferably is chosen as Δ ω/10 ⁿ, n ∈ z is beneficial to precision analysis.Here note the corresponding relation of Δ ω and analog frequency resolution, promptly Δ ω corresponding simulating frequency resolution is Δ f=f _s/ N is promptly determined jointly by sample frequency and analysis of spectrum exponent number.

Experiment 3

Two-frequency signal parameter measurement in the above-mentioned experiment all is to carry out under the little situation of the amplitude difference of two kind frequency contents, when the amplitude of two frequency contents differs greatly, the conventional vibration amplitude spectrum leakage of strong composition can occur and cover weak signal composition phenomenon, causes the big strong signal of amplitude to be easy to cover the parameter information of the little weak signal of amplitude.

For this reason, estimate the frequency and the amplitude parameter value of following signal

$y\left(n\right)=A_1{e}^{j({\mathrm{\&omega;}}_{1}n+{\mathrm{\&theta;}}_1)}+A_2{e}^{j({\mathrm{\&omega;}}_{2}n+{\mathrm{\&theta;}}_2)}{e}^{,}$ z(n)＝A ₁cos(ω ₁n+θ ₁)+A ₂cos(ω ₂n+θ ₂)

Table 3 has write down with methods analyst of the present invention not to be had under the situation of making an uproar, A ₁, A ₂Frequency and amplitude testing result with bispectrum signal of different amplitude differences.

Table 3 does not have the parameter measurements (scanning step: Δ ω/10 of this method under the situation of making an uproar ³, Δ ω=2 π/N, N=16)

Can find out from table 3, for the double frequency complex exponential signal, A no matter ₂Be A ₁2 times, 100 times or even 1000 times, this method can both very accurately estimate frequency and amplitude parameter value; For the double frequency cosine signal, because influencing each other between two sidebands of cosine signal, than the complex exponential situation, its estimated accuracy decreases, and wherein the weak signal composition parameter is lower than the parameter estimation precision of strong signal content: at A ₂Be A ₁1～15 times scope in, its weak composition frequency estimation precision can be controlled in the 0.1 Δ ω, and the Amplitude Estimation error can be controlled at about 10%, can meet the engineering practice requirement substantially, and the basic bias free of the Frequency Estimation of its strong signal content, the Amplitude Estimation error is also very little.

Why the present invention can still can finish parameter in assessing under the very big situation of two intensive composition capacity volume variances, be because the cause that the present invention estimates according to phase spectrum, has avoided classic method only according to the defective of spectral amplitude when the Testing of Feeble Signals.

Experiment 4

More than discuss and only limit to not have the situation of making an uproar.For there being the noise situation, be example with double frequency complex exponential signal and double frequency cosine signal, still carry out parameter detecting and error analysis to testing 1 signal with the inventive method.Wherein the apDTFT phase place spectrogram of the intensive spectrum signal y of double frequency (n) as shown in figure 11.

Experiment institute plus noise is a white Gaussian noise, signal to noise ratio (S/N ratio) 24dB.As can be seen from Figure 11, this noise has produced influence to the phase place spectrogram, but transitional zone still clearly, can read relevant frequency deviation value and initial phase from Figure 11.Experiment has also carried out having the single-frequency complex exponential signal under the condition of making an uproar, the parameter estimation of double frequency cosine signal, inserts table 4 after the estimated result statistics to three kinds of signal with different type.

The parameter estimation result of this method under table 4 noise situations (scanning step is a pi/2 000, signal to noise ratio snr=24dB)

Can find out that from table 4 this measuring method has certain anti-noise ability, show in the estimation to frequency component.For the double frequency complex exponential signal, when spectrum exponent number N=8 (number of samples is 15), its weak composition estimated frequency error is about 0.07 Δ ω, and the estimated frequency error of strong composition is because self-energy is strong, and it is big to resist the noise ability, and its error only is 0.003 Δ ω.When spectrum exponent number N=16 (number of samples is 31), the Frequency Estimation precision of strong and weak signals further improves.For the double frequency cosine signal, owing to have the phase mutual interference of two sidebands, thereby under same exponent number, Frequency Estimation ratio of precision double frequency complex exponential signal is lower.

In addition, can find out also from table 4 that the Amplitude Estimation ratio of precision Frequency Estimation precision of all situations is all low, this is easy to from formula (4) as can be seen, and amplitude is calculated by Frequency Estimation, thereby the error of Frequency Estimation can be updated in the Amplitude Estimation and goes.

Give simple declaration to implementing hardware of the present invention below.As shown in figure 12, sampling obtains sample sequence x (n) to measured signal through A/D (analog-to-digital conversion device), and 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 13 (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 12 ₂Clock frequency), makes it the frequency resolution Δ ω under this sampling rate condition, satisfy actual needs as best one can;

Export to during (3) with frequency, amplitude and phase estimation fructufy and drive and display module.

Need point out, owing to adopted digitized method of estimation, thereby determined the complexity of Figure 12 system, in real time the principal element of degree and degree of stability is not that the periphery of DSP device among Figure 12 is connected, but the core algorithm for estimating that the DSP internal program memory is stored.

Be pointed out that involved in the present invention is, and short sample parameter is measured, thus do not need very big storage space storage of collected to sample data (only needing the space of the data internal RAM of DSP).Why Figure 12 connects external RAM is because under sample size situation seldom, expect enough accurate spectrogram, the frequency sweeping point must be obtained very intensively, so just needs the data after the certain storage space of consumption is used to deposit the apDTFT processing.

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

The present invention implants " the intensive spectrum method for parameter estimation of novel short sample " this core algorithm for estimating that is proposed in the DSP device, finishes high precision, low complex degree, frequency, amplitude and phase parameter estimation efficiently based on this.

Figure 13 flow process is divided into following several steps:

(1) at first need according to concrete application requirements (as the concrete measurement requirement of medical science and military affairs etc.), the concrete needs of guestimate signal frequency resolution ax/ω and basis are set the phase accuracy requirement.This step is to propose real needs from the engineering aspect, so that follow-up flow process is handled targetedly.

(2) then, the CPU primary controller enters internal RAM from I/O port reads sampled data.

(3) follow-up " DC processing " is in order to eliminate the influence of the flip-flop in the measured signal.Otherwise the existence of flip-flop can reduce measuring accuracy.Flip-flop is easy to measure, and the mean value that only need calculate sampling point can obtain.

(4) carrying out phase measurement by Fig. 2 processing procedure of the present invention is the most crucial part of DSP algorithm, move this algorithm after, can obtain frequency, amplitude and the phase measurement of intensive composition.

(5) judge whether the inventive method satisfies engineering demand, if do not satisfy, program is returned, again setpoint frequency resolution as requested.

(6) meet engine request until measurement result, the output bus by DSP exports outside display drive device to then, frequency, amplitude and the phase measurement of intensive composition 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 (4)

1. a measurement method of parameters of lacking sample dense frequencies signal is characterized in that, comprises the steps:

1) at first with sample frequency f _sSimulating signal x (t) to input carries out equal interval sampling acquisition 2N-1 discrete sample x (n);

2) x (n) is carried out obtaining complex sequences X after the apFFT conversion _a(k), getting its mould is worth searching out spectrum peak position k=k behind the spectral amplitude ₀

3) for shortening the frequency domain region of search, get with k ₀Δ ω is frequency range ω the ∈ [(k at center ₀-1) Δ ω, (k ₀+ 2) Δ ω], Δ ω=2 π/N does the apDTFT conversion with x (n) on this frequency range, get the phase spectrum of this conversion;

4) observe the phase spectrum curvilinear characteristic,, then be judged as the single-frequency situation, at this moment do not have intensive spectrum discrimination problem if the phase spectrum curve presents straight distribution characteristics; If the phase spectrum curve presents continuous oscillation, and have wide transitional zone feature, then be judged as the double frequency distribution situation; In addition then be into mark more than or equal to 3 multifrequency situation;

5) if be judged as the double frequency situation, then further estimate the frequency, amplitude and the phase value that two intensive compositions.

2. the measurement method of parameters of short sample dense frequencies signal according to claim 1 is characterized in that, described apFFT conversion is with long convolution window w for (2N-1) _cTo the input data weighting, the data that will be spaced apart N then superpose in twos, except the neutral element, and form data y (0), y (1) ..., y (N-1) carries out FFT to this N new data again and promptly obtains apFFT X as a result _a(0), X _a(1) ..., X _a(N-1), w wherein _cForm for the rectangular window convolution of N by two.

3. the measurement method of parameters of short sample dense frequencies signal according to claim 1 is characterized in that, described apDTFT conversion is that 2N-1 input data are divided into N sub-segmentation x ₀, x ₁, x ₂..., x _N-1, its neutron segmentation x _m=[(m), x is (m+1), for x ..., x (m+N-1)] ^T, m=0,1 ..., N-1; And then the DTFT that asks for each sub-segmentation respectively composes

, m=0,1 ..., N-1; Again to each son spectrum X _m(j ω) summation that adds up promptly obtains apDTFT phase spectrum X _a(j ω).

4. the measurement method of parameters of short sample dense frequencies signal according to claim 1 is characterized in that, described apFFT conversion, be with input signal x (n) through after the quarter window weighting respectively with e ^{-j ω n}, n ∈ [N+1, N-1] multiplies each other, and the accumulation summation promptly obtains the apDTFT phase spectrum X of equivalence again _a(j ω).

Images