CN102215187A - 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 frequency domain correlation technique of Frequency Estimation

Technical field

The present invention relates to a kind of multistage stationary signal frequency domain correlation technique in signal processing field, particularly Frequency Estimation.

Background technology

The multistage stationary signal is meant the constant coherent signal of plurality of sections frequency 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 handled, after the low-pass filtering, can be constituted multistage with sinusoidal signal (being the identical sinusoidal signal of multistage frequency) frequently through the past.From the information theory viewpoint, 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 that improves the signal processing precision that the multistage signal is carried out fusion treatment, can be widely used in Frequency Estimation, LFM (Linear Frequency Modulated) signal parameter estimates that 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) the 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. phase association method of multisection sampled data [J]. 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 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; Though 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 handle the isometric signal of multistage, can not handle 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, needs to propose higher, the better frequency estimating methods of universality of a kind of precision.

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 the integer greater than 1; Secondly, by spectrum peak search weighting accumulation spectral matrix, obtain optimum weighting accumulation frequency spectrum; Then, the frequency spectrum that adds up that optimum weighting is accumulated frequency spectrum and M section normalization frequency spectrum carries out related operation, obtains the relevant spectrum of frequency domain; At last, the relevant spectrum of spectrum peak search frequency domain obtains high-precision frequency estimation.

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

1, a kind of multistage stationary signal frequency domain correlation technique of Frequency Estimation, applicable object is a M section stationary signal, and M is the integer greater than 1, 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 this method may further comprise the steps:

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

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

In the formula, m ∈ [1, M], f _ARepresent a sequence, by frequency f to be estimated ₀Span f _ScopeLinear five equilibrium (A-1) part generation, A is the integer greater than 1, f _Scope=[f _Min, f _Max], f _A(a) expression f _AIn a element, a ∈ [1, A]; D (m) represents x _mWith x ₁Between difference on the frequency, x _mM segment signal in the expression 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,HSA00000510556000042.png"\; state="\{\{state\}\}">j</figure-callout>}{2\mathrm{\&pi;n}}_{m}C(m,a)/f_s}$

In the formula, X _m[f _A(a)], N _m, n _mAnd f _sRepresent x respectively _mNormalization frequency spectrum, sampling number, time series number and sample frequency, N _mBe positive integer, n _mBe natural number, f _sValue satisfies nyquist sampling theorem;

It is characterized in that:

2) M section normalization frequency spectrum is added up, frequency spectrum X[f obtains adding up _A(a)], X[f _A(a)] calculate by following formula,

$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 the formula, abs represents the mould of calculated complex;

3) design weighted factor e ^{-jD (m, a, b)}M section normalization frequency spectrum in the step 1) is weighted accumulation, obtains weighting accumulation spectral matrix X ' _b[f _A(a)], D in the weighted factor (m, a b) are 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)] calculate by following formula

${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 the formula, a _b∈ [1, A], b ∈ [1, B], B is the integer greater than 1, θ _z(m, a _b)=angle{X _m[f _A(a _b)]+W[f _A(a _b)], angle represents 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) expression sequence f _BIn b element, f _BBy frequency f to be estimated ₀Span f _ScopeLinear five equilibrium (B-1) part generation, W[f _A(a _b)] represent 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) accumulation of the weighting in spectrum peak search step 3) spectral matrix, its peak value element column are optimum weighting accumulation frequency spectrum;

5) with step 2) in add up that optimum weighting accumulation frequency spectrum carries out related operation in frequency spectrum and the step 4), obtain the relevant spectrum of frequency domain;

6) the relevant spectrum of spectrum peak search step 5) frequency domain is calculated f as follows ₀Estimated value

In the formula, H _MaxThe abscissa of the relevant spectrum peak element of expression frequency domain, K represents the length of the relevant spectrum of frequency domain.

Frequency estimating methods precision height, the universality that the present invention relates to are good, are applicable to multistage length stationary signal arbitrarily.

Description of drawings

The present invention is further elaborated with concrete enforcement with reference to the accompanying drawings below.With M=4, signal type is that example describes with the sinusoidal signal.

Fig. 1 is the algorithm basic thought.

Among the figure: 1 expression M section sinusoidal signal comprises that mainly the M section is with the isometric signal of frequency, M the section not isometric signal of frequency together, the different isometric signals frequently of M section, the not isometric signals of M section difference frequency; 2 expression M section normalization frequency spectrums; 3 represent to add up frequency spectrum, 4 expression weightings accumulation spectral matrix; The optimum weighting accumulation of 5 expressions spectral matrix; The relevant spectrum of 6 expression frequency domains; 7 expression frequency estimations; 8 expression normalization processes; 9 expression cumulative process; 10 expression weighting cumulative process; 11 expression spectrum peak search processes; 12 expression related operation processes.

Fig. 2 is an algorithm flow chart.

Fig. 3 is that the M section is with the isometric signal time domain schematic diagram of frequency; Among the figure: 13 expression M sections are with the isometric signal of frequency, and 14-17 represents 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; Among the figure: 18 expression M sections are with the not isometric signal of frequency, and 19-22 represents 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; Among the figure: the different isometric signals frequently of 23 expression M sections, 24-27 represents first to fourth segment signal in the different isometric signals (23) frequently of M section respectively.

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

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 are provided with as under table 1 situation, and the present invention and the isometric signal fused algorithm of multistage frequency reducing carry out the comparison diagram as a result of Frequency Estimation to 4 sections isometric signals of different frequencies; Among the figure: 33 and 34 represent isometric signal fused algorithm of multistage frequency reducing and Frequency Estimation result of the present invention respectively.

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

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

Figure 10 is at experiment parameter setting: SNR=-5dB, [N ₁, N ₂, N ₃, N ₄]=[5,10,90,95], all the other ginsengs are provided with as under 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; Among the figure: 39 expressions Frequency Estimation result of the present invention, the actual value of 40 expression signal frequencies.

Figure 11 for to equal respectively at experiment parameter setting: SNR-15dB ,-13dB ,-11dB ,-9dB ,-7dB ,-5dB ,-3dB, 1dB, 3dB, 5dB, [N ₁, N ₂, N ₃, N ₄]=[5,10,90,95], all the other parameters are provided with as under 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 4 sections of the table 2 different group of 1-10 shown in the not isometric signals frequently signal lengths respectively, all the other parameters are provided with as under 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 frequencies.

Embodiment

The basic thought of a kind of multistage stationary signal frequency domain correlation technique 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 is that example describes with the different not isometric signals frequently with multistage of the different isometric signals frequently of multistage, when estimate multistage with isometric signal or multistage frequently with frequently not during the frequency of isometric signal, be equivalent to estimate different isometric signals frequently of multistage or the different not isometric signals frequently of multistage at d (m)=0 o'clock.

1) rule of thumb, or by conventional frequency spectrum analysis method such as DFT (Discrete Fourier Transform), FFT (Fast Fourier Transform) obtains frequency f to be estimated ₀Roughly span f _Scope, f _Scope=[f _Min, f _Max]; With 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 greater than 1, f _A(a) expression sequence f _AIn a element, a ∈ [1, A]; f _B(b) expression 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 the M section stationary signal _mWith first segment signal x ₁Between difference on the frequency d (m), m ∈ [1, M], according to following formula calculate among the 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,HSA00000510556000042.png"\; state="\{\{state\}\}">j</figure-callout>}{2\mathrm{\&pi;n}}_{m}C(m,a)/f_s}$

Wherein, X _m[f _A(a)] expression x _mThe normalization frequency spectrum; N _mBe positive integer, expression x _mSampling number; n _mBe natural number, expression x _mTime series number; f _sExpression x _mSample frequency, f _sValue satisfies nyquist sampling theorem.

4) calculate the frequency spectrum X[f that adds up 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 the formula, abs represents the mould of calculated complex;

5) design weighted factor e ^{-jD (m, a, b)}M section normalization frequency spectrum in the step 3) is weighted accumulation, obtains weighting accumulation spectral matrix X ' _b[f _A(a)], the D in the weighted factor (m, a b) are 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)] calculate by following formula,

${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 the formula, θ _z(m, a _b)=angle{X _m[f _A(a _b)]+W[f _A(a _b)], angle represents the phase place of calculated complex; W[f _A(a _b)] represent 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 relevant spectrum of frequency domain;

8) the relevant spectrum of spectrum peak search frequency domain is calculated f as follows ₀Estimated value

In the formula, H _MaxThe abscissa of the relevant spectrum peak element of expression frequency domain, K represents the length of the relevant spectrum of frequency domain.

According to present embodiment the present invention is carried out emulation, with the present invention with at present only a kind of can handle 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 1000 Monte Carlo experiments, institute's plus noise is additive white Gaussian noise in the emulation experiment, the first phase of arbitrary signal satisfies the Gaussian Profile that amplitude is 2 π, and all the other experiment parameters are provided with as table 1.

The table 1 experiment parameter value of setting

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 at 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 handled the not isometric signal of multistage, so can only in experiment 1., be analyzed with the present invention;

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

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

By Fig. 8-Figure 12 simulation result as can be known, 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 frequency domain correlation technique of Frequency Estimation, applicable object is a M section stationary signal, and M is the integer greater than 1, 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 this method may further comprise the steps:

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

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

In the formula, m ∈ [1, M], f _ARepresent a sequence, by frequency f to be estimated ₀Span f _ScopeLinear five equilibrium (A-1) part generation, A is the integer greater than 1, f _Scope=[f _Min, f _Max], f _A(a) expression f _AIn a element, a ∈ [1, A]; D (m) represents x _mWith x ₁Between difference on the frequency, x _mM segment signal in the expression 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^{-j{2\mathrm{\&pi;n}}_{m}C(m,a)/f_s}$

In the formula, X _m[f _A(a)], N _m, n _mAnd f _sRepresent x respectively _mNormalization frequency spectrum, sampling number, time series number and sample frequency, N _mBe positive integer, n _mBe natural number, f _sValue satisfies nyquist sampling theorem;

It is characterized in that:

2) M section normalization frequency spectrum is added up, frequency spectrum X[f obtains adding up _A(a)], X[f _A(a)] calculate by following formula,

$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 the formula, abs represents the mould of calculated complex;

3) design weighted factor e ^{-jD (m, a, b)}M section normalization frequency spectrum in the step 1) is weighted accumulation, obtains weighting accumulation spectral matrix X ' _b[f _A(a)], D in the weighted factor (m, a b) are 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)] calculate by following formula,

${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 the formula, a _b∈ [1, A], b ∈ [1, B], B is the integer greater than 1, θ _z(m, a _b)=angle{X _m[f _A(a _b)]+W[f _A(a _b)], angle represents 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) expression sequence f _BIn b element, f _BBy frequency f to be estimated ₀Span f _ScopeLinear five equilibrium (B-1) part generation, W[f _A(a _b)] represent 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) accumulation of the weighting in spectrum peak search step 3) spectral matrix, its peak value element column are optimum weighting accumulation frequency spectrum;

5) with step 2) in add up that optimum weighting accumulation frequency spectrum carries out related operation in frequency spectrum and the step 4), obtain the relevant spectrum of frequency domain;

6) the relevant spectrum of spectrum peak search step 5) frequency domain is calculated f as follows ₀Estimated value

In the formula, H _MaxThe abscissa of the relevant spectrum peak element of expression frequency domain, K represents the length of the relevant spectrum of frequency domain.

Images