CN103441975B - A kind of Coded Signals parameter estimation method based on power spectrum - Google Patents

Abstract

The invention discloses a kind of Coded Signals parameter estimation method based on power spectrum, the method includes: obtain subpulse width τ corresponding during Coded Signals different code length and three dB bandwidth B_3dBFit correlation；Detect the power spectrum of Coded Signals to be valuated, described power spectrum is carried out amplitude normalization and obtains P (f), to described M point of P (f) uniform sampling, become sequence P (m) of a length of M, described P (m) is done discrete cosine transform and obtains Y_pre(P (m)), to described Y_pre(P (m)) carries out threshold process and obtains Y (P (m))；Seek the coenvelope of described Y (P (m)), calculate peak value number n of described coenvelope, obtain described signal code length N according to the relation of signal code length N Yu described peak value number n；Described Y (P (m)) is carried out inverse discrete cosine transform, obtains P_iM (), calculates described P_iThe carrier frequency f of (m)_c, three dB bandwidth B_3dBWith subpulse width τ.The present invention overcomes in prior art Coded Signals valuation result big and computationally intensive defect affected by noise under low signal-to-noise ratio environment, it is achieved the valuation of the Coded Signals parameter under the conditions of unknown prior information.

Description

A kind of Coded Signals parameter estimation method based on power spectrum

Technical field

The present invention relates to a kind of Coded Signals parameter estimation method based on power spectrum discrete cosine transform.

Background technology

Phase-coded signal, its phase modulation function is discrete limited state, belongs to discrete codes pulse compression letter Number, owing to phase code uses pseudo-random sequence, therefore it is also called pseudorandom encoded signal.Phase-coded signal is usually used in Doppler Change less occasion.Pseudorandom encoded signal can be classified by phase place value number, if phase place only takes two values, and referred to as two Encode signal mutually.

Coded Signals refers to that signal phase modulation function is discrete two values, and phase code typically uses pseudorandom Sequence, time-domain expression is:

$s\left(t\right)=v\left(t\right)\⊗\underset{k=0}{\overset{N-1}{\Σ}}{c}_k\mathrm{\δ}(t-k\mathrm{\τ})$

Wherein, τ is subpulse width, and N is code length, subpulse function

Coded Signals power spectrum expression formula is:

$\begin{array}{c}P\left(f\right)={\mathrm{\τ}}^2\sin c^2\[\mathrm{\τ}(f-f_c)\]\{{\[\underset{k=0}{\overset{N-1}{\Σ}}{c}_{k}cos\left(2\mathrm{\π}fk\mathrm{\τ}\right)\]}^2+{\[\underset{k=1}{\overset{N-1}{\Σ}}{c}_{k}sin\left(2\mathrm{\π}fk\mathrm{\τ}\right)\]}^2\}\\ ={\mathrm{\τ}}^2\sin c^2\[\mathrm{\τ}(f-f_c)\]\{N-1+\frac{\sin \[2N\mathrm{\π}\mathrm{\τ}(f-f_c)\]}{\sin \[2\mathrm{\π}\mathrm{\τ}(f-f_c)\]}\}\end{array}$

Coded Signals power spectrum feature is as follows:

1. power spectrum is symmetrical about mid frequency；

2. the peaks or valleys number in multimodal, code length and bandwidth is presented in bandwidth equal.

The signal more general due to the Parameter Estimation Problem of Coded Signals is complicated, therefore becomes the difficult point of research.Two The parameter encoding signal common mutually is as follows:

1. carrier frequency f_c；

2. three dB bandwidth B_3dB；

3. subpulse width τ, width time namely.

Existing Coded Signals parameter estimation method is all based on greatly the signal analysis of time domain, conventional parameter estimation side Method has time domain Phase difference, time frequency analysis and circulation spectrometry etc..In patent CN101984613A disclosed in 9 days March in 2011 Propose a kind of low rate BPSK burst bit rate method of estimation, utilize the method can obtain the code speed of Coded Signals Rate.2011, Cui Weiliang et al. was published in document " the Cyclic Spectrum estimation fast algorithm of improvement and the property of " electronics and information journal " Can analyze ".2012, Xu Huifa et al. was published in document " several exemplary phase based on FRFT coding letter of " war industry's journal " Number Detection and Parameter Estimation ".The amount of calculation of Coded Signals parameter estimation based on time frequency analysis and Cyclic Spectrum is the biggest, it is impossible to Solve the contradiction between estimated accuracy and operand, and be not suitable for low signal-to-noise ratio environment.

Summary of the invention

Present invention aim at: the valuation of the Coded Signals parameter under the conditions of the unknown prior information of realization, overcome existing There is in technology Coded Signals valuation result big and computationally intensive defect affected by noise under low signal-to-noise ratio environment.

The technical scheme is that and the invention provides a kind of biphase coding based on power spectrum discrete cosine transform letter Number parameter estimation method, it is characterised in that the method includes:

Step 1, obtain subpulse width τ corresponding during Coded Signals different code length N and three dB bandwidth B_3dBMatching Relation, it specifically includes:

Step a, set Coded Signals code length N span as { N₁,N₂,…,N_l,…,N_r, l ∈ 1,2 ..., r}, Set l=1；

Step b, set described Coded Signals code length as N_l, subpulse width τ constant interval is [τ₁,τ₂], change step A length of Δ τ；

Step c, in described τ constant interval, calculate the described Coded Signals power spectrum that different τ-value is corresponding, if institute Stating power spectrum maximum is P_max, search for first value of described power spectrum for P_maxFrequency f corresponding at/2₁With last value it is P_maxFrequency f corresponding at/2₂, according to definition B_3dB=f₂-f₁Calculate three dB bandwidth B that each described power spectrum is corresponding_3dB, then Obtain τ-B one to one_3dBData point；

Step d, selection exponential functionTo described τ-B_3dBData point carries out data matching；

Step e, obtain coefficient a, b that fitting function is respectively corresponding, obtain described subpulse width τ and described three dB bandwidth B_3dBFit correlation formula, wherein, fitting function refers to the fitting function that in step d, data matching obtains；

Step f, l is made to add 1；

Step g, it is judged that whether l is equal to r, if l is equal to r, performs step h, if l is not equal to r, returns step b；

Step h, draws described subpulse width τ the most corresponding for different code length N and described three dB bandwidth B_3dBMatching close It it is formula.

Step 2, detect the power spectrum of Coded Signals to be valuated, described power spectrum is carried out amplitude normalization and obtains P (f), to M point of described P (f) uniform sampling, becomes sequence P (m) of a length of M, described P (m) is done discrete cosine transform and obtains To Y_pre(P (m)), to described Y_pre(P (m)) carries out threshold process and obtains Y (P (m))；

Step 3, seek the coenvelope of described Y (P (m)), calculate peak value number n of described coenvelope, according to signal code length N with The relation of described peak value number n obtains described signal code length N；

Step 4, described Y (P (m)) is carried out inverse discrete cosine transform, obtain P_iM (), calculates described P_iThe carrier frequency f of (m)_c、 Three dB bandwidth B_3dBWith subpulse width τ, it specifically includes:

A, set described P_iM the maximum of () is P_max', search for described P_iM () first value is P_maxFrequency corresponding at '/2 f’₁It is P with last value_maxFrequency f corresponding at '/2₂', according to definitionCalculate described P_iThe load of (m) Frequently f_c；

B, according to definition B_3dB=f₂’-f₁' calculate described P_iThree dB bandwidth B of (m)_3dB；

C, according to described subpulse width τ and described three dB bandwidth B_3dBFit correlation formula obtain described subpulse width τ.

The invention has the beneficial effects as follows: the invention discloses a kind of biphase coding based on power spectrum discrete cosine transform letter Number parameter estimation method, the method includes: obtain subpulse width τ and 3dB band corresponding during Coded Signals different code length Wide B_3dBFit correlation；Detect the power spectrum of Coded Signals to be valuated, described power spectrum is carried out amplitude normalization and obtains To P (f), to M point of described P (f) uniform sampling, become sequence P (m) of a length of M, described P (m) is done discrete cosine transform Obtain Y_pre(P (m)), to described Y_pre(P (m)) carries out threshold process and obtains Y (P (m))；Seek the coenvelope of described Y (P (m)), meter Calculate peak value number n of described coenvelope, obtain described signal code length N according to the relation of signal code length N Yu described peak value number n； Described Y (P (m)) is carried out inverse discrete cosine transform, obtains P_iM (), calculates described P_iThe carrier frequency f of (m)_c, three dB bandwidth B_3dBAnd son Pulse width τ.The amount of calculation of the inventive method is mainly the complex multiplication amount of a FFT and twice real number dct transform Amount of calculation, overcomes the defect that the most conventional Time-Frequency Analysis Method is computationally intensive, under low signal-to-noise ratio environment, it is achieved unknown first Test the accurate valuation of the Coded Signals parameter of information condition.

Accompanying drawing explanation

Fig. 1 a kind of Coded Signals parameter estimation method flow diagram based on power spectrum disclosed by the invention；

Fig. 2 is disclosed by the invention determines Coded Signals subpulse width τ and three dB bandwidth B_3dBRelation flow chart；

Fig. 3 selected threshold disclosed by the invention flow chart；

Subpulse width τ and three dB bandwidth during Fig. 4 Coded Signals disclosed by the invention code length N={5,7,11,13} B_3dBThe matched curve of relation；

Wherein, when upper left is code length N=5；When upper right is code length N=7；When lower-left is code length N=11；Bottom right is code length N When=13；

Signal root-mean-square error curve chart after noise reduction during different threshold value in Fig. 5 embodiment of the present invention；

Wherein, when representing threshold value T=0.01 max (Y (m))；When * representing threshold value T=0.02 max (Y (m))；△ When representing threshold value T=0.03 max (Y (m))；Zero when representing threshold value T=0.04 max (Y (m))；Dotted line represents and does not take threshold value Time；

The root-mean-square error curve chart of signal after noise reduction during different parameters in Fig. 6 embodiment of the present invention；

Wherein,Representation parameter N=5, during the μ s of τ=0.01；Representation parameter N=5, during the μ s of τ=0.1；Zero representation parameter N When=7, τ=0.01 μ s；Representation parameter N=7, during the μ s of τ=0.1；Representation parameter N=11, during the μ s of τ=0.01；△ represents Parameter N=11, during the μ s of τ=0.1；* representation parameter N=13, during the μ s of τ=0.01；Representation parameter N=13, during the μ s of τ=0.1；

Coded Signals power spectrum signal graph after discrete cosine transform in Fig. 7 embodiment of the present invention；

Coded Signals power spectrum envelope diagram of signal after discrete cosine transform in Fig. 8 embodiment of the present invention；

The accuracy rate figure that in Fig. 9 embodiment of the present invention, Coded Signals code length is estimated；

Wherein, when * represents code length N=5；When △ represents code length N=7；When ◇ represents code length N=11；Zero represents code length N= When 13.

In Figure 10 embodiment of the present invention Coded Signals through experiment process before and after power spectrum chart with without noise power spectrum Figure；

Wherein, upper figure be not Noise time power spectrum chart；Middle figure is power spectrum chart during SNR=0dB；Figure below is upper figure Power spectrum chart after the inventive method processes.

The root-mean-square error figure of Coded Signals subpulse width valuation in Figure 11 embodiment of the present invention；

Wherein, when * represents pulsewidth τ=0.01 μ s；When representing τ=0.05 μ s；Zero when representing pulsewidth τ=0.1 μ s；

The root-mean-square error figure of Coded Signals carrier frequency valuation in Figure 12 embodiment of the present invention.

Wherein, when * represents pulsewidth τ=0.01 μ s；When representing τ=0.05 μ s；Zero when representing pulsewidth τ=0.1 μ s.

Detailed description of the invention

Illustrate hereinafter with reference to-12 pairs of embodiments of the present invention of Fig. 1.

As it is shown in figure 1, the embodiment of the present invention carries out under Coded Signals parameter estimation method based on power spectrum includes Row step:

Step 1, obtain subpulse width τ corresponding during Coded Signals different code length N and three dB bandwidth B_3dBMatching Relation, as in figure 2 it is shown, it specifically includes:

Step a, set Coded Signals code length N span as { N₁,N₂,…,N_l,…,N_r, l ∈ 1,2 ..., r}, Set l=1；

Step b, set described Coded Signals code length as N_l, subpulse width τ constant interval is [τ₁,τ₂], change step A length of Δ τ；

Step c, in described τ constant interval, calculate the described Coded Signals power spectrum that different τ-value is corresponding, if institute Stating power spectrum maximum is P_max, search for first value of described power spectrum for P_maxFrequency f corresponding at/2₁With last value it is P_maxFrequency f corresponding at/2₂, according to definition B_3dB=f₂-f₁Calculate three dB bandwidth B that each described power spectrum is corresponding_3dB, then Obtain τ-B one to one_3dBData point；

Step d, selection exponential functionTo described τ-B_3dBData point carries out data matching；

Step e, obtain coefficient a, b that fitting function is respectively corresponding, obtain described subpulse width τ and described three dB bandwidth B_3dBFit correlation formula；

Step f, l is made to add 1；

Step g, it is judged that whether l is equal to r, if l is equal to r, performs step h, if l is not equal to r, returns step b；

Step h, draws described subpulse width τ the most corresponding for different code length N and described three dB bandwidth B_3dBMatching close It it is formula.

Step 2, detect the power spectrum of Coded Signals to be valuated, described power spectrum is carried out amplitude normalization and obtains P (f), to M point of described P (f) uniform sampling, becomes sequence P (m) of a length of M, described P (m) is done discrete cosine transform and obtains To Y_pre(P (m)), to described Y_pre(P (m)) carries out threshold process and obtains Y (P (m))；

Wherein, in step 2, obtain described Y (P (m)) by following steps:

In Coded Signals power spectrum expression formula, making τ is subpulse width, and N is code length, f_cFor carrier frequency, two-phase is compiled Code power spectrum signal expression formula is:

$\begin{array}{c}P\left(f\right)={\mathrm{\τ}}^2\sin c^2\[\mathrm{\τ}(f-f_c)\]\{{\[\underset{k=0}{\overset{N-1}{\Σ}}{c}_{k}cos\left(2\mathrm{\π}fk\mathrm{\τ}\right)\]}^2+{\[\underset{k=1}{\overset{N-1}{\Σ}}{c}_{k}sin\left(2\mathrm{\π}fk\mathrm{\τ}\right)\]}^2\}\\ ={\mathrm{\τ}}^2\sin c^2\[\mathrm{\τ}(f-f_c)\]\{N-1+\frac{\sin \[2N\mathrm{\π}\mathrm{\τ}(f-f_c)\]}{\sin \[2\mathrm{\π}\mathrm{\τ}(f-f_c)\]}\}\end{array}$

To M point of P (f) uniform sampling, become sequence P (m) of a length of M, in wherein said discrete cosine transform formula

Described discrete cosine transform formula is:

And the signal actually received is that band is noisy, the signal after discrete cosine transform also can be by noise shadow Ring, therefore set threshold value T, by described Y_pre(P (m)) value is that the part zero setting of (-T, T) just obtains described Y (P (m)).

It should be noted that as it is shown on figure 3, the selecting step of threshold value is as follows:

A, described Y (P (m)) is done inverse discrete cosine transform obtain signal P_iM (), sets described P_i(m) and template signal P_nnM the root-mean-square error threshold value of () is σ_r；

B, set described threshold value T thick scope (0, Δ₁), described thick in the range of with 10^-nΔ₁For step-length, wherein n=1；

C, calculate described root-mean-square error σ respectively_i, i=1,2,3 ...；

D, select min{ σ_iThreshold value T corresponding to }_a；

E, judge described min{ σ_iWhether less than described σ_r, it is then to perform step h, no, then perform step f；

F, reset threshold range (T_a-10^-nΔ₁,T_a+10^-nΔ₁), step-length is 10^-(n+1)Δ₁；

G, make n add 1, return and perform step c；

H, choose T_aFor optimal threshold.

Wherein, in step a, P described in patent of the present invention_iM () is calculated by inverse discrete cosine transform, its Expression formula is as follows:

$P_i\left(m\right)=\sqrt{\frac{2}{M}}c\left(k\right)\underset{k=0}{\overset{M-1}{\Σ}}Y\left(P\right(m\left)\right)cos\frac{(2m+1)k\mathrm{\π}}{M},k,m=0,1,...,M-1;$

It should be noted that template signal P described in patent of the present invention_nnM () is defined as signal during not Noise；

It should be noted that P described in patent of the present invention_i(m) and described template signal P_nnThe root-mean-square error definition of (m) As follows:

The most described P of root-mean-square error_i(m) and described template signal P_nnThe square root of the variance between (m), its definition For:

$RMSE=\sqrt{\underset{m=0}{\overset{M-1}{\Σ}}{\[P_{nn}\left(m\right)-P_i\left(m\right)\]}^2/M};$

Step 3, seek the coenvelope of described Y (P (m)), calculate peak value number n of described coenvelope, according to signal code length N with The relation of described peak value number n obtains described signal code length N；

Wherein, in step 3, signal code length N with the relation of described Y (P (m)) coenvelope peak value number n is:

N=2n+1

It should be noted that coenvelope is defined as follows in the present invention:

Peak point in described Y (P (m)) is linked to be a curve, is the coenvelope of described Y (P (m))；

It should be noted that the peak value number determination methods of coenvelope is as follows:

Each point of described coenvelope is judged, when estimative point be simultaneously greater than on the left of it 5 continuous print points and On the right side of it during 5 continuous print points, just can determine that this point is peak point；

Step 4, described Y (P (m)) is carried out inverse discrete cosine transform, obtain P_iM (), calculates described P_iThe carrier frequency f of (m)_c、 Three dB bandwidth B_3dBWith subpulse width τ, it specifically includes:

Step 4, described Y (P (m)) is carried out inverse discrete cosine transform, obtain P_iM (), calculates described P_iThe carrier frequency f of (m)_c、 Three dB bandwidth B_3dBWith subpulse width τ, it specifically includes:

A, set described P_iM the maximum of () is P_max', search for described P_iM () first value is P_maxFrequency corresponding at '/2 f’₁It is P with last value_maxFrequency f corresponding at '/2₂', according to definitionCalculate described P_iThe load of (m) Frequently f_c；

B, according to definition B_3dB=f₂’-f₁' calculate described P_iThree dB bandwidth B of (m)_3dB；

C, according to described subpulse width τ and described three dB bandwidth B_3dBFit correlation formula obtain described subpulse width τ.

For embodiment, the present invention is described in detail with Coded Signals for the present invention；The specific embodiment of the invention contracts Write BPSK and represent Coded Signals, represent discrete cosine transform with abbreviation DCT, represent discrete cosine inversion with abbreviation IDCT Change, represent signal to noise ratio with abbreviation SNR.Simulation parameter is arranged: code length N={5,7,11,13}, subpulse width τ=0.05 μ s, carries Frequency is 500MHz.

The embodiment of the present invention may apply to military electronic antagonism and civilian frequency spectrum supervision field, at low signal-to-noise ratio environment Under estimate signal parameter exactly, can reflect function and the purposes of radar, improve the positioning precision of sonar tracking system with Improve tracking effect significant.

One, the τ-B that different code length is corresponding is calculated_3dBFit correlation formula；

As shown in Figure 4, and the goodness of fit is all higher than 99.5% to fit correlation curve, and fitting formula is:

1. 5 bpsk signal τ-B_3dBFitting formula:

τ=0.18604 B_3dB ^-0.9922

2. 7 bpsk signal τ-B_3dBFitting formula:

τ=0.82959 B_3dB ^-0.99148

3. 11 bpsk signal τ-B_3dBFitting formula:

τ=0.86469 B_3dB ^-0.98822

4. 13 bpsk signal τ-B_3dBFitting formula:

τ=0.21173 B_3dB ^-0.98466

Two, the choosing of threshold value T；

If the signal maximum after DCT is that (P (m}), sets described σ to max{Y_r=0.1, Δ₁=0.1 max{Y (P (m)) }, according to threshold value selection principle, when threshold value is taken as 0.02 max{Y (P (m)) time can meet threshold value and choose condition, if Determining bpsk signal parameter is code length N=13, SNR when-5dB～20dB change, described root-mean-square error curve as it is shown in figure 5, When Parameters variation, described root-mean-square error curve is as shown in Figure 6.

From Fig. 5-6, signal can be made to have noise reduction when taking threshold process, and when threshold value is taken as 0.02 max{Y (P (m)) } time, described root-mean-square error is respectively less than 0.1, and patent the most of the present invention takes threshold value T=0.02 max{Y(P(m))}。

Three, detect bpsk signal power spectrum P (f) to be valuated, and described bpsk signal is carried out code length estimation；

Described bpsk signal power spectrum is carried out dct transform, and signal after conversion is as it is shown in fig. 7, and do threshold process and obtain Y (P (m)), seeks the described coenvelope of signal after threshold process, as shown in Figure 8, obtains peak value number n, then obtains described The estimation of code length N of bpsk signal, when SNR excursion is [-5dB～20dB], carries out 100 Monte Carlo tests.

As it is shown in figure 9, when SNR is more than-5dB, bpsk signal code length estimates that accuracy rate is 100%, therefore can enter One step realizes the subpulse width to described bpsk signal and the accurate valuation of carrier frequency.

Four, described bpsk signal is carried out carrier frequency and the valuation of subpulse width；

Described Y (P (m)) is carried out idct transform, obtains P_iM (), sets SNR=0dB, as shown in Figure 10, through processing After signal closer to the signal of not Noise.

Bpsk signal when being [-5dB, 20dB] to SNR in interval carries out 100 MonteCarlo tests, calculates described The root-mean-square error (RMSE) of bpsk signal carrier frequency and subpulse width, result is as is illustrated by figs. 11 and 12.

The inventive method when SNR is [-5dB, 20dB] estimation difference of carrier frequency as shown in figure 11, and at SNR >-5dB Time carrier frequency estimation difference be respectively less than 0.5MHz.

The inventive method when SNR is [-5dB, 20dB] estimation difference of subpulse width as shown in figure 12, Er Qie During SNR >-5dB, the estimation difference of subpulse width is respectively less than 0.032 μ s.

Claims (1)

1. a Coded Signals parameter estimation method based on power spectrum, it is characterised in that the method includes:

Step 1, obtain subpulse width τ corresponding during Coded Signals different code length N and three dB bandwidth B_3dBFit correlation, It specifically includes:

Step a, set Coded Signals code length N span as { N₁,N₂,…,N_l,…,N_r, l ∈ 1,2 ..., and r}, set l =1；

Step b, set described Coded Signals code length as N_l, subpulse width τ constant interval is [τ₁,τ₂], change step be Δτ；

Step c, in described τ constant interval, calculate the described Coded Signals power spectrum that different τ-value is corresponding, if described two Coding power spectrum signal maximum is P mutually_max, search for first value of described Coded Signals power spectrum for P_maxAt/2 corresponding Frequency f₁It is P with last value_maxFrequency f corresponding at/2₂, according to definition B_3dB=f₂-f₁Calculate each described two-phase to compile Three dB bandwidth B that code power spectrum signal is corresponding_3dB, then τ-B one to one is obtained_3dBData point；

Step d, selection exponential functionTo described τ-B_3dBData point carries out data matching；

Step e, obtain coefficient a, b that fitting function is respectively corresponding, obtain described subpulse width τ and described three dB bandwidth B_3dB's Fit correlation formula, wherein, fitting function refers to the fitting function that in step d, data matching obtains；

Step f, l is made to add 1；

Step g, it is judged that whether l is equal to r, if l is equal to r, performs step h, if l is not equal to r, returns step b；

Step h, draws described subpulse width τ the most corresponding for different code length N and described three dB bandwidth B_3dBFit correlation formula；

Step 2, detect the power spectrum of Coded Signals to be valuated, described power spectrum is carried out amplitude normalization and obtains P F (), to described M point of P (f) uniform sampling, becomes sequence P (m) of a length of M, described P (m) is done discrete cosine transform and obtains To Y_pre(P (m)), to described Y_pre(P (m)) carries out threshold process and obtains Y (P (m))；

Step 3, seek the coenvelope of described Y (P (m)), calculate peak value number n of described coenvelope, according to signal code length N with described The relation of peak value number n obtains described signal code length N, wherein: N=2n+1；

Step 4, described Y (P (m)) is carried out inverse discrete cosine transform, obtain P_iM (), calculates described P_iThe carrier frequency f of (m)_c、3dB Bandwidth B_3dBWith subpulse width τ, it specifically includes:

A, set described P_iM the maximum of () is P_max', search for described P_iM () first value is P_maxFrequency f corresponding at '/2 '₁With Last value is P_maxFrequency f corresponding at '/2₂', according to definitionCalculate described P_iThe carrier frequency f of (m)_c；

B, according to definition B_3dB=f₂’-f₁' calculate described P_iThree dB bandwidth B of (m)_3dB；

C, according to described subpulse width τ and described three dB bandwidth B_3dBFit correlation formula obtain described subpulse width τ.