CN102215187B - Multistage stationary signal frequency domain related method of frequency estimation - Google Patents

Abstract

The invention relates to the field of signal processing, and in particular relates to a multistage stationary signal frequency domain related method of frequency estimation. The method is applicable to M stages of stationary signals, wherein M is an integer greater than 1, the length of each stage of signal is arbitrary, the frequency is constant, and the frequency difference of any two stages of signals is known. The method in the invention comprises the following steps: firstly, designing the weighting factor to perform weighting accumulation on M stages of normalized spectrums to obtain a weighting accumulation spectrum matrix; secondly, searching for the weighting accumulation spectrum matrix through a peak to obtain the best weighting accumulation spectrum; then, conducting corresponding operations on the best weighting accumulation spectrum and the accumulation spectrum of M stages of normalized spectrums, to obtain related spectrums of the frequency domain; and finally, searching related spectrums of the frequency domain according to the peak, to obtain the high-accuracy frequency estimation value. The frequency estimation method in the invention has the advantages of high accuracy and favorable universality.

Description

A kind of multistage stationary signal domain related method of frequency of Frequency Estimation

Technical field

The present invention relates to a kind of multistage stationary signal domain related method of frequency of signal process field, particularly Frequency Estimation.

Background technology

The multistage stationary signal refers to the constant coherent signal of some band frequencies that same target or related object are repeatedly sampled and obtained.For example, the intermediate-freuqncy signal of a plurality of piecewise linearity frequency modulated continuous wave radars was tiltedly processed, after low-pass filtering, can be formed multistage with sinusoidal signal (being the sinusoidal signal that the multistage frequency is identical) frequently through the past.From information view, the multistage signal has the amount of information that is several times as much as the single hop signal, be easy to obtain, it is one of effective way improved the signal processing accuracy that the multistage signal is carried out to fusion treatment, can be widely used in Frequency Estimation, the estimation of LFM (Linear Frequency Modulated) signal parameter, the numerous areas such as Frequency Hopping Signal parameter Estimation and VCO (Voltage Controlled Oscillator) nonlinearity correction, have important research meaning and using value.

At present, frequency estimating methods based on the multistage stationary signal mainly contains: and (1) card cloth distribution (list of references [1]: Becker K.New algorithm for frequency estimation from short coherent pulses of a sinusoidal signal[C] .Radar and Signal Processing, IEEE Proceedings F, Aug 6,1990:(4): 283-288).(2) Yule-Walker equation (list of references [2]: Sill J A, Black Q R.Frequency estimation from short pulses of sinusoidal signals[C] .Military Communications Conference, Conference Proceedings, IEEE, 21-24 Oct, 1996, McLean, VA, USA, 1996 (3): 979-983).(3) Spectrum Averaging Method (list of references [3]: Liu Liangbing. the information fusion method of Frequency Estimation and application thereof [D]. Chongqing: Logistics Engineering College's doctorate paper, in June, 2008: 12).(4) the phase accumulation method (list of references [4]: build the Meng. the phase association technology [J] of block sampling signal. system engineering and electronic technology, 2004,26 (12): 1784-1786,1797).(5) the isometric signal fused algorithm of multistage frequency reducing (list of references [4]: Liu Liangbing, Tu Yaqing. a kind of frequency estimating methods [J] based on the isometric signal fused of multistage frequency reducing. information and control, 2008,37 (4): 403-407).But the problems such as the ubiquity universality is poor, estimated accuracy is lower.For example, card cloth distribution, Yule-Walker equation, Spectrum Averaging Method and phase accumulation method all can not be estimated the not frequency of homogenous frequency signal of multistage, can only estimate the frequency of multistage homogenous frequency signal, and estimated accuracy are not high, and noise immunity is poor; Although the isometric signal fused algorithm of multistage frequency reducing can be estimated the not frequency of homogenous frequency signal of multistage, and noise immunity is better, precision is higher, can only process the isometric signal of multistage, can not process the not isometric signal of multistage.

To sum up, multistage signal fused processing method has important research meaning and using value, but there are problems in existing method, and a kind of precision of needs proposition is higher, the better frequency estimating methods of universality.

Summary of the invention

The objective of the invention is to propose that a kind of precision is higher, the better frequency estimating methods of universality, solve the subject matter that existing multistage signal fused processing method exists, be applicable to multistage length stationary signal arbitrarily, expand its range of application.

In order to achieve the above object, design of the present invention is: at first, the design weighted factor is weighted accumulation to M section normalization frequency spectrum, obtains weighting accumulation spectral matrix, and M is greater than 1 integer; Secondly, by spectrum peak search weighting accumulation spectral matrix, obtain optimum weighting accumulation frequency spectrum; Then, the cumulative frequency spectrum of optimum weighting accumulation frequency spectrum and M section normalization frequency spectrum is carried out to related operation, obtain the spectrum correlation spectrum; Finally, spectrum peak search spectrum correlation spectrum, obtain high-precision frequency estimation.

Based on above-mentioned design, the present invention adopts following technical scheme:

1, a kind of multistage stationary signal domain related method of frequency of Frequency Estimation, applicable object is M section stationary signal, and M is greater than 1 integer, and the length of every segment signal is any, frequency is constant, and the difference on the frequency of any two segment signals is known, and the method comprises the following steps:

1) structure normalization matrix C carries out respectively normalized to the frequency spectrum of M section stationary signal, obtains M section normalization frequency spectrum; In normalization matrix C the (m, the Elements C of a) locating (m, a) calculated by following formula,

C(m，a)＝f _A(a)+d(m)

In formula, m ∈ [1, M], f _amean a sequence, by frequency f to be estimated ₀span f _scopelinear decile (A-1) part generates, and A is greater than 1 integer, f _scope=[f _min, f _max], f _a(a) mean f _ain a element, a ∈ [1, A]; D (m) means x _mwith x ₁between difference on the frequency, x _mmean m segment signal in M section stationary signal;

The normalization frequency spectrum is calculated by following formula,

$X_m\left[f_A\left(a\right)\right]=0.5\underset{n_m=1}{\overset{N_m}{\mathrm{\&Sigma;}}}{x}_m\left(n_m\right)e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000510556000011.png"\; state="\{\{state\}\}">j</figure-callout>}{2\mathrm{\&pi;n}}_{m}C(m,a)/f_s}$

In formula, X _m[f _a(a)], N _m, n _mand f _smean respectively x _mnormalization frequency spectrum, sampling number, time series number and sample frequency, N _mfor positive integer, n _mfor natural number, f _svalue meets nyquist sampling theorem;

It is characterized in that:

2) M section normalization frequency spectrum is added up, obtained cumulative frequency spectrum X[f _a(a)], X[f _a(a)] by following formula, calculated,

$X\left[f_A\left(a\right)\right]=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{abs}\left\{X_m\right[f_A\left(a\right)\left]\right\}$

In formula, abs means the mould of calculated complex;

3) design weighted factor e ^{-jD (m, a, b)}to step 1) in M section normalization frequency spectrum be weighted accumulation, obtain weighting accumulation spectral matrix X ' _b[f _a(a)], D in weighted factor (m, a, b) is calculated by following formula,

D(m，a，b)＝θ _z(m，a _b)-θ _z(1，a _b)-g _b(a)[(2m-1)N _m+1]

X ' _b[f _a(a)] by following formula, calculated

${X^{\&prime;}}_b\left[f_A\left(a\right)\right]=\mathrm{abs}\left\{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\right\{e^{-\mathrm{jD}(m,a,b)}{X}_m\left[f_A\left(a\right)\right]\left\}\right\}$

In formula, a _b∈ [1, A], b ∈ [1, B], B is greater than 1 integer, θ _z(m, a _b)=angle{X _m[f _a(a _b)]+W[f _a(a _b)], angle means the phase place of calculated complex, f _a(a _b) be illustrated in sequence f _ain with f _b(b) be worth immediate element, f _b(b) mean sequence f _bin b element, f _bby frequency f to be estimated ₀span f _scopelinear decile (B-1) part generates, W[f _a(a _b)] mean that noise is at Frequency point f _a(a _b) the positive frequency part of frequency spectrum, g _b(a)=π [f _b(b)-f _a(a)]/f _s;

4) the weighting accumulation spectral matrix spectrum peak search step 3), its peak value element column is optimum weighting accumulation frequency spectrum;

5) by step 2) in cumulative frequency spectrum and step 4) in optimum weighting accumulation frequency spectrum carry out related operation, obtain spectrum correlation and compose;

6) spectrum peak search step 5) the frequency domain Correlated Spectroscopy, calculate f as follows ₀estimated value

In formula, H _maxthe abscissa that means spectrum correlation spectrum peak element, K means the length of spectrum correlation spectrum.

The frequency estimating methods precision the present invention relates to is high, universality good, is applicable to multistage length stationary signal arbitrarily.

The accompanying drawing explanation

Below with concrete enforcement, the present invention is further elaborated with reference to the accompanying drawings.With M=4, signal type be take sinusoidal signal and is described as example.

Fig. 1 is the algorithm basic thought.

In figure: 1 means M section sinusoidal signal, mainly comprises that the M section is with the isometric signal of frequency, the same not isometric signal frequently of M section, the different isometric signals frequently of M section, the not isometric signals of M section difference frequency; 2 mean M section normalization frequency spectrum; 3 mean cumulative frequency spectrum, and 4 mean weighting accumulation spectral matrix; 5 mean optimum weighting accumulation spectral matrix; 6 mean the spectrum correlation spectrum; 7 mean frequency estimation; 8 mean the normalization process; 9 mean cumulative process; 10 mean the weighting cumulative process; 11 mean the spectrum peak search process; 12 mean the related operation process.

Fig. 2 is algorithm flow chart.

Fig. 3 is that the M section is with the isometric signal time domain schematic diagram of frequency; In figure: 13 mean that the M section is with the isometric signal of frequency, and 14-17 means that respectively the M section is with first to fourth segment signal in the isometric signal of frequency (13).

Fig. 4 is that the M section is with the not isometric signal time domain schematic diagram of frequency; In figure: 18 mean that the M section is with the not isometric signal of frequency, and 19-22 means that respectively the M section is with first to fourth segment signal in the not isometric signal of frequency (18).

Fig. 5 is the different isometric signal time domain schematic diagrames frequently of M section; In figure: 23 mean the different isometric signals frequently of M section, and 24-27 means respectively first to fourth segment signal in the different isometric signals (23) frequently of M section.

Fig. 6 is the different not isometric signal time domain schematic diagrames frequently of M section; In figure: 28 mean the different not isometric signals frequently of M section, and 29-32 means respectively first to fourth segment signal in the different not isometric signals (28) frequently of M section.

Fig. 7 is in the experiment parameter setting: signal to noise ratio (Signal Noise Ratio, SNR) SNR=-5dB, single hop signal length [N ₁, N ₂, N ₃, N ₄]=[50,50,50,50], all the other parameters arrange as in table 1 situation, the present invention and the isometric signal fused algorithm of multistage frequency reducing to 4 sections Bu Tong the isometric signal of frequency carry out the comparison diagram as a result of Frequency Estimation; In figure: 33 and 34 mean respectively the isometric signal fused algorithm of multistage frequency reducing and Frequency Estimation result of the present invention.

Fig. 8 for equal respectively-15dB of experiment parameter setting: SNR ,-13dB ,-11dB ,-9dB ,-7dB ,-5dB ,-3dB, 1dB, 3dB, 5dB, [N ₁, N ₂, N ₃, N ₄]=[50,50,50,50], all the other parameters arrange as in table 1 situation, the present invention and the isometric signal fused algorithm of multistage frequency reducing to 4 sections Bu Tong the isometric signal of frequency carry out root-mean-square error (the MSE:Mean Square Error) comparison diagram of Frequency Estimation; In figure: 35 and 36 mean respectively the isometric signal fused algorithm of multistage frequency reducing and Frequency Estimation root-mean-square error of the present invention.

Fig. 9 is at experiment parameter setting: SNR=-5dB, [N ₁, N ₂, N ₃, N ₄] equal respectively 4 sections of the table 2 different group of 1-10 shown in isometric signals frequently signal lengths, all the other parameters arrange as in table 1 situation, the present invention and the isometric signal fused algorithm of multistage frequency reducing to 4 sections Bu Tong the isometric signal of frequency carry out the root-mean-square error comparison diagram of Frequency Estimation; In figure: 37 and 38 mean respectively the isometric signal fused algorithm of multistage frequency reducing and Frequency Estimation root-mean-square error of the present invention.

Figure 10 is at experiment parameter setting: SNR=-5dB, [N ₁, N ₂, N ₃, N ₄]=[5,10,90,95], all the other ginsengs arrange as in table 1 situation, and the present invention carries out the figure as a result of Frequency Estimation to 4 sections different not isometric signals frequently; In figure: 39 mean Frequency Estimation result of the present invention, and 40 mean the actual value of signal frequency.

Figure 11 for equal respectively-15dB of experiment parameter setting: SNR ,-13dB ,-11dB ,-9dB ,-7dB ,-5dB ,-3dB, 1dB, 3dB, 5dB, [N ₁, N ₂, N ₃, N ₄]=[5,10,90,95], all the other parameters arrange as in table 1 situation, and the present invention carries out the root-mean-square error figure of Frequency Estimation to 4 sections different not isometric signals frequently.

Figure 12 is at experiment parameter setting: SNR=-5dB, [N ₁, N ₂, N ₃, N ₄] equaling respectively 4 sections of the table 2 different group of 1-10 shown in not isometric signals frequently signal lengths, all the other parameters arrange as in table 1 situation, and the present invention carries out the root-mean-square error figure of Frequency Estimation to 4 sections not isometric signals of different frequency.

Embodiment

The basic thought of a kind of multistage stationary signal domain related method of frequency of Frequency Estimation as shown in Figure 1.

Algorithm flow of the present invention as shown in Figure 2.

A kind of embodiment is as follows, with M=4, the Frequency Estimation object be take the not isometric signals of multistage different isometric signals frequently frequency different from multistage and is described as example, when estimate multistage with isometric signal or multistage frequently with frequently not during the frequency of isometric signal, be equivalent to estimate the different isometric signals frequently of multistage or the different not isometric signals frequently of multistage of at d (m)=0 o'clock.

1) rule of thumb, or by conventional frequency spectrum analysis methods such as DFT (Discrete Fourier Transform), FFT (Fast Fourier Transform) obtain frequency f to be estimated ₀roughly span f _scope, f _scope=[f _min, f _max]; By f _scopelinearity is divided into (A-1) part and (B-1) part, formation sequence f respectively _aand f _b, A and B are the integer that is greater than 1, f _a(a) mean sequence f _ain a element, a ∈ [1, A]; f _b(b) mean f _bin b element, b ∈ [1, B]; f _a(a _b) be illustrated in sequence f _ain with f _b(b) be worth immediate element, a _b∈ [1, A].

2) according to m segment signal x in M section stationary signal _mwith first segment signal x ₁between difference on the frequency d (m), m ∈ [1, M], according to following formula calculate in normalization matrix C the (m, the Elements C of a) locating (m, a),

C(m，a)＝f _A(a)+d(m)

3) calculate the normalization frequency spectrum according to following formula,

$X_m\left[f_A\left(a\right)\right]=0.5\underset{n_m=1}{\overset{N_m}{\mathrm{\&Sigma;}}}{x}_m\left(n_m\right)e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000510556000011.png"\; state="\{\{state\}\}">j</figure-callout>}{2\mathrm{\&pi;n}}_{m}C(m,a)/f_s}$

Wherein, X _m[f _a(a)] mean x _mthe normalization frequency spectrum; N _mfor positive integer, mean x _msampling number; n _mfor natural number, mean x _mtime series number; f _smean x _msample frequency, f _svalue meets nyquist sampling theorem.

4) calculate the cumulative frequency spectrum X[f of normalization frequency spectrum according to following formula _a(a)],

$X\left[f_A\left(a\right)\right]=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{abs}\left\{X_m\right[f_A\left(a\right)\left]\right\}$

In formula, abs means the mould of calculated complex;

5) design weighted factor e ^{-jD (m, a, b)}to step 3) in M section normalization frequency spectrum be weighted accumulation, obtain weighting accumulation spectral matrix X ' _b[f _a(a)], the D in weighted factor (m, a, b) is calculated by following formula,

D(m，a，b)＝θ _z(m，a _b)-θ _z(1，a _b)-g _b(a)[(2m-1)N _m+1]

X ' _b[f _a(a)] by following formula, calculated,

${X^{\&prime;}}_b\left[f_A\left(a\right)\right]=\mathrm{abs}\left\{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\right\{e^{-\mathrm{jD}(m,a,b)}{X}_m\left[f_A\left(a\right)\right]\left\}\right\}$

In formula, θ _z(m, a _b)=angle{X _m[f _a(a _b)]+W[f _a(a _b)], angle means the phase place of calculated complex; W[f _a(a _b)] mean that noise is at Frequency point f _a(a _b) the positive frequency part of frequency spectrum; g _b(a)=π [f _b(b)-f _a(a)]/f _s;

6) spectrum peak search weighting accumulation spectral matrix X ' _b[f _a(a)], its peak value element column is optimum weighting accumulation frequency spectrum

7) the frequency spectrum X[f that will add up _a(a)] with optimum weighting accumulation frequency spectrum

carry out related operation, obtain the spectrum correlation spectrum;

8) spectrum peak search spectrum correlation spectrum, calculate f as follows ₀estimated value

In formula, H _maxthe abscissa that means spectrum correlation spectrum peak element, K means the length of spectrum correlation spectrum.

According to present embodiment, the present invention is carried out to emulation, by the present invention with at present only a kind of can process the different frequency of multistage stationary signals method---the isometric signal fused algorithm of multistage frequency reducing compares, emulation experiment is divided into two parts, every partial simulation experiment comprises Monte Carlo Experiment 1000 times, in emulation experiment, institute's plus noise is additive white Gaussian noise, the first phase of arbitrary signal meets the Gaussian Profile that amplitude is 2 π, and all the other experiment parameters arrange as table 1.

Table 1 experiment parameter settings

Parameter name

f 0

fs

fmin

fmax

M

d(m)

A

B

Set point

10MHz

40MHz

9MHz

11.5 MHz

4

[0 17 23 29]KHz

120

120

Experiment is 1. partly for the different Frequency Estimation contrast experiments of isometric sinusoidal signal frequently of multistage; Because the isometric signal fused algorithm of multistage frequency reducing can not be processed the not isometric signal of multistage, therefore can only be analyzed with the present invention in experiment 1.;

Experiment 2. part for the different Frequency Estimation contrast experiments of isometric sinusoidal signal not frequently of multistage;

4 sections, table 2 not in homogenous frequency signal the single hop signal length table is set

From Fig. 8-Figure 12 simulation result, in the frequency of estimating the different isometric signals frequently of M section, estimated accuracy of the present invention improves a lot than the isometric signal fused algorithm of multistage frequency reducing; Estimating that the M section is different frequently not in the frequency of isometric signal, the present invention has higher precision.

Claims (1)

1. a kind of multistage stationary signal domain related method of frequency of Frequency Estimation, applicable object is M section stationary signal, and M is greater than 1 integer, and the length of every segment signal is any, frequency is constant, and the difference on the frequency of any two segment signals is known, and the method comprises the following steps:

1) structure normalization matrix C carries out respectively normalized to the frequency spectrum of M section stationary signal, obtains M section normalization frequency spectrum; In normalization matrix C the (m, the Elements C of a) locating (m, a) calculated by following formula,

C(m，a)＝f _A(a)+d(m)

In formula, m ∈ [1, M]; f _a(a) mean f _ain a element, a ∈ [1, A]; f _amean a sequence, A is greater than 1 integer, when A>2, and f _aby frequency f to be estimated ₀span f _scopelinear decile (A-1) part generates, when A=2, and f _a=f _scope; f _scope=[f _min, f _max], f _minmean f ₀the lower limit of value, f _maxmean f ₀the upper limit of value; D (m) means x _mwith x ₁between difference on the frequency, x _mmean m segment signal in M section stationary signal;

The normalization frequency spectrum is calculated by following formula,

$X_m\left[f_A\left(a\right)\right]=0.5\underset{n_m=1}{\overset{N_m}{\mathrm{\&Sigma;}}}{x}_m\left(n_m\right)e^{-j2\mathrm{\&pi;n}}{m}^C(m,a)/f_s$

In formula, X _m[f _a(a)], N _m, n _mand f _smean respectively x _mnormalization frequency spectrum, sampling number, time series number and sample frequency, N _mfor positive integer, n _mfor natural number, f _svalue meets nyquist sampling theorem;

It is characterized in that:

2) M section normalization frequency spectrum is added up, obtained cumulative frequency spectrum X[f _a(a)], X[f _a(a)] by following formula, calculated,

$X\left[f_A\left(a\right)\right]=\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{abs}\left\{X_m\right[f_A\left(a\right)\left]\right\}$

In formula, abs means the mould of calculated complex;

3) design weighted factor e ^{-jD (m, a, b)}to step 1) in M section normalization frequency spectrum be weighted accumulation, obtain weighting accumulation spectral matrix X ' _b[f _a(a)], D in weighted factor (m, a, b) is calculated by following formula,

D(m，a，b)＝θ _z(m，a _b)-θ _z(1，a _b)-g _b(a)[(2m-1)n _m+1]

X ' _b[f _a(a)] by following formula, calculated,

$X_b^{\&prime;}\left[f_A\left(a\right)\right]=\mathrm{abs}\left\{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\right\{e^{-\mathrm{jD}(m,a,b)}{X}_m\left[f_A\left(a\right)\right]\left\}\right\}$

In formula, a _b∈ [1, A], b ∈ [1, B], B is greater than 1 integer, θ _z(m, a _b)=angle{X _m[f _a(a _b)]+W[f _a(a _b)], angle means the phase place of calculated complex, f _a(a _b) be illustrated in sequence f _ain with f _b(b) be worth immediate element, f _b(b) mean sequence f _bin b element, when B>2, f _bby frequency f to be estimated ₀span f _scopelinear decile (B-1) part generates, when B=2, and f _b=f _scope; W[f _a(a _b)] mean that noise is at Frequency point f _a(a _b) the positive frequency part of frequency spectrum, g _b(a)=π [f _b(b)-f _a(a)]/f _s;

4) the weighting accumulation spectral matrix spectrum peak search step 3), its peak value element column is optimum weighting accumulation frequency spectrum;

5) by step 2) in cumulative frequency spectrum and step 4) in optimum weighting accumulation frequency spectrum carry out related operation, obtain spectrum correlation and compose;

6) spectrum peak search step 5) the frequency domain Correlated Spectroscopy, calculate f as follows ₀estimated value

$\hat{f}=(H_{\max }-1)\&times;(f_{\max }-f_{\min })/(K-1)+f_{\min }$

In formula, H _maxthe abscissa that means spectrum correlation spectrum peak element, K means the length of spectrum correlation spectrum.

Images