Roll angular rate extraction based on modified spline-kernelled chirplet transform

Hui Zhao ,Zhong Su ,*,Qing Li ,Fu-chao Liu ,Ning Liu

aSchool of Automation,Beijing Institute of Technology,Beijing,100081,China

bBeijing Key Laboratory of High Dynamic Navigation Technology,Beijing Information Science and Technology University,Beijing,100192,China

Keywords:High spinning projectile Roll angular rate Time-frequency analysis Spline-kernelled chirplet transform Chirp Z-transform

ABSTRACT The roll angular rate is much crucial for the guidance and control of the projectile.Yet the high-speed rotation of the projectile brings severe challenges to the direct measurement of the roll angular rate.Nevertheless,the radial magnetometer signal is modulated by the high-speed rotation,thus the roll angular rate can be achieved by extracting the instantaneous frequency of the radial magnetometer signal.The objective of this study is to find out a precise instantaneous frequency extraction method to obtain an accurate roll angular rate.To reach this goal,a modified spline-kernelled chirplet transform(MSCT)algorithm is proposed in this paper.Due to the nonlinear frequency modulation characteristics of the radial magnetometer signal,the existing time-frequency analysis methods in literature cannot obtain an excellent energy concentration in the time-frequency plane,thereby leading to a terrible instantaneous frequency extraction accuracy.However,the MSCT can overcome the problem of bad energy concentration by replacing the short-time Fourier transform operator with the Chirp Z-transform operator based on the original spline-kernelled chirplet transform.The introduction of Chirp Z-transform can improve the construction accuracy of the transform kernel.Since the construction accuracy of the transform kernel determines the concentration of time-frequency distribution,the MSCT can obtain a more precise instantaneous frequency.The performance of the MSCT was evaluated by a series of numerical simulations,high-speed turntable experiments,and real flight tests.The evaluation results show that the MSCT has an excellent ability to process the nonlinear frequency modulation signal,and can accurately extract the roll angular rate for the high spinning projectiles.

1.Introduction

Guidance reformation of the conventional projectile is an effective way to improve the hit accuracy and reduce the collateral damage[1].However,the high-speed rotation of the projectile poses a profound challenge to achieve accurate navigation and guidance during the flight[2,3].For example,the howitzer's projectiles can reach a maximum rotation speed of 300 r/s[4].Under such a high rotation speed condition,the commonly used navigation methods will operate abnormally.Specifically,the Global Navigation Satellite System(GNSS)may encounter loss of signal[5].The inertial navigation system(INS),due to the saturation of the gyroscope range,will completely fail[6].Thus,the INS must wait until the rotation speed is reduced into the gyroscope range before starting to power up.However,it will bring great difficulties to the initial alignment of the INS.Moreover,the roll angle and the roll angular rate are essential in terms of the performance of the guidance and control.To solve the above-mentioned problems introduced by the high-speed rotation,many schemes are designed to estimate the roll angle and the roll angular rate of the spinning projectiles with low-cost sensors.

According to the principle of centrifugal force,multiple accelerometers were employed to calculate the roll angular rate in Ref.[7].However,this scheme is not feasible due to the large volume and the influence of external force.Additionally,the solar azimuth sensor[8]and the infrared sensor[9]were also used to measure the roll angular rate,but they might suffer from the environment and the weather.

To improve the estimation performance of the roll angular rate,multi-source fusion filters were designed in many studies.Christophe et al.[10] proposed using a 6-degree-of-freedom ballistic model, combined with a three-axis magnetometer and a two-axis accelerometer, to estimate the roll angular rate.Raul et al.[11]used accelerometers,GNSS and,semi-active laser quadrant photodetector, also assisted by the ballistic models to obtain the gravity vectors, the velocity vectors and, the line-of-sight vectors, thereby solving the roll angle of the projectile by multi-vector attitude determination algorithm.Long et al.[12]also adopted a three-axis magnetometer,a two-axis gyroscope and,a GNSS to determine the attitude of the projectile.However, the multi-sensor integration scheme counteracts the advantages of low-cost and small-size sensors.Moreover, the parameter setting of the multi-source fusion filter is also challenging work.

Many pieces of literature treated the magnetometer as an effective alternative to perceive the attitude changes of the projectile.Xiang et al.[13]presented an attitude measurement method for high spinning projectiles using only a three-axis magnetometer.According to their description, the roll angle measurement accuracy can reach 0.5.However, magnetometers are susceptible to external ferromagnetic materials and strong current fields and also need a complex pre-calibration process.Li et al.[14] designed a three-axis orthogonal magnetic induction coil, combined with the zero-phase threshold processing technique, to measure the roll angular rate.Similar to Li's plan,Shang et al.[15]extracted the roll angular rate by estimating the instantaneous frequency (IF) of the radial axis magnetometer signal.Besides, Deng et al.[16] adopted the phase-locked loop (PLL) technique to track the phase angle of the radial magnetometer signal,thereby obtaining the roll angle of the projectile.

The objective of this paper is to propose a roll angular rate extraction method for the high spinning projectiles,and also solve the problem of navigation and guidance failure caused by gyroscope range saturation.First of all,the high spin we define refers to the situation where the rotation speed is greater than 20 r/s,that is,7200/s.The rotation speed of the spinning projectiles may be distributed between a few revolutions per second to a few hundred revolutions per second.However, if the rotation speed is greater than 20 r/s, the range of the high-precision gyroscopes will be saturated.Even if there are some large-range gyroscopes, their measurement accuracy is poor and can only make a rough estimate of the rotation speed.Therefore, we adopted the rotation speed of 20 r/s to define the high spin.

According to the literature, the radial magnetometer scheme can realize an accurate estimation of the roll angular rate in a harsh environment.However, the key to the scheme is a high-precision instantaneous frequency (IF) estimation method.Traditional IF estimation methods include short-time Fourier transform (STFT)[17], continuous wavelet transform (CWT) [18], Wigner-Ville distribution(WVD)[19],and etc.Among them,the STFT and the CWT are typical linear transforms, which may be limited by the Heisenberg uncertainty principle.For nonlinear frequency modulated signals, the linear transforms cannot accurately describe the local time-frequency distribution [20].Although the WVD is a typical quadratic term transform, it suffers from the cross-term interference introduced by the multicomponent signals [21].Researchers had put forward many improved forms of the WVD, for example,smoothed pseudo-Wigner-Ville distribution (SPWVD), but the time-frequency resolution of the SPWVD is worse than the original WVD.

Recently, Peng et al.[22] presented a class of parameterized time-frequency analysis (PTFA) methods, including polynomial chirplet transform (PCT) [23], spline-kernelled chirplet transform(SCT) [24], and generalized warblet transform (GWT) [25].The essence of the PTFA is describing the intrinsic characteristics of the signal by the constructed transform kernel.According to the literature, the PTFA can obtain more concentrated time-frequency distribution (TFD) than the traditional methods for nonlinear frequency modulated signals.However, our preliminary test results show that the IF estimation accuracy of the PTFA cannot meet the actual roll angular rate extraction requirement.Especially, 0.01 Hz IF estimation error will result in a roll angular rate error of 3.6/s,which will seriously deteriorate the accuracy of subsequent navigation updating.Through a deeply investigated PTFA,we found that its shortcoming lay in the construction of the transform kernel.

In this paper,we introduce the chirp Z-transform(CZT)into the initialization process of the SCT transform kernel, thereby enhancing the IF estimation accuracy.The proposed algorithm is named modified SCT (MSCT).It can provide a more concentrated time-frequency distribution (TFD) with nonlinear frequency modulated signals than that of the SCT algorithm.The contribution of this paper is presenting a high-precision IF estimation method,which can be used to solve the problem of accurate extraction of the roll angular rate for high spinning projectiles.Compared with the other methods provided by the related literature, the innovations of our method are as follows.(1) Minimal sensor resources are required.We only employed a three-axis magnetometer to estimate the roll angular rate.(2)Small-size,lowcost, and low power consumption hardware.(3) Our proposed method will not be affected by the weather.(4) No strict precalibration process is required, and the workflow is simple.(5)The roll angular rate estimation accuracy is higher than the traditional methods.

Based on the above content, this paper is organized as follows.The roll modulation of strap-down radial magnetometer signal is described in detail in Section 2.The principle of MSCT is derived in Section 3.Numerical validation is carried out in Section 4.In Section 5,the performance of the MSCT algorithm is further verified by the high-speed turntable experiment and the actual flight tests.The conclusion is provided in Section 6.

2.Strapdown radial magnetometer roll-modulated signal analysis

For the convenience of description,three coordinate systems are defined in the first, which are the local geographic coordinate system[26],the local geomagnetic coordinate system[27],and the projectile coordinate system.Among them, the local geographic coordinate system (North-East-Ground) is selected as the navigation reference coordinate system.Theaxis points to the North,theaxis points to the east,and thepoints vertically to the ground.The projectile coordinate system adopts the center of mass of the projectile as the origin point.Thepoints to the warhead along the axial direction of the projectile.Theaxis is along the radial direction of the projectile and perpendicular to theaxis.Theis along the radial direction of the projectile and perpendicular to theplane.To accurately measure the attitude of the projectile, the installation position of the three-axis magnetometer coincides with the projectile coordinate system, as shown in Fig.1.

The geomagnetic field,similar to the gravity field,is an intrinsic property of the earth.Usually,seven elements are used to represent the geomagnetic information at a certain point on the earth.The seven elements are magnetic declination(),magnetic inclination(), horizontal intensity (), north component (), east component (), vertical component (), and total field () [28],respectively.

The relationships between the seven elements[29]are shown in Eq.(1).

Fig.1.The coordinate systems.

According to our previous research work[6],we can obtain the ideal output of two radial magnetometers.Their specific expressions can be written as follows:

If we set

Eq.(2) can be rewritten as

where,ωis the roll angular rate of the projectile.

Since the projectile has a slower change in pitch and yaw during flight, then we can obtain

wheredenotes the rotation frequency.

Likewise, we can obtain

According to Eq.(7) and Eq.(8), if we can extract the instantaneous frequencyof the Y-axis and Z-axis magnetometer signal,then we will acquire the roll angular rate ωof the projectile.

3.The modified spline-kernelled chirplet transform (MSCT)

3.1.Spline-kernelled chirplet transform (SCT)

The parameterized time-frequency analysis method can provide a time-frequency distribution with excellent energy concentration for nonstationary signals by constructing the appropriate transform kernel.Typical parameterized time-frequency analysis method mainly includes linear chirplet transform (LCT) [30], PCT, SCT, and GWT.The main difference between them lies in the structure of the transform kernel.Especially, the SCT employs the cubic spline fitting model to construct the transform kernel, which has higher fitting accuracy for complex frequency curves.Besides,the SCT is an extension of the STFT, LCT, and PCT.More details on the timefrequency analysis method and the SCT can be found in Ref.[25].Next,we will briefly discuss the transformation principle of the SCT.

For a signal of（）∈（）, the spline-kernelled chirplet transform (SCT) is defined as

In Eq.(10),tis the-th node of the spline kernel and there exists-∞=<<…<<= +∞.=｛｝,(=1，…,=1，…-1),denotes the local polynomial coefficient matrix of the spline kernel.When= 0, the SCT degenerates into the STFT.When the total number of nodes is two, that is= 2, the SCT degenerates into the PCT.When=2 and= 2, the SCT degenerates into the LCT.

The basic principle of the SCT can be illustrated in the following steps:

3) Finally, the obtained signal is processed by the STFT with a Gaussian window function.

Under the ideal condition where the constructed kernel function can perfectly characterize the real IF law of the signal, the frequency resolution of the SCT will be a constant for the wholetime span.Finally, a simple STFT operator can be used to analysis the pseudo non-modulation signals and get a high concentration in time-frequency distribution (TFD).

1) The similarity between the spline kernel function and the real IF law of the signal;

2) The sampling frequency;

3) The time width of window function σ.

The measures to improve the SCT frequency resolution include finding the best spline kernel function, decreasing the sampling frequency, and increasing the time width of the window function.

According to Ref.[22],a better transform kernel can be obtained by repeating the execution of the SCT.In other words, the SCT should be implemented at least twice.For the first iteration,= 0,the SCT degenerates into the STFT.A course TFD can be obtained by the STFT.Then a rough IF law of the signal can be extracted from the course TFD.Finally, the spline fitting method is adopted to approximate the rough IF law and acquire the coefficient matrix of the spline kernel,i.e.,C.The constructed coefficient matrix C will be used for the second iteration.As mentioned above, if the constructed kernel function equals the real IF law of the signal (=), the SCT will get the best frequency resolution.

To summarize,the first iteration of the SCT is very important for the TFD and the frequency resolution.However,the first iteration of the SCT is simply an STFT operation.Therefore, we can choose another operation with a higher frequency resolution than the STFT to improve the overall performance of the SCT.In this paper, we select the Chirp Z-transform (CZT) with a Gaussian window function to replace the STFT.

3.2.Chirp Z-transform (CZT)

The Chirp-Z transform (CZT) is a special kind of Z-transformation.Its essence is the Z-transformation with an equal-angle sampling of a spiral trajectory on the Z-plane.The CZT is an efficient spectrum analysis method, which is usually used to calculate a specific scope of frequency.That is, the CZT can perform a spectral zoom analysis focusing on a smaller band of frequencies and break through the limitation of the STFT[31].

For example,（）is a signal sequence of length N,（）is the Ztransformation of（）.By using the CZT, the Z-transformation at pointis

At the same time, we have

whereis the radius of the initial sample point.θdenotes the starting angular frequency which is defined as θ= 2π（/）.A is used to determines the position of initial sampling point of the transform.is the stretch rate of the spiral trajectory.φis the angular frequency difference between the adjacent sampling points and can be expressed as φ= 2π（-）/（-1）.［，］is the scope of refined frequency band.andrepresent the minimum frequency and the maximum frequency, respectively.M is the number of sampling points.is the sampling frequency.

When the values of A, W, and M are set, a specific refined frequency band will be determined.If=1,the sampling trajectory will be an arc with radius of.At this moment, if= 1, the sampling trajectory will be an arc on the unit circle.The initial sampling point and the stop sampling point of this arc depend on θand φ.

It is well known that the FFT of the signal is equivalent to the Z transformation with an equal-interval sampling on the unit circle.When=1,=1,the CZT will be equivalent to the FFT on the specified frequency band.However,the frequency resolution of the CZT is higher than that of the FFT.Therefore,we can adopt the CZT to replace the FFT in the SCT to obtain a higher frequency resolution and improve the overall performance of the SCT.

Compared with the FFT, the CZT has the following two outstanding advantages.

1) The angular frequency difference between the adjacent sampling points, φ, can be set arbitrarily.Therefore, the frequency resolution of the CZT can also be arbitrary in theory.

2) The initial sampling pointcan be selected arbitrarily, that is,the CZT can analyze a signal using any initial frequency.

3.3.The modified spline-kernelled chirplet transform (MSCT)

Considering that the interested frequencies of the signal are distributed in a narrow band, the modified spline-kernelled chirplet transform (MSCT) is proposed in this paper.The essence of the MSCT is to replace the FFT operator in the SCT with the CZT operator.The CZT's excellent frequency band subdivision ability endows the MSCT with a better TFD concentration and a frequency resolution.The principle and execution process of the MSCT is described in detail below.

According to Eq.(9), we can get the discrete form of the SCT.

The discrete form of the STFT is

The discrete form of the CZT is

By comparing Eq.(13) with Eq.(14), we can find that the SCT mainly contains three operators.They are frequency-rotate operator, frequency-shift operator, and STFT operator, respectively.If replacing the STFT operator with the CZT operator in the SCT, we will obtain the expression of the MSCT, its discrete form is as follows.

The implementation steps of the MSCT are identical to those of the SCT.However, the SCT usually requires multiple iterations to find the best spline kernel function.The proposed MSCT only requires two iterations,and the first iteration is used to initialize the parameters of the spline transform kernel.In the second iteration,the MSCT can achieve a more concentrated TFD.The detailed implementation steps of the MSCT are as follows.

For a signal sequence（）with a length of N.

1) Perform the Hilbert transform on（） and obtain（）.

2) Use the Gaussian window functionσ（）with a time-width of σ to intercept signal（）.At the same time,zero padding is carried out on both sides of the signal（）.Then we can obtain

3) Process the intercepted signal fragments by the CZT.The initial frequency, the stop frequency, and the number of sampling points M are set based on the priori information of the signal and the accuracy requirement.Then we can get a rough TFD of the signal.

4) Extract the ridge curve from the rough TFD by finding the peak at every moment.The ridge curve is close to the IF law of the signal.

whereis the number of the b-splines and cis the control node.All the control nodes form a knot sequence.（）denotes the j-th n-order b-spline for a knot sequence, which can be computed as

Rewrite B-spline into a piecewise polynomial form,we can get

where,（） is the coefficients of the spline transform kernel.

7) Utilize the frequency rotation operator to process（） and obtain the rotated signal（）.

8) Utilize the frequency-shift operator to process（）and obtain the shifted signal（）.

9) Repeat steps (2)-(4) on the signal（）, then the refined TFD and a more accurate IF can be obtained.

For the convenience of understanding, we provide a flow chart of the MSCT, as shown in Fig.2.

To further demonstrate the advancement of the MSCT, we analyzed the frequency resolution of three algorithms, including the STFT, the SCT, and the MSCT.

If the signal in the sliding window is an ideal stationary signal,the frequency resolution of the STFT is

where,fis the sampling frequency,N is the length of the window.

If the constructed spline kernel model exactly matches with the real IF law of the signal, the frequency resolution of the SCT and MSCT are

According to Eqs.(27) and (28), under ideal conditions, the frequency resolution of the STFT and the SCT is depended on the sampling frequency and the window length.However, those two parameters cannot be set arbitrarily due to the limitation of practical conditions.As shown in Eq.(29), the frequency resolution of the MSCT is determined by the initial frequency, stop frequency,sampling frequency, and the number of sampling points in the narrow band.Among them, the frequency range can be set as narrow as possible, the number of sampling points in the narrow band can be set arbitrarily.In this way, the MSCT has an arbitrary frequency resolution under the ideal conditions.Therefore, the MSCT is very suitable for IF estimation for the high spinning projectiles.

Fig.2.The flow chart of the MSCT.

4.Numerical studies

To evaluate the performance of the MSCT on IF estimation, a simulated nonstationary signal in Ref.[24], is used, which is given by

The instantaneous frequency of the above signal is

The instantaneous frequency (IF) is shown in Fig.3.

The generated signal is artificially added with additive Gaussian noise with a variance of 0.64 and a mean of zero.The SNR(signalnoise ratio)of the signal is 1.7335 dB.The sampling frequency is set to be 100 Hz, the contaminated signal and the original signal are shown in Fig.4.

The MSCT is used to analyze the signal's time-frequency distribution and extract its instantaneous frequency.Specifically, the sampling frequency is set to be 100 Hz, the window length is 512,the number of splines is 25,and the initial parameter matrix of the SCT is 0.Since the instantaneous frequency of the above curve is distributed between 4 Hz and 16 Hz, the start frequency of the narrowband is set to 0 Hz, the stop frequency is 20 Hz, and the number of spectrum subdivision points is 2000.

The TFD obtained by the MSCT is compared with the TFDs obtained by the Short-time Fourier transform (STFT), the Smoothed Pseudo Wigner-Ville Distribution (SPWVD), and the original SCT.The window length of the STFT is 512.The window sliding step size is 1.The time smoothing window size of SPWVD is 1000, the frequency smoothing window size is 512.The selected window function is Kaiser windows with shape factor β = 20.The window length of the SCT is 512, the number of splines is 25.For the iteration number of the SCT and MSCT, both of them are 2.The first iteration is to initial the parameters of the transform kernel.The second iteration is to generate a more concentrated TFD.The TFDs generated by the STFT,SPWVD,SCT,and MSCT are shown in Fig.5.

As shown in Fig.5(a), it is difficult for the STFT to identify the underlying time-frequency pattern of the nonlinear frequency modulated signal.But it can characterize the stable part.The TFD generated by the SPWVD,as shown in Fig.5(b),shows poor time or frequency resolution because of the interferes of multicomponent cross terms.At the same time, the higher the noise intensity, the more serious the interference.However, in Fig.5(c) and (d), both the SCT and the MSCT produce the TFDs with an excellent concentration,it is difficult to distinguish the difference between these two TFDs.Then we can acquire the instantaneous frequency by extracting the ridge curve of the TFD.The ridge curve is extracted by finding the peak point of the TFD.It is defined as

Fig.3.The instantaneous frequency of the simulation signal.

Fig.4.The contaminated simulation signal and the original signal.

According to the four TFDs shown in Fig.5,We can get four sets of estimates of instantaneous frequencies.The estimated instantaneous frequency and its corresponding estimation errors are described in Fig.6 and Fig.7.

As shown in Figs.6 and 7,the STFT can accurately extract the IF when the frequency is stable, but it can only obtain a rough estimation of the IF when the frequency changes nonlinearly.If without the cross-term interference, the SPWVD can accurately extract the IF.Yet it can be seen in Figs.6 and 7 that the cross-term interference caused considerable estimation errors of the IF.However, both the SCT and the MSCT can obtain high-precision IF estimation.It is difficult to determine which one has a higher accuracy only from Figs.6 and 7.Therefore, we adopt some statistical indicators to evaluate the performance of the above four algorithms,including the maximum error,the Mean Absolute Error(MAE),and the Root Mean Square Error (RMSE).The comparison of these performance indicators is shown in Table 1.

As shown in Table 1,the performance of the STFT and WVD is far worse than that of the SCT and MSCT in the case of the non-linear frequency modulation signal.Especially for the SPWVD, the continuous occurrence of extremely large frequency extraction errors is unacceptable for practical applications.Though the maximum error of the MSCT is slightly higher than that of the SCT,the overall performance of the MSCT is better than that of the SCT.The RMSE term indicates that the IF extracted by the MSCT is the closest to the real IF.Therefore, the MSCT proposed in this paper can be applied to the IF estimation of a high spinning projectile in practical applications.

The above simulation tests are implemented by MATLAB Version 9.5.0 (R2018b) on a PC with Intel (R) Core (TM) i7-8550U CPU @1.80 GHz and 4 GB RAM.The calculation time of the STFT,the SPWVD, the SCT, and the MSCT are listed in Table 2.

Fig.5.The TFD of the signal (a) STFT, (b) SPWVD, (c) SCT, and (d) MSCT.

Fig.6.The extracted IF by four methods.

As shown in Table 2,the calculation time of the SCT and MSCT in one iteration is more than that of the STFT and SPWVD.Especially,the calculation time of the MSCT is 0.5100599 s and is the most among those four algorithms.Since the data sampling frequency is 100 Hz and the sampling time is 10 s, we can get 1000 sampling points.At each sampling moment, we will run the MSCT once to obtain the instantaneous frequency for the current moment.We have to run the MSCT 1000 times to obtain the instantaneous frequencies for every moment.After a simple calculation, the calculation time of the MSCT for each sampling point is 0.5100599 ms.At the same time, we set the window length of the MSCT as 512.In other words, the calculation time of the MSCT with a window length of 512 is 0.5100599 ms.If we decrease the window length,the calculation time will be less.To summarize, the real-time performance of the MSCT is not as good as that of the STFT, but it can meet the needs of most practical applications.

Table 1 The comparison of performance indicators.

Table 2 The comparison of calculation times.

5.Experimental validation

To further verify the performance of the proposed MSCT,a highspeed turntable experiment and spinning projectile flight tests were carried out.We consider the accuracy of the IF estimation to be the primary evaluation criterion.Thus, we introduce the reference value of IF through a high-precision gyroscope and a ballistic monitoring radar.Besides, we also adopt a three-axis magnetometer to collect the magnetic field intensity during the test.

5.1.High-speed turntable experiments

Fig.7.The errors of IF by four methods.

In this experiment, we adopted the Ellipse2-N equipment produced by SBG-France to collect the roll angular rate and magnetic field intensity.As shown in Fig.8, the Ellipse2-N were fixed in the center of the single-axis high-speed turntable by screws.The data sampling frequency was 200 Hz.The roll angular rate collected by the gyroscope and the magnetic field intensity collected by the magnetometer were both recorded.The rotation mode of the single-axis high-speed turntable was set through the console.

According to Section II, the magnetometer will output in harmonic waveform under the rotating condition.During the test,we set the turntable to rotate 100 rounds (36000) with a constant angular acceleration of 20/s.The rotation time was 85 s.After 42.5 s of constant acceleration, the rotation speed of the turntable reached 850/s.Then after a constant deceleration of 42.5 s, the turntable stopped rotating.During the test, the rotation frequency of the turntable was distributed between 0 Hz-2.36 Hz.The actual measured magnetic field intensity by the radial axis magnetometer is shown in Fig.9.

The rotation time of the turntable test was 85 s.Since the turntable accelerated from a standstill,we intercepted a section of the data with a higher rotation frequency for analysis and discarded the part with a lower rotation frequency.As can be seen in Fig.9(a),the acquired signal quality is relatively excellent.At the same time,the instantaneous frequency of the collected signal is piecewise linearly varying.The frequency band of the signal is narrow.Therefore, the accurate extraction of the instantaneous frequency has stricter requirements on the energy concentration of the timefrequency distribution.Similar to the numerical studies part, we used the STFT, the SPWVD, the SCT, and the MSCT to obtain the time-frequency distribution of the signal, then obtained the IF by extracting the ridge curve of the time-frequency distribution.

Specifically,for the MSCT,the window length is 512,the number of splines is 25, and the initial parameter matrix is 0.Since the instantaneous frequency of the above curve is distributed between 1.4 Hz and 2.4 Hz,the start frequency of the narrowband is set to be 0 Hz, the stop frequency is 3 Hz, and the number of spectrum subdivision points is 2000.For the STFT,the window size is also 512.The window sliding step size is 1.For the SPWVD, the time smoothing window size is 5000.The frequency smoothing window size is 512.The window function is Kaiser windows with shape factor β = 20.For the SCT, the window size is 512.The number of splines is 25.For the iteration time of the SCT and the MSCT,both of them are 2.The calculation process of the instantaneous frequency is the same as that of numerical studies.Firstly, we used the STFT,SPWVD,SCT,and MSCT to generate the TFDs.Then we adopted the TFD maxima-based method, as shown in Eq.(20), to extract the ridge curve of the TFDs.Finally,we could obtain the instantaneous frequency of the signal,i.e.,rotation frequency.The generated TFDs and calculated rotation frequencies are shown in Fig.10 and Fig.11.

Fig.8.High-speed turntable and measuring sensors.

As shown in Fig.10, all four methods can produce TFDs with a better concentration for linear frequency modulation signals.The SCT and the MSCT have the best TFD concentration.As shown in Fig.10, the IF obtained by the STFT is the worst matched with the real IF.The IF estimation accuracy of the SCT is comparable to that of the SPWVD.This phenomenon proves once again that the SPWVD can provide valuable time-frequency resolution for singlefrequency component signals with a high signal-to-noise ratio.Besides,the estimated IF of the MSCT almost coincides with the real IF,and the estimated IF curve is pretty smooth.Since the MSCT has a better frequency resolution, it can detect microscopic frequency changes.

To better evaluate the IF estimation accuracy of the above algorithms,we calculated the IF estimation errors and error statistical characteristics of each algorithm.The details are presented in Fig.12 and Table 3.

Fig.12 presents the IF estimation errors of four algorithms.All four methods can obtain the IF with acceptable accuracy.However,after comparing the error indexes in Table 3, it can be concluded that the proposed MSCT outperforms all the other considered algorithms in IF estimation.This result is due to the introduction of the CZT in MSCT.The CZT makes the time-frequency distribution generated by the MSCT have a higher concentration so that more accurate instantaneous frequency can be acquired.However, the MSCT can't deal with the sudden change of frequency.As shown in Fig.11,there is a sudden change in frequency near 10 s.According to the IF estimation errors presented in Fig.12, all methods can't estimate the IF accurately when the IF changes sharply.The authors are currently working on extending the MSCT to handle this kind of problem.

5.2.Low speed rotation flight test

In this section, we carried out a flight test to further prove the MSCT's superiority in IF estimation accuracy.The three-axis magnetometer and the three-axis gyroscope are strapped on a spinning projectile to measure the roll angular rate and the magnetic field intensity during the projectile flight.Then the abovementioned four methods were used to perform time-frequency analysis on the output from the radial magnetometer and estimate the roll angular rate of the projectile.Here we firstly demonstrate the relationship between the IF and the roll angular rate.The details can refer to Eq.(33).

where ωdenotes the roll angular rate, which can be obtained from the gyroscopes.Its unit is ?/s.

In the flight test, we used self-developed high dynamic measuring instrumentation to collect the angular velocity and the magnetic field strength.The high dynamic measuring instrumentation mainly consists of sensor unit, data acquisition unit, data storage unit, and data process unit.The high dynamic measuring instrumentation is shown in Fig.13.

The selected magnetic sensor is Honeywell's HMC1043,and the gyroscope is ADI's ADXR649.The specific parameters of the gyroscope and the magnetometers are shown in Table 4.

The projectile was launched from a barrel.A launch engine was installed at the tail of the projectile.Four fins were evenly distributed on the outside of the launch engine.Before launch,the wings were folded.After launch,the fins were unfolded.The unfolded fins could provide steering torque for the projectile during flight so that the projectile rotated quickly around its longitudinal axis.

Fig.9.(a) Radial magnetometer data (b) Rotation frequency.

Fig.10.The generated TFDs, (a) STFT, (b) SPWVD, (c) SCT, (d) MSCT.

The weather was clear with light breezes during the test.The lateral wind speed was about 5 m/s.The longitudinal wind speed was about 8 m/s,and the air humidity was 60%.The high dynamic measuring instrumentation was powered on after launch and was controlled by the overload switch.The data sampling frequency was 1 kHz.

The collected gyroscope data and the magnetometer data are shown in Fig.14 and Fig.15, respectively.

The STFT, SPWVD, SCT, and MSCT were used to analyze the collected radial magnetic signal during the projectile flight.Specifically, for the MSCT, the window length is 1024, the number of splines is 25, and the initial parameter matrix is 0.Since the instantaneous frequency of the above curve is distributed between 1.8 Hz and 2.2 Hz,the start frequency of the narrowband is set to be 1 Hz, the stop frequency is 3 Hz, and the number of spectrum subdivision points is 2000.For the STFT, the window size is also 1024.The window sliding step size is 1.For the SPWVD, the time smoothing window size is 5000,the frequency smoothing window size is 5000, the window function is Kaiser windows with shape factor β=20.For the SCT,the window size is 1024.The number of splines is 25.For the iteration time of SCT and MSCT,both of them are 2.The calculation process of the instantaneous frequency is the same as that of the turntable experiments.The generated TFDs and calculated rotation frequencies are presented in Fig.16.

Fig.11.The extracted IF, (a) STFT, (b) SPWVD, (c) SCT, (d) MSCT.

The analysis of narrowband signals requires a higher concentration of the TFD.As can be observed in Fig.16, the MSCT has the best TFD concentration and the highest IF estimation accuracy.In Fig.16 (a), the STFT can barely reveal the inherent time-frequency pattern of the narrowband nonlinear frequency modulation signal.As shown in Fig.16(b),the TFD generated by the SPWVD is too blur to reveal the IF trajectory.In Fig.16(c),the TFD generated by the SCT scatters the energy around the IF due to its coarse frequency resolution, and it is also inaccurate to estimate the IF.As shown in Fig.16 (d), the MSCT outperforms the STFT, the SPWVD,and the SCT as it reveals the true time-frequency pattern of the signal.To better evaluate the IF accuracy of these four algorithms,we calculated the estimation errors,as shown in Fig.17.Meanwhile,the statistical characteristics of errors are given in Table 6.

We found an interesting phenomenon from Fig.17 and Table 6 that the IF estimation accuracy of the SPWVD is worse than that of the STFT.At the same time,the signal of the magnetometer had a good quality.Therefore,we inferred that the magnetometer signal must contain other frequency components.To verify this doubt,we performed an FFT on the magnetometer signal.The spectrum of the above-mentioned signal is shown in Fig.18.

Table 3 The comparison of performance indicators.

Table 4 The specific parameters of sensors.

As can be seen from Fig.18, the main frequency is about 4 Hz.However, the real rotation frequency is distributed between 1 Hz and 2 Hz.There is a clear difference between these two frequencies.According to Ref.[32], when the nutation angle is large, the measurement signal contains not only the rotation frequency of the projectile but also the compound frequency of vibration and rotation.It is as follows.

wheredenotes the compound frequency,denotes the rotation frequency,denotes the vibration frequency.According to Eq.(34), the main frequency shown in Fig.18 means that the compound frequency of vibration and rotation is about 4 Hz.

Fig.18 shows that the magnetometer signal indeed contains multifrequency components.Other frequency components may be introduced by the vibration of the projectile.But it is difficult for the SPWVD to precisely differentiate the true IF trajectory from the spurious frequency contents introduced by the cross terms.What's more,the IF estimation accuracy of the MSCT drops sharply both in the beginning and end stages of the signal.There are a few possible reasons.Among them, the key is the MSCT has a windowing operation during the execution process, which is bound to affect the accuracy of the time-frequency analysis at the beginning and end of the signal.Therefore, follow-up work should be conducted on this issue.

Fig.12.The extracted IF errors.

According to Ref.[24], the SCT is focused on mono-component signals.For multi-component signals with multiple frequencies,the SCT cannot produce a time-frequency distribution with good concentration.Unfortunately,our proposed MSCT also suffers from this dilemma.However,we can set the frequency band range in the CZT to filter other frequency components, thereby realizing accurate analysis of specified frequency components.Based on this principle, we can accurately extract the rotation frequency in the above flight test.

5.3.High speed rotation flight test

Fig.13.High dynamic measuring instrumentation.

Fig.14.Roll angular rate measured by gyroscope.

In the above section,we carried out a low-speed rotation flight test to verify the instantaneous frequency estimation accuracy of our proposed algorithm.In the case of low rotation speed, the measurement results of the gyroscope are more reliable.In this section, we will conduct a high rotation speed test to verify the robustness of our algorithm.

Fig.15.Radial axis magnetic strength information during flight.

A three-axis magnetometer was strapped on the high spinning projectile to measure the magnetic field intensity during the projectile flight.Similarly,the magnetometer was installed in the radial axis of the projectile.Due to the low measurement accuracy of the large-range gyroscope, we chose to use high-precision ballistic monitoring radar to obtain the rotation frequency of the projectile.The sampling frequency of the magnetometer is 1 KHz, and the sampling time is 40 s.The collected magnetometer data is shown in Fig.19.

The obtained rotation frequency is shown in Fig.20.

We used the STFT, SPWVD, SCT, and our proposed MSCT to process the collected magnetometer data.Specifically, for the MSCT, the window length is 1024, the number of splines is 3, and the initial parameter matrix is 0.Since the instantaneous frequency of the above curve is distributed between 5 Hz and 22 Hz,the start frequency of the narrowband is set to be 4 Hz,the stop frequency is 24 Hz,and the number of spectrum subdivision points is 5000.For the STFT,the window size is also 1024.The window sliding step size is 1.For the SPWVD,the time smoothing window size is 5000,the frequency smoothing window size is 5000,the window function is Kaiser windows with shape factor β=21.For the SCT,the window size is 1024.The number of splines is 25.For the iteration numbers of SCT and MSCT, both of them are 2.After calculation, we can obtain the estimated instantaneous frequencies of four algorithms.They are shown in Fig.21.

Based on the rotation frequency obtained from the highprecision ballistic monitoring radar, we can acquire the estimation errors of four algorithms.They are shown in Fig.22.Meanwhile, the statistical characteristics of errors are given in Table 7.

As can be seen from Figs.21 and 22, our proposed MSCT can accurately estimate the rotation frequency using the collected radial magnetometer data.Compared with the other three algorithms, our algorithm has the highest frequency estimation accuracy.Based on Table 7, the mean average error of our proposed method is 0.0478 Hz,the mean average errors of other methods are all greater than 0.1 Hz.0.1 Hz rotation frequency estimation error means the angular rate measurement error of 36/s.Such a large error is unacceptable for strap-down attitude calculation and trajectory correction.The accuracy of our proposed algorithm can meet the above requirements.

Table 5 The specific parameters of projectile.

Table 6 The comparison of performance indicators.

Table 7 The comparison of performance indicators.

Fig.16.Four methods generated TFDs and estimated IFs for flight test.

In general, through the above-mentioned various simulation and experiments, it can be concluded that the algorithm we proposed can efficiently and accurately estimate the rotation frequency for the high spinning flying body.In practical applications,if there are extremely high requirements for real-time performance, this paper recommends the use of a simple STFT.If there are strict requirements for frequency estimation accuracy, the MSCT proposed in this paper will be a better choice.

Fig.17.Frequency estimation errors for flight test.

Fig.18.Spectrum of magnetic signal.

Fig.19.High speed rotation flight magnetometer output.

Fig.20.The rotation frequency of high spinning projectile.

Fig.21.The estimated rotation frequency of high spinning projectile.

Fig.22.The estimation errors of rotation frequency.

6.Conclusion

The high-speed rotation of the projectile brings great difficulties for the direct measurement of the roll angular rate.The radial magnetometers scheme is approved in this paper.However, the existing time-frequency analysis methods have inadequate abilities in processing the nonlinear frequency modulation signals.Extraction accuracy of instantaneous frequency by the existing methods cannot meet the requirement of guidance and control of the high spinning projectile.

To solve the rotation frequency estimation problem, this paper proposed a modified spline-kernelled chirplet transform (MSCT)algorithm.For the MSCT,the construction accuracy of the transform kernel is improved by replacing the STFT operator with the CZT operator.Since the construction accuracy of the transform kernel determines the concentration of time-frequency distribution, we can obtain a more precise instantaneous frequency by the MSCT.The performance of the MSCT was evaluated by a series of numerical simulations, high-speed turntable experiments, and real flight tests.Besides, the MSCT was compared with the STFT, the SPWVD, and the SCT for the estimation accuracy of rotation frequency.

The MSCT demonstrates noticeable superiority in timefrequency distribution concentration and instantaneous frequency estimation accuracy.Its extraction accuracy of the roll angular rate is comparable to that of the high-precision gyroscope.Above all, the MSCT has a wider range of the roll angular rate extraction than that of the high-precision gyroscope.Although the MSCT needs more computation and time consumption, it is not difficult to implement on a smart projectile with extreme functional hardware.In the future, our work will focus on the improvement of the real-time performance for the MSCT.

No conflict of interest exists in the submission of this manuscript, and this manuscript has been approved by all authors for publication.

The authors would like to acknowledge National Natural Science Foundation (NNSF) of China under Grant 61771059, National Natural Science Foundation (NNSF) of China under Grant 61471046,and Beijing Natural Science Foundation under Grant 4172022 to provide fund for conducting experiments.