Sinusoidal Frequency Modulation Fourier-Bessel Series for Multicomponent SFM Signal Estimation and Separation

Multicomponent sinusoidal frequency modulated (SFM) signals are widely used in radar, acoustics, and biomedicine. The instantaneous frequency (IF) characterizes important physical parameters of the real applications. In this paper, a sinusoidal frequencymodulation Fourier-Bessel (SFMFB) series is defined for IF estimation. It provides the signal decomposition on the Bessel function basis with a finer resolution, which proposes an extension of the performance and the applicability of the classic FourierBessel transform (FBT). Based on the property analysis of the SFMFB series, an algorithm of IF estimation and signal separation is introduced. Unlike the existing estimation methods which apply sliding windows to make an instantaneous approximation, the proposed method uses the global data, which provides a longer period gain, therefore achieving a better estimation performance. Moreover, considering that most estimation methods are invalid in multicomponent separation, the individual signals are well separated by the proposed algorithm, which facilitates the further monocomponent analysis. A performance comparison between the proposed method, the FBT, and another recently proposed sinusoidal frequency modulation Fourier transform (SFMFT) is also provided. Simulation results indicate that the proposed method outperforms the existing methods in estimation precision and computation load, and it is free of interference which exists in SFMFT.


Introduction
Multicomponent sinusoidal frequency modulated (SFM) signals have been widely used in various real-life applications such as radar target recognition, speech analysis, and biomedical disease diagnosis [1][2][3][4].One of the primary motivations for studying multicomponent SFM signals is from micro-Doppler parameter estimation of group targets.It is known that the frequency modulation induced by mechanical rotation or vibration of the group targets is time dependent and is multisinusoid overlapped in the time-frequency plane.As the spectral contents of the SFM signals vary with time, the frequency at a particular time is well described by the instantaneous frequency (IF).The IF characterizes important physical parameters and provides significant information for real applications.Therefore, it is desirable to have effective methods for IF estimation of the multicomponent SFM signals.Meanwhile, the individual signals separated from the original multicomponent signals will facilitate the further analysis of monocomponent if necessary.
For decades, IF estimation of the SFM signals has remained one of the most popular topics, and numerous IF estimation approaches have been proposed.Some of the most popular estimators are time-frequency distribution (TFD) based [5][6][7], polynomial based [8][9][10], cyclostationary based [11,12], or demodulation based [13].Considering that the SFM signal is one of the typical nonstationary signals, the IF estimation approaches based on the TFD have been widely used in the past few decades.Instead of constructing parametric models or equations, the TFD based methods analyze the time-frequency features of a signal first and the IF is estimated by tracking their periodical features [5][6][7].In these methods, the length of the time window is one of the most important parameters.To distinguish the time-varying frequency, the time window should be much less than the period of the SFM signals.At the same time, in order to obtain accurate IF, the window length should be enlarged to contain samples as many as possible.As the SFM signal can be categorized into the nonlinear polynomial frequency modulated (FM) signal, the IF estimation of the SFM signal can be realized based on the discrete polynomial-phase transform (DPT) [8].Combining with the Wigner-Ville distribution (WVD), the polynomial Wigner-Ville distribution (PWVD) [9] and its improvement [10] are also proposed as the IF estimators for the SFM signals.PWVD depicts the IF variance, and the frequency is estimated according to the peak of PWVD.To avoid the loss of the time-frequency resolution, the cyclostationary based methods, such as the autocorrelation functions [11,12], are introduced into IF estimation.Because of their simplicity and low computation burden, the cyclostationary based methods are widely used in IF estimation of the SFM signals.However, these methods are subjected to multiperiod errors because of the peak periodicity.In [13], a method for SFM signal demodulation is shown to be efficient in the IF estimation of the monocomponent SFM signal as well.
However, the above IF estimators are most only valid for the monocomponent signal.The correct IF estimation results can hardly be obtained when dealing with multicomponent SFM signals.For IF estimation of multicomponent, some literature presented that the IF can be extracted by projecting the phase term of the signals on some orthogonal function basis, such as the sinusoidal function basis, referred to as the sinusoidal frequency modulation Fourier transform (SFMFT) [14].However, interference terms occur in the SFMFT spectrum when there are more than two components, and the estimation precision of SFMFT is bounded by the sampling frequency.In nonstationary signal processing, another function basis widely used in signal decomposition is the Bessel function.The Bessel functions form the orthogonal basis and decay over the time, so that the signals which do not overlap in both the time and the frequency domain, including single frequency signals and linear frequency modulated (LFM) signals, can be represented well using the Fourier-Bessel transform (FBT) or the Fourier-Bessel (FB) series expansion [15][16][17][18].Because of the one-to-one relationship between the order of FB series and the frequency content of a signal, FB series expansion performs well in multicomponent signal separation.However, as the frequency content of the signal corresponds to a positive root of Bessel functions, the frequency resolution of Bessel function basis is fixed and imprecise, so that the estimation precision is bounded by the Bessel functions.Moreover, the damping property of Bessel functions causes that the projection of frequency content on the adjacent orders of a Bessel function basis is not zero, so the estimation error is enlarged periodically in some frequency interval.The periodic estimation errors cannot be ignored in the IF estimation when using Bessel based estimators because of its damping property.
Considering the existing problems of the IF estimation on the Bessel domain, this paper proposes an extension of the estimation performance as well as the applicability of the classic FBT and the FB series in IF estimation of multicomponent SFM signals.A sinusoidal frequency modulation Fourier-Bessel (SFMFB) series is defined, the properties of the SFMFB series are derived, and a novel IF estimation and signal separation approach for multicomponent SFM signals is introduced via the property analysis.The main contributions of this paper are in the following four aspects: (1) A metric -resolution is introduced into the Bessel function basis, by which the frequency resolution is refined.A SFMFB series is newly defined, and the properties of the SFMFB series are derived.
(2) The IF estimation error by SFMFB series is analyzed.An algorithm is proposed to reduce the periodic estimation error, which exists in the Bessel based estimators.The algorithm can be extended to use in other estimators whose function basis has the damping property as well.
(3) Some problems in IF estimation of the discrete signal are discussed, including the phase shift ambiguity revision and the maximal computation order of SFMFB series.
(4) A comparison of the IF estimation performance between the proposed method, the FB-based method, and the SFMFT-based method is provided.
The reminder of this paper is organized as follows.In Section 2, the definition of SFMFB series is given, and the properties related to the IF estimation are derived.In Section 3, the periodic estimation error in the Bessel based estimators is analyzed, and an SFMFB expansion based algorithm of multicomponent SFM signal IF estimation and signal separation is proposed.In Section 4, simulation results are discussed.Section 5 concludes the paper.[15].With the metric -resolution introduced into the Bessel functions, over the finite time interval (0, ), the phase term of a continuous signal () can be expressed as a weighted sum of a finite number of Bessel functions:

Definitions and Basic Properties
where   is the th SFMFB coefficient,   (⋅) is the th Bessel function,   is the th positive root of the th Bessel function in ascending order, and  is the resolution of Bessel function basis.The SFMFB series with the -resolution is called the -SFMFB series for simplicity.The th positive root   of the th Bessel function can be computed by the Newton-Raphson method [19], which is added in Appendix A of the paper.Similar to the calculation of FB coefficients, considering the phase term of the signal (), the th -SFMFB coefficient   of the signal over the interval (0, ) is computed as the projection on the Bessel functions with -resolution   (  /): where (⋅) is the Dirac delta function.Only when  = , the inner product equals zero.The orthogonality of the function basis is proved as follows.
Proof.Note the equation in the handbook [20] Let  = /; then It is proved that the function basis   (  /) is orthogonal with   (  /) with the power function .

Quasi-Periodicity.
The quasi-periodicity of the th positive root   of the th Bessel function presents as where  = 0, 1, . ... Based on Bessel function's property in [21], when  → ∞, the difference between the two consecutive positive roots tends to : When  ≥ 7, the two consecutive roots satisfy |( +1 −   ) − | < 10 −3 .Namely, the difference between the two consecutive positive roots of   () = 0 is extremely close to , where the bias is smaller than 10 −3 .It should also be emphasized that the bias approaches zero rapidly as the order  ascends.Moreover, the quasi-periodicity can help to compute the  + th positive root of th Bessel function by the th one.

Magnitude Property.
The relationship between the IF, -resolution, and the order of SFMFB series is reflected on the magnitude of the series, which is called the magnitude property.Assume that a SFM signal presents as where  is the amplitude,  is the IF, and  and  are modulation indexes of the sine term and the cosine term, respectively.The relationship between the IF, -resolution, and the SFMFB coefficient presents as follows: where  max = arg max  {|Re(  )|} is the order of the SFMFB coefficient with the maximal magnitude of the real part.
Proof.Substitute the SFM signal () in ( 10) into the SFMFB series in (2) and consider that  = 2, so the SFMFB series of signal () is given as where Based on the indefinite integral with the Bessel function in [20] ∫   () d =  +1 () , 1 in ( 14) is derived as Considering another indefinite integral with the Bessel function in [20] ∫  +1 () d = −  () , 1 in ( 17) is further derived as Then  1 is simplified as In the same way,  2 is obtained as Consider the equation in [18] Bring the parameter  into (22); the real part and the imaginary part of ( 22) can be obtained as Substitute ( 23) into ( 20), and substitute (24) into (21);  1 and  2 are, respectively, derived as Substitute (25) into (12); the SFMFB series of the signal () is obtained as where  0 is imaginary, tan   = ( 2  −  2  2  2 )/, and tan  = /.Since the frequency  and the time interval  are both constant, the magnitude of the SFMFB coefficient varies with the order .When  →   /, the absolute value of the real part of the SFMFB series |Re(  )| presents as a maximal coefficient.Therefore, the magnitude property in (11) is proved.

Frequency Resolution.
We suppose that the operator ≺ is slightly less than.Considering (8) in quasi-periodicity and (11) in magnitude property, the difference of the two consecutive theoretical frequencies of the -SFMFB series is calculated as In the analysis of magnitude property, the relationship between the IF with the -resolution, the signal interval, and the order of -SFMFB series is derived.As (11) stated, when   → 2, the magnitude of the real part of the -SFMFB coefficient reaches the maximum, where  corresponds to the IF of the signal.Thus, the IF can be estimated according to the maximal -SFMFB coefficient where  max = arg max{|Re(  )|}.Suppose the estimation error equals the absolute value of the difference between the estimated IF and the real IF: where   max is the estimated IF and  is the real IF.Thereby, the maximal estimation error of the -SFMFB series is bounded by the frequency resolution Δ  , which satisfies max   ≺ 1/2.With  increasing, the frequency resolution is refined, and the estimation error can be reduced.However, only when the signal interval is large enough will the resolution capability be improved.

Parameter Estimation and Separation Algorithm
As stated in Section 2, the IF estimation can be obtained by the SFMFB coefficient with the maximal real part.In this section, the estimation error is analyzed in detail.Based on the analysis, an approach is introduced to reduce the estimation error and to separate multicomponent signals.

Central Estimation Error.
In Fourier transform, the estimation error is bounded by half of the bin size.Therefore, the maximal estimation error equals half of the frequency resolution (or the bin size) of two adjacent function bases, when the real frequency equals the mean value of the two sinusoidal basis.In this case, the frequency content only projects on a single order of series or a single spectral line of a spectrum, such as the FFT spectrum.Different from the FFT spectrum, in SFMFB series, a single frequency content projects on several adjacent orders.In other words, the theoretical frequency   reaches the maximum at the order , while the projections of the frequency   on several adjacent orders are not zero.
In SFMFB series,   /2 and  +1 /2 correspond to two adjacent orders  and  + 1, respectively.We suppose that the absolute value of the real part of the th SFMFB coefficient   is equivalent to that of the  + 1th coefficient With some approximation, (30) is derived as (see Appendix A) Based on the basic inequality, the real IF  satisfies where   ̸ =  +1 .When |Re(  )| equals |Re( +1 )|, the real IF will be more than the average of the theoretical frequencies of the two adjacent orders.Namely, the IF estimated by the SFMFB coefficient of the maximal real part is smaller than the real IF.On the one hand, the estimation error is bounded by the frequency resolution; on the other hand, unlike the Fourier transform, whose maximal estimation error is bounded by half of the frequency resolution, the maximal estimation error of SFMFB series is larger than half of the frequency resolution: Therefore, a solution is introduced to reduce the maximal estimation error below Δ  /2.An effective way is to calculate the center estimation error  cen and take the revised estimated IF as the difference between the original estimated IF and  cen .The center estimation error is defined as the mean value of the maximal estimation error and the minimal estimation error.Because of the quasi-periodicity of the Bessel function basis,  cen varies slightly with different order of series.The estimated IF revision process is summarized in Algorithm 1.

Discussion
. The phase shift ambiguity occurs when the total phase shift of the discrete signal is beyond 2.Considering that the phase shift ambiguity exists in discrete signal sampling generally, the phase shift ambiguity is discussed and revised in this subsection.Moreover, to avoid the superfluous computation amount in the calculation of the SFMFB series, the maximal order of the SFMFB coefficients of the discrete signal is also discussed.

Phase Shift Ambiguity.
As the phase measurements are only possible in (−, ), when the modulation index and the signal sequence are large, the phase shift ambiguity occurs when the total phase shift is more than 2 when sampling, which may lead to wrong estimation result of the IF.Assume that the total phase shift is the integer times of 2, where the integer is called the number of phase shift ambiguities.The phase shift ambiguity can be revised by comparing the adjacent samples.A solution of the phase revising function is where pha() = Im[ln ()] and phar() is the difference between the phase shift and the measurements.However, the phase can be revised when the maximal phase shift between two adjacent samples is smaller than .

Maximal Order of Series.
For discrete signal, the Nyquist sampling theory limits that the maximal frequency presented by sampling is half of the sampling frequency.Meanwhile, as (11) states, there is a one-to-one relationship between the frequency content and the order of series.Therefore, the maximal order  of SFMFB series is limited by the Nyquist sampling theory as well.The maximal order depends on two factors, that is, the length of sampling sequence and the -resolution, which is derived as (see Appendix B) Specially, the maximal order is  ≈  approximately when the -SFMFB series take the zeroth Bessel functions as the function basis.

The Proposed Algorithm.
We adopt the CLEAN algorithm in the signal separation.We suppose that there does not exist any signal component when the energy is below a threshold .Since there is a linear relationship between the frequency content of a signal and the order of the series, the energy of a signal can be represented by the square of its corresponding order of coefficients.The IF estimation and the multicomponent separation algorithm of SFM signals is illustrated in Algorithm 2.

Simulation Results
In this section, the estimation performance of the proposed method is simulated and discussed.Moreover, the estimation performance of the proposed method, the SFMFT-based method, and the FB-based IF estimation method is also compared.

Comparison of the Estimation Precision before and after Revision.
The simulation is conducted to analyze and compare the estimation precision before and after revision with different -resolution in this part.Considering the quasi-periodicity, the IF in a quasi-period range is simulated to illustrate the estimation performance.The IF of the SFM signal is discretized into a series  () ∈ [9.8753, 10.3753] Hz with the step Δ () = 0.02 Hz between each IF content.A metric to evaluate the estimation performance is the absolute estimation error  ()   = | f() −  () |, which equals the absolute difference between the estimated IF f() and the real IF  () .The estimation error  ()   before and after revision with different -resolution is shown in Figure 1.
By comparing Figures 1(a) and 1(b), it shows that the estimation precision is improved distinctly by the revision algorithm.In both Figures 1(a) and 1(b), the estimation error takes positive proportion to the -resolution.The more the  values, the smaller the estimation error.The estimation precision is improved with the -resolution increasing when the signal length is large enough, which coincides with the analysis in ( 27).The -resolution decides the resolution unit.However, the resolution capability depends on both the signal length and the -resolution.

Estimation Performance with Different Window Size.
To present the effect of the window size  for the estimation performance, the estimation error before and after revision with different window size and the same -resolution is simulated.The simulation results are shown in Figure 2.
As shown in Figure 2, with the same -resolution, the estimation error takes inverse proportion to the window size  of the signal, which agrees with the conclusion in [22] as well.

Robustness under Different SNR Conditions.
In order to present the robustness of the IF estimation approach under different signal-to-noise ratio (SNR) conditions, white Gaussian noise is added to the simulations.200 times of Monte Carlo simulations are conducted under SNR = 8 dB to −4 dB with different -resolution.Simulation results are shown in Figure 3.
From Figures 3(a) to 3(d) we can grasp that the estimation error is smaller than 0.1 Hz when SNR > 0 dB.Moreover, with the -resolution increasing under the same SNR condition, the estimation precision as well as the robustness is improved.However, the quasi-periodicity of the estimation error is contaminated and the performance decays distinctly when SNR = −4 dB.Because the magnitude of the several SFMFB coefficients adjacent to the maximal one is relatively large, the maximal coefficient is submerged by the adjacent coefficients in strong noise.

Estimation Performance with Different Modulation
Index.As (26) denotes, with the modulation index increasing, the maximal SFMFB coefficient is increasingly distinct.Hence, the more the modulation index is, the more robust the estimation performance is.To analyze the estimation performance with different modulation index, 200 times of Monte Carlo simulations are conducted under different SNR conditions.One metric to evaluate the estimation performance is the normalized root-mean-square error (NRMSE).NRMSE values the normalized Euclidean distance between Mathematical Problems in Engineering the estimated IF f and the real one , where the IF is  = 4 Hz in the simulation.NRMSE is defined as where ‖ ⋅ ‖ denotes the  2 -norm operator.The NRMSE curves under the conditions where SNR = 8 dB to −4 dB are shown in Figure 4.As shown in Figure 4, the modulation index of the signal affects the robustness of the estimation approach based on the SFMFB series.With the modulation index increasing in a certain range, the NRMSE of the estimated IF decreases, and the estimation performance and the robustness of the proposed approach improve under the same SNR condition.Nevertheless, the robustness cannot be ideal for the signal whose modulation index is small.

Comparison of the Methods.
In this subsection, the proposed method and the SFMFT-based method are simulated on the same signal model.Additionally, the estimation performance of the IF by the proposed method is compared with the FB-based method in single frequency estimation.

The Proposed Method.
The multicomponent SFM signals are modeled as (10), which consist of three components.The signal length is  = 256; the sampling length is   = 256 Hz and  = 10.The amplitude   , the modulation index   , and the IF   of the three components are set as follows, respectively: signal 1: From Figure 6 we can see that the TFD of the multicomponent SFM signal can hardly illustrate the IF of the signal.For the separated individual components, except for a little instability on the edges of the TF images, the TFD of an individual component can reflect the signal well.To guarantee the comparability of estimation performance of these two estimators, the FB-based and the proposed method is simulated on the signal with the same length  = 256.
The estimation errors of these two conditions are shown in Figure 7.
In a period of the frequency content, the maximal estimation error of the FB-based estimator is about 0.5 Hz, which is much larger than that of the proposed method.However, the estimation error decreases distinctly when the -resolution of the proposed method ascends.

The SFMFT-Based Method.
The essential of the proposed method and the SFMFT-based method [14] is to decompose the IF of the signals on some function basis, while the former one is the Bessel functions, and the latter one is the sinusoidal functions.The SFMFT is a parametric approach for the SFM signal estimation proposed in recent years.It provides the spectrum of multicomponent SFM signals.In this part, we compare estimation performance of the SFMFTbased method with the proposed method.As stated in [14], the key parameter that affects the estimation performance of SFMFT is the signals length .Hence, the IF of the SFM signals with different length (i.e.,  = 128,  = 256,  = 512, and  = 1024) are estimated by the SFMFT, respectively.The spectra of the signals with different length calculated by SFMFT are shown in Figure 8.
From Figure 8, we can grasp that when  = 128 the bin size of the spectrum is 1 Hz, when  = 256, the bin size is refined to 0.5 Hz, and when  = 1024 the bin size is further refined to 0.125 Hz.The estimation precision becomes finer as the signal length increases.However, Figure 8(d) shows that except the 5.25 Hz, 19.63 Hz, and 34.38 Hz, there are many interference terms in the spectra.
The frequency resolution of the function basis is an important metric that affects the estimation performance.When the SFMFT and the proposed method are of the same frequency resolution, the computation time of 1000 times simulation of these two methods is shown in Figure 9.
As shown in Figure 9, in conditions where the frequency resolution is less than 0.25 Hz, the proposed method obtains the less computation amount, in contrast with the SFMFTbased method.When finer estimation precision is required, the computation amount of the proposed method is much less than the SFMFT-based method.However, the computation amount of the proposed method increases exponentially with -resolution arising, so the SFMFB series with large resolution may not meet the requirement of real-time signal processing conditions.
From the theoretical analysis and the above simulations, the major subjects that contribute to the IF estimation performance of the proposed method, the FB-based method, and the SFMFT-based method are shown in Table 2.

Conclusion and Discussion
In this paper, an approach to the IF estimation of the multicomponent SFM signal and signal separation is developed.A newly defined SFMFB series is introduced to illustrate the complex modulation directly, which provides an accurate IF estimation of the multicomponent SFM signals.Based on the property analysis of the series and its function basis, an algorithm is proposed to reduce the periodic estimation error, which exists in the Bessel based estimators.The algorithm can be extended to use in other estimators whose function basis has the damping property as well.In addition, the phase shift ambiguity revision and the maximal computation order of SFMFB series in IF estimation of the discrete signal are considered.
The advantages and the limitations of the modifications in SFMFB series are listed as follows.The advantages are that (1) the metric -resolution reduces the frequency resolution of the Bessel function basis less than 1/2 Hz.In conventional Fourier-Bessel series expansion, the frequency resolution of the Bessel function basis (1/2 Hz) is fixed and imprecise.By introducing the metric -resolution, the frequency resolution is refined as 1/2 Hz. (2) The estimation revision reduces  From Figure 10 we can see that when  > 8, the ratio is less than 1.2, and when  → ∞, the ratio approaches 1.So when  → ∞,   ≈ −1.Hence, (B.

1 .
The Continuous Signal.The conventional FB series decomposes a signal into a set of Bessel functions   (  /)

Figure 1 :Figure 2 :
Figure 1: Estimation error before and after revision with different -resolution: (a) before revision and (b) after revision.
and  1 = 5.19 Hz; signal 2:  2 = 3,  2 =  2 = 1, and  2 = 19.64Hz; signal 3:  3 = 2,  3 =  3 = 1, and  3 = 34.38Hz.The SFMFB coefficients are calculated and the SFMFB spectrum is shown in Figure 5.In Figure 5, there are three distinct IF components.The estimated IFs and absolute errors of the multicomponent signals are listed in Table 1.The absolute estimation errors are within the maximal estimation error when  = 10.The multicomponent SFM signal is shown in Figure 6(a), and the monocomponent signals after separation are shown in Figures 6(b)-6(d), respectively.

Figure 10 :
Figure 10: The ratio   varying with order .
According to the definition of the -SFMFB series, the inner product of the th and the th function basis   (  /) and   (  /) is computed as ) Obtain the revised estimated IF f =   max −  cen .Output: The revised estimated IF f.While ∑  (   ) 2 / ∑  =1 (  ) 2 > .Output: The number of signal components , estimated IFs f and separated individual components   () (or   ()), where  = 1, . . ., .

Table 1 :
Estimated results of the proposed method.The classic FBT and the FB series is an IF estimator of the single frequency signal, but the frequency resolution of FB-based method is fixed and imprecise.In SFMFB series, the frequency resolution is dependent on both the signal length  and the -resolution.

Table 2 :
Major subject of the three IF estimation methods.

The Calculation of the Maximal Order of Series in Section 3.2
5)is obtained approximately as1 √ ( 2  −  2  2  2 ) 2 +  2  2  2  2  2 ) 2 +  2  2  2 .Suppose the sampling frequency is   , the maximal order corresponding to half of the sampling frequency is , and the sampling length is .Therefore,   can be calculated by the theoretical frequency of the first order  1 and the frequency resolution Δ  ,  =  1 + ( − 1) Δ  =  1    1 ≈ 2.4048 when taking  1 as the first root of the zeroth Bessel function. Hnce, the maximal order is derived approximately as  ≈ .